Labs/Machine Learning/linear-regression

Machine LearningBeginner

~20 min

Linear Regression

Fit hyperplanes to multidimensional feature spaces by minimizing quadratic prediction residuals via Ordinary Least Squares and gradient descent.

#Ordinary Least Squares (OLS)#Mean Squared Error (MSE)#Normal Equation#Gradient Descent Optimization#Multicollinearity#Outlier Sensitivity

Overview

The Problem Statement

Modern engineering systems continuously predict continuous real-valued targets—server latency, disk IOPS utilization, financial risk, or cloud compute costs—from observable feature telemetry. Linear regression establishes the optimal affine transformation mapping feature vectors to target values by minimizing scalar residual error.

System Invariant:The optimal Ordinary Least Squares parameter vector w* ensures that the residual error vector e = y - ŷ is strictly orthogonal to the column space of the feature design matrix X (i.e., Xᵀ(y - Xw*) = 0).

✓ When To Use

Predicting continuous metrics where feature effects are approximately additive and monotonic.
High-interpretability domains where each coefficient represents an audit-traceable physical marginal cost.
Ultra-low-latency production scoring (evaluating a dot product wᵀx + b takes < 15 nanoseconds).
Baseline modeling to establish empirical performance lower bounds before deploying deep neural networks.

✕ When NOT To Use

Complex non-linear relationships with high-order feature interactions without explicit polynomial feature transforms.
Classification problems with bounded categorical outcomes (use Logistic Regression or Softmax).
High-dimensional datasets with severe multicollinearity without L1/L2 regularization (Ridge/Lasso).
Datasets containing heavy-tailed Cauchy noise or extreme leverage outliers without Huber / robust loss formulation.

Why It Exists

Catastrophic Outage Scenario

An automated cloud autoscaler uses unregularized linear regression on CPU and request rates to scale container replicas. A single batch scraper injects an extreme leverage outlier (low CPU, massive socket connection count). Because Ordinary Least Squares minimizes squared residuals, the fitted hyperplane pivots violently toward the outlier. The autoscaler predicts zero resource demand during peak shopping hours, shutting down 80% of cluster nodes and inducing a complete site outage.

Downstream System Degradation:

Quadratic penalty amplification causes extreme vulnerability to data corruptions and anomalies.
Singular matrix inversion failure when collinear features create zero eigenvalues in XᵀX.
Sub-optimal capacity forecasting causing thrashing and severe autoscaling latency cascades.
Silent model drift when underlying system relationships shift from linear to saturated non-linear regimes.

How It Works

Linear regression models the relationship between dependent scalar y and d-dimensional feature vector x as: ŷ = wᵀx + b = w₁x₁ + w₂x₂ + ... + w_d x_d + b.

The optimization objective is the Mean Squared Error (MSE) loss function: J(w, b) = (1 / 2n) ∑ (yᵢ - (wᵀxᵢ + b))². Because J is a convex quadratic paraboloid, it possesses a unique global minimum with zero local minima traps.

The parameter vector can be solved either analytically in closed form via the Normal Equation: w* = (XᵀX)⁻¹ Xᵀy, or iteratively via Gradient Descent: w ← w - η ∇_w J.

Single-Node vs Distributed Reality:

Small to medium datasets (n < 100,000, d < 1,000) fit comfortably in server RAM and are solved in sub-second time via Cholesky or SVD decomposition of the Normal Equation. Massive streaming datasets (billions of events) cannot afford the O(d³) matrix inversion cost and must use Stochastic Gradient Descent (SGD) or distributed parameter servers (AllReduce).

Normal Equation (Closed-Form Analytical Solution)

Directly computes w* = (XᵀX)⁻¹ Xᵀy by setting the gradient of MSE to zero. Requires no learning rate tuning or iterations.

+ Exact mathematical solution in a single step; no hyperparameters to tune.

- Matrix inversion scales with O(d³); fails if XᵀX is singular (collinear features) or d > 10,000.

Batch Gradient Descent (Iterative)

Iteratively steps parameters down the loss surface using the full dataset gradient: w ← w - η (1/n) Xᵀ(Xw - y).

+ Scales efficiently to millions of features d; memory footprint is O(nd).

- Requires computing gradients over the entire dataset per step; sensitive to learning rate η.

Mini-Batch / Stochastic Gradient Descent (SGD)

Approximates the gradient over small randomized batches (e.g. 64 or 256 samples), allowing continuous online parameter updates.

+ Fast convergence; constant memory footprint O(batch_size · d); handles streaming real-time telemetry.

- Noisy gradient updates cause parameters to oscillate around the minimum rather than landing precisely at zero error.

Ridge Regression (L2 Regularization / Tikhonov)

Adds quadratic parameter penalty λ ||w||² to loss: J_ridge = MSE + λ ∑ wᵢ². Solved via (XᵀX + λI)⁻¹ Xᵀy.

+ Guarantees (XᵀX + λI) is strictly invertible; suppresses parameter explosion caused by multicollinearity.

