A mathematical and architectural derivation of Simple, Multiple, and Polynomial Linear Regression, Ordinary Least Squares (OLS) Normal Equations, Gradient Descent dynamics, Ridge (L_2), Lasso (L_1), Elastic Net penalty topologies, and rigorous metric evaluation.
Linear Regression models the relationship between a continuous scalar dependent variable y in R and one or more independent predictor features x in R^d assuming a linear parametric structure with additive stochastic error.
For a single predictor x, the regression model is parameterized by slope w_1 and intercept w_0:
Minimizing the sum of squared residuals yields the closed-form Ordinary Least Squares (OLS) solutions:
When extended to d input features, each observation x_i = [xi1, xi2, dots, xid]^T is augmented with a dummy bias feature xi0 = 1:
Where X in RN × (d+1) is the design matrix, w in R(d+1) is the parameter vector, and y in R^N is the target response vector.
The objective is to minimize the Residual Sum of Squares (RSS):
J(w) = (1 / 2) Xw - y₂² = (1 / 2) (Xw - y)^T (Xw - y)w:
∇w J(w) = X^T (Xw - y) = X^TXw - X^Ty∇w J(w) = 0) establishes the Normal Equation:
N-dimensional sample space. Generally does not lie in the column space of X.col(X) to y. Residual vector e = y - y_predicted is perpendicular to col(X) (X^T e = 0).y → y_predicted = Hy. Symmetric (H=H^T) and idempotent (H^2 = H).Real-world relationships frequently exhibit non-linear curvature (e.g., accelerating aerodynamic drag, diminishing marginal returns). Polynomial Regression captures non-linear trends by projecting input features into a higher-dimensional polynomial feature space:
In statistical machine learning, linearity is defined with respect to the unknown parameters w, NOT the input variables x. Because the prediction equation y_predicted = Sum j=0^p w_j φ_j(x) is a linear combination of weights w, the cost function remains strictly quadratic and convex, allowing exact optimization via standard linear OLS and gradient descent.
As polynomial degree p increases, feature count explodes combinatorially across interactions (binom{d+p}{p}). High-degree polynomials oscillate wildly near interval boundaries (Runge's phenomenon), causing severe overfitting and necessitating regularization.
While the analytical Normal Equation w^* = (X^TX)-1X^Ty is exact, inverting the (d+1) × (d+1) matrix has a computational complexity of O(d^3). When features exceed d > 10,000 (e.g., text representations, genomics), matrix inversion becomes computationally intractable. We employ Gradient Descent.
| Gradient Descent Variant | Batch Size Per Step | Computational Efficiency | Convergence Properties |
|---|---|---|---|
| Batch Gradient Descent | All N training samples |
Slow on large N; memory intensive |
Smooth, monotonic decrease directly to global minimum. |
| Stochastic Gradient Descent (SGD) | Single random sample (N=1) |
Extremely fast; minimal memory footprint | Noisy, erratic updates; can jump out of shallow local minima. |
| Mini-Batch Gradient Descent | Mini-batch (e.g., B in [32, 256]) |
Optimal GPU vectorized SIMD parallelism | Standard industry compromise: stable, fast, robust convergence. |
When features are correlated (multicollinearity) or exceed the sample size (d gg N), the matrix X^TX becomes singular or ill-conditioned (det(X^TX) ≈ 0). Inverting an ill-conditioned matrix causes weight estimates to explode to astronomical positive and negative magnitudes, yielding high variance. Regularization penalizes large weight magnitudes.
L_2 Penalty)
JRidge(w) = MSE(w) + λ w₂² = MSE + λ Sum j=1^d w_j^2wRidge = (X^TX + λ I)-1X^Tyλ I guarantees invertibility by forcing all eigenvalues to be strictly positive.
L_1 Penalty)
JLasso(w) = MSE(w) + λ w₁ = MSE + λ Sum j=1^d w_jw_j = 0). Acts as an automatic embedded feature selector.
JEN(w) = MSE + λ [ ρ w₁ + (1 - ρ / 2) w₂² right]w_1^2 + w_2^2 ≤ C is smooth and round.w_1 ≠ 0 and w_2 ≠ 0.|w_1| + |w_2| ≤ C has sharp corners on coordinate axes.w_j = 0).
| Evaluation Metric | Mathematical Formulation | Optimal Value | Critical Trade-Off & Architectural Use |
|---|---|---|---|
| Mean Absolute Error (MAE) | (1 / N)Sum i=1^N |y_i - y_predicted_i| |
0.0 |
Linear penalty. Robust to outliers; reflects median absolute error in natural target units. |
| Mean Squared Error (MSE) | (1 / N)Sum i=1^N (y_i - y_predicted_i)^2 |
0.0 |
Differentiable loss function; heavily penalizes rare catastrophic errors quadratically. |
| Root Mean Squared Error (RMSE) | √(MSE) |
0.0 |
Standard industry benchmark; expresses quadratic error in identical physical units as y. |
R-Squared (R^2) |
1 - (Sum (y_i - y_predicted_i)^2 / Sum (y_i - Mean_y)^2) |
1.0 |
Proportion of variance explained. Flawed: Monotonically increases when adding useless features. |
Adjusted R-Squared (R^2adj) |
1 - [ ((1 - R^2)(N - 1) / N - d - 1) right] |
1.0 |
Penalizes addition of non-informative features d. Increases only if new feature improves model beyond chance. |
In Ordinary Least Squares, adding any feature to design matrix X—even pure Gaussian random white noise—will never decrease raw R^2. Every additional column expands the subspace col(X), allowing OLS to project y slightly closer to the training points. Adjusted R^2 corrects for degrees of freedom (N - d - 1). If a feature does not reduce unexplained variance by more than (1 / N - d - 1), Adjusted R^2 declines.
The following standalone script demonstrates the Normal Equation vs. Scikit-learn OLS, Polynomial Regression, Ridge, Lasso sparsity, Elastic Net, and comprehensive metric calculations:
# =====================================================================
# DevAscent Academy: Verified Comprehensive Regression Algorithms Suite
# Demonstrating: Analytical Normal Equation, Ridge, Lasso, Elastic Net,
# and Complete Metric Calculations (MSE, RMSE, MAE, R2, Adj R2)
# =====================================================================
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
def compute_regression_metrics(y_true, y_pred, n_features):
"""Calculates all essential regression performance metrics."""
n_samples = len(y_true)
mae = mean_absolute_error(y_true, y_pred)
mse = mean_squared_error(y_true, y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(y_true, y_pred)
# Adjusted R2 Formula: 1 - [(1 - R2) * (N - 1) / (N - d - 1)]
adj_r2 = 1.0 - ((1.0 - r2) * (n_samples - 1.0) / (n_samples - n_features - 1.0))
return {
"MAE": mae,
"MSE": mse,
"RMSE": rmse,
"R2": r2,
"Adj_R2": adj_r2
}
def run_analytical_vs_sklearn_ols():
print("=== 1. ANALYTICAL NORMAL EQUATION VS SCIKIT-LEARN OLS ===")
np.random.seed(42)
N, d = 150, 3
X = np.random.randn(N, d)
true_weights = np.array([3.5, -2.0, 1.8])
true_bias = 4.2
y = X @ true_weights + true_bias + np.random.normal(0, 0.5, size=N)
# 1. Analytical Normal Equation: w* = (X_aug^T @ X_aug)^(-1) @ X_aug^T @ y
X_aug = np.hstack([np.ones((N, 1)), X])
w_analytical = np.linalg.inv(X_aug.T @ X_aug) @ X_aug.T @ y
# 2. Scikit-Learn LinearRegression
lr = LinearRegression(fit_intercept=True)
lr.fit(X, y)
print(f"Analytical Weights: Intercept={w_analytical[0]:.4f}, Slopes={w_analytical[1:].round(4)}")
print(f"Scikit-Learn OLS: Intercept={lr.intercept_:.4f}, Slopes={lr.coef_.round(4)}")
print(f"Verification Match: {np.allclose(w_analytical, np.insert(lr.coef_, 0, lr.intercept_))}\n")
def run_regularization_comparison():
print("=== 2. REGULARIZATION SUITE: OLS VS RIDGE VS LASSO VS ELASTIC NET ===")
np.random.seed(42)
# Generate dataset with 10 features, but only 3 are truly informative
N, d = 200, 10
X = np.random.randn(N, d)
true_w = np.array([5.0, -3.5, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
y = X @ true_w + 2.5 + np.random.normal(0, 1.2, size=N)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)
# Scale features (mandatory for regularized models)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# Initialize models
models = {
"OLS (Unregularized)": LinearRegression(),
"Ridge (L2, alpha=2.0)": Ridge(alpha=2.0),
"Lasso (L1, alpha=0.35)": Lasso(alpha=0.35),
"ElasticNet (L1+L2, alpha=0.4, l1_ratio=0.5)": ElasticNet(alpha=0.4, l1_ratio=0.5)
}
results = []
coefficients = {}
for name, model in models.items():
model.fit(X_train_scaled, y_train)
y_pred = model.predict(X_test_scaled)
metrics = compute_regression_metrics(y_test, y_pred, n_features=d)
metrics["Model"] = name
results.append(metrics)
coefficients[name] = model.coef_
results_df = pd.DataFrame(results)[["Model", "MAE", "RMSE", "R2", "Adj_R2"]]
print(results_df.to_string(index=False))
print("\n=== Learned Feature Coefficients (Checking Sparsity) ===")
coef_df = pd.DataFrame(coefficients, index=[f"Feat_{i+1}" for i in range(d)])
print(coef_df.round(4))
print("\nInsight: Notice how Lasso drives non-informative features 4-10 strictly to 0.0000!")
if __name__ == "__main__":
run_analytical_vs_sklearn_ols()
run_regularization_comparison()E[ε] = 0), are uncorrelated (Cov(ε_i, ε_j) = 0), and have constant homoscedastic variance (Var(ε_i) = σ^2), the Ordinary Least Squares estimator w^* = (X^TX)-1X^Ty is BLUE: Best Linear Unbiased Estimator. "Best" means that among all possible linear unbiased estimators, OLS achieves the strictly minimum sampling variance (smallest possible parameter covariance matrix). Notably, this theorem holds without requiring the errors to follow a normal distribution.
x_2 = 2.5 x_1), the columns of the design matrix X are linearly dependent. Consequently, the Gram matrix X^TX in R(d+1) × (d+1) is rank-deficient (rank(X^TX) < d+1) and its determinant is strictly zero (det(X^TX) = 0). A matrix with zero determinant is singular and has no inverse (X^TX)-1, rendering the analytical Normal Equation undefined. Geometrically, infinitely many weight vectors achieve the exact same minimum RSS. Ridge Regression resolves this by adding λ I, which shifts all eigenvalues by +λ, guaranteeing non-zero determinant and invertibility.
λ Sum w_j^2 or λ Sum |w_j| applies an identical, uniform penalty multiplier across all model weights. If Feature A has a tiny numerical scale (e.g., salary in Millions: 0.05 to 0.5) and Feature B has a massive scale (e.g., credit score: 300 to 850), Feature A requires an enormous weight coefficient to produce a meaningful prediction, while Feature B requires a tiny weight. The regularizer will aggressively penalize Feature A simply because of arbitrary physical units, shrinking its coefficient toward zero regardless of true predictive importance. Standardization (μ=0, σ=1) puts all features on equal footing so regularization reflects true predictive contribution.
L_1) produces sparse solutions (weights =0) while Ridge (L_2) only shrinks weights toward zero.
L_2 budget region is a smooth hypersphere (Sum w_j^2 ≤ C). The gradient of the L_2 penalty is 2λ w_j, which diminishes proportionally as w_j → 0 (decaying to zero as the weight shrinks). In contrast, the L_1 budget region is a hyper-rhombus (polytope) with sharp corners positioned exactly on the coordinate axes where orthogonal weights are zero (w_j = 0). Furthermore, the derivative of |w_j| is the sign function sign(w_j), which exerts a constant penalty force ± λ regardless of how small w_j becomes, pulling the parameter all the way to zero and locking it there.
R^2 decrease when adding a new feature, whereas raw R^2 can only increase or stay constant?
R^2 = 1 - (SSres / SStot). In OLS, adding any feature column expands the column space of X; because projection onto a larger subspace can never increase distance to target y, SSres can never increase, so raw R^2 artificially creeps upward even for random noise. Adjusted R^2 penalizes degrees of freedom:
R^2adj = 1 - [ ((1 - R^2)(N - 1) / N - d - 1) right]d → d + 1), the denominator (N - d - 1) decreases, which scales up the penalty term. Unless the reduction in (1 - R^2) outweighs this penalty scaling factor, Adjusted R^2 will decline, alerting engineers that the new feature adds more model complexity than predictive power.
d > N): Lasso can select at most N features before becoming saturated, whereas Elastic Net can select all relevant features.L_2 component provides a grouping effect, shrinking the correlated group together while its L_1 component removes entirely unrelated features.(y - y_predicted)^2), a single extreme prediction outlier (e.g., an error of 100 vs typical error of 2) contributes 10,000 to the sum, dominating the objective function and distorting the regression line away from typical inliers. MAE penalizes all errors linearly (|y - y_predicted|), providing a median-centric, robust measure of expected error in true problem domain units.
w.
J(w) = (1 / 2N) (Xw - y)^T (Xw - y). Expanding:
J(w) = (1 / 2N) Big( w^TX^TXw - w^TX^Ty - y^TXw + y^Ty Big)w^TX^Ty is a scalar, its transpose equals itself: w^TX^Ty = y^TXw. Thus:
J(w) = (1 / 2N) Big( w^T(X^TX)w - 2y^TXw + y^Ty Big)∇w(w^TAw) = 2Aw (for symmetric A = X^TX) and ∇w(b^Tw) = b:
∇w J(w) = (1 / 2N) Big( 2X^TXw - 2X^Ty Big) = (1 / N) X^T(Xw - y)