- Introduces regularization hyperparameter λ; shrinks coefficients toward zero without performing feature selection.

Engineering Prediction Challenge

Prove It: Predict the System Behavior

Reason before simulating

System State Given

A dataset has 10 points generated by y = 2x + 1. An extreme outlier is introduced at (x=10, y=0).

Prediction Question:

When minimizing Mean Squared Error (MSE) to fit the line y = wx + b, how does this single outlier affect the model parameters?

The outlier is automatically ignored because 10 points easily outvote 1 point.

The regression line is aggressively pulled downward toward the outlier, destroying accuracy across the 10 legitimate points.

The slope w increases to compensate for the disturbance.

Interactive Visualizer

Interact directly with this distributed system primitive. Experiment with fault injection and observe state transitions.

Linear Regression & Residual Minimization

Model: ŷ = w·x + b. Drag points or tweak learning rate to inspect convergence.

MSE Loss18.7437

Weight (Slope w)0.100

Bias (Intercept b)0.500

Iterations0

Data Points8 (Draggable)

Data Point (Draggable)

Fitted Line ŷ = wx + b

Residual Error (y - ŷ)

Learning Rate (α):0.03

Dataset Noise (σ):0.5

Failure & Experiment Scenarios:

Observe: Increasing α beyond the curvature boundary causes the residual error to oscillate wildly instead of minimizing.

Build It Step-by-Step

The engineering blueprint. Expand each sequential milestone to examine requirements, architectural impacts, common pitfalls, and blueprint pseudocode.

Formulate Design Matrix & Target Vector

To unify slope weights and bias into a single vector dot product, augment the n×d feature matrix with a column of ones: X_aug = [1, X], yielding parameters θ = [b, w₁, ..., w_d]ᵀ.

def prepare_design_matrix(X):
  n = X.shape[0]
  return np.hstack([np.ones((n, 1)), X])

Implement Analytical Normal Equation Solver

Directly computing (XᵀX)⁻¹ via naive matrix inversion is numerically unstable. Production libraries solve the least squares problem via QR decomposition or SVD (e.g. scipy.linalg.lstsq) to handle rank-deficient matrices.

def solve_ols(X, y):
  # Using Moore-Penrose pseudo-inverse for numerical stability
  return np.linalg.pinv(X.T @ X) @ X.T @ y

Implement Vectorized Batch Gradient Descent

Vectorized gradient computation eliminates Python loops: residuals = Xw - y, gradient = (1/n) Xᵀ · residuals, update: w ← w - η · gradient.

for epoch in range(max_epochs):
  residuals = (X @ w) - y
  grad = (1 / n) * (X.T @ residuals)
  w -= learning_rate * grad
  if np.linalg.norm(grad) < tolerance:
    break

Add Metric Telemetry: R² Score and RMSE

Root Mean Squared Error (RMSE) provides error in the target unit (e.g. milliseconds). The Coefficient of Determination (R²) measures the fraction of variance explained by the model: R² = 1 - SS_res / SS_tot.

def compute_r2(y_true, y_pred):
  ss_res = np.sum((y_true - y_pred) ** 2)
  ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
  return 1.0 - (ss_res / (ss_tot + 1e-9))

Code Implementation

Production-ready, idiomatic reference implementations across Go, TypeScript, Python, and Java with sticky language switching, copy-to-clipboard, and source download.

Core Architectural Design Decisions

Invariants & Mechanics

Decision 01

Used np.linalg.pinv instead of naive np.linalg.inv(X.T @ X) to prevent LinAlgError crashes on collinear matrices.

Decision 02

Kept weights and bias separate in interface while unifying them internally in design matrix.

Algorithmic Complexity:SVD Fit Time: O(n·d² + d³). GD Fit Time: O(epochs · n · d). Inference Time: O(d) sub-microsecond dot product.

linear_regression.py•65 LOC

import numpy as np
from typing import Tuple, Optional

class LinearRegression:
    """
    Production-grade Ordinary Least Squares Linear Regression.
    Supports both analytical Normal Equation (SVD) and Gradient Descent solvers.
    """

    def __init__(self, solver: str = "svd", learning_rate: float = 0.01, max_iter: int = 1000):
        self.solver = solver
        self.learning_rate = learning_rate
        self.max_iter = max_iter
        self.weights: Optional[np.ndarray] = None
        self.bias: float = 0.0

    def fit(self, X: np.ndarray, y: np.ndarray) -> "LinearRegression":
        """
        Fit linear model to training data (X: [n_samples, n_features], y: [n_samples]).
        """
        X = np.asarray(X, dtype=np.float64)
        y = np.asarray(y, dtype=np.float64).reshape(-1, 1)
        n_samples, n_features = X.shape

        if self.solver == "svd":
            # Add bias column of ones
            X_aug = np.hstack([np.ones((n_samples, 1)), X])
            # Moore-Penrose pseudo-inverse handles rank-deficient systems
            theta = np.linalg.pinv(X_aug) @ y
            self.bias = float(theta[0, 0])
            self.weights = theta[1:, 0]
        elif self.solver == "gradient_descent":
            self.weights = np.zeros(n_features, dtype=np.float64)
            self.bias = 0.0

            for _ in range(self.max_iter):
                # Forward prediction: y_hat = Xw + b
                y_pred = (X @ self.weights.reshape(-1, 1)) + self.bias
                error = y_pred - y  # [n_samples, 1]

                # Gradient computations
                grad_w = (1.0 / n_samples) * (X.T @ error).ravel()
                grad_b = float((1.0 / n_samples) * np.sum(error))

                # Gradient descent step
                self.weights -= self.learning_rate * grad_w
                self.bias -= self.learning_rate * grad_b
        else:
            raise ValueError(f"Unknown solver: {self.solver}")

        return self

    def predict(self, X: np.ndarray) -> np.ndarray:
        """Score unseen feature vectors: ŷ = Xw + b."""
        if self.weights is None:
            raise RuntimeError("Model must be fitted before calling predict().")
        X = np.asarray(X, dtype=np.float64)
        return (X @ self.weights.reshape(-1, 1)).ravel() + self.bias

    def score(self, X: np.ndarray, y: np.ndarray) -> float:
        """Returns R² coefficient of determination."""
        y_pred = self.predict(X)
        ss_res = np.sum((y - y_pred) ** 2)
        ss_tot = np.sum((y - np.mean(y)) ** 2)
        return float(1.0 - (ss_res / (ss_tot + 1e-12)))

Implementation Strategy: Vectorized Python implementation leveraging NumPy BLAS/LAPACK routines. The SVD solver uses Moore-Penrose pseudo-inversion for unconditional numerical stability against rank-deficient collinear data.

SVD Fit Time: O(n·d² + d³). GD Fit Time: O(epochs · n · d). Inference Time: O(d) sub-microsecond dot product.

Edge Cases & Failure Modes

Investigate real-world production incidents. Expand each failure mode to analyze symptoms, root causes, and defensive engineering mitigations.

Scenario: Multicollinear Feature Columns (Rank Deficiency)

Consequence: Two feature columns are linearly dependent (e.g. latency_ms and latency_seconds). Matrix XᵀX has zero determinant; naive inversion throws singular matrix exception.

Engineering Solution: Use SVD pseudo-inverse (pinv), remove collinear features via Variance Inflation Factor (VIF) filtering, or add L2 Ridge regularization (XᵀX + λI) which guarantees strictly positive eigenvalues.

Scenario: Extreme Leverage Outlier

Consequence: Because MSE squares the error e², an outlier with residual 100 has 10,000x the weight of a point with residual 1, pivoting the hyperplane and ruining legitimate predictions.

Engineering Solution: Adopt Huber Loss (smooth transition from quadratic to linear penalty for |e| > δ) or replace with RANSAC / Quantile Regression.

Scenario: Zero Feature Variance

Consequence: A feature column has identical constant values across all rows. Gradient descent produces zero gradients; analytical solver divides by zero.

Engineering Solution: Implement variance threshold preprocessing to drop zero-variance features prior to model ingestion.

Production Considerations

Production checklist covering observability metrics, horizontal sharding, reliability fail-open patterns, and client identity security.

Scaling to 10x / 100x Traffic:

For massive batch datasets (100M+ rows), avoid computing (XᵀX) on single nodes; stream through MapReduce or Spark MLlib using BlockMatrix multiplication.
Export fitted weights [w, b] to SIMD-accelerated C++ / Rust binaries or ONNX runtime for sub-10ns scoring directly inside API gateway proxies.

Telemetry & Observability:

PROMETHEUS METRICS:

➔ model_prediction_latency_nanoseconds (gauge)

➔ model_residual_mse_rolling (gauge)

➔ model_feature_drift_kl_divergence (histogram)

DISTRIBUTED TRACES & LOGS:

➔ OpenTelemetry span for 'model.infer' capturing input feature vector hash and execution time.

Trade-offs

Select an algorithm below to dynamically compare memory footprints, burst tolerance, boundary precision, and distributed suitability.

Approach	Advantages	Disadvantages	Use When
Ordinary Least Squares (OLS)	Zero hyperparameters; closed-form mathematical guarantees; fast exact solution.	O(d³) matrix inversion scaling; hyper-sensitive to outliers.	Clean datasets with d < 1,000 features where interpretability is paramount.
Ridge Regression (L2)	Guaranteed invertible; stabilizes parameters under multicollinear features.	Requires tuning regularization strength λ.	Correlated telemetry metrics (e.g. disk read IOPS and disk write IOPS).
Lasso Regression (L1)	Drives non-informative feature weights to exact zero; performs automatic feature selection.	Non-differentiable at w=0; requires coordinate descent solver.	Sparse high-dimensional problems with thousands of noisy telemetry signals.