linear_regression.tex

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{siunitx}
\usepackage{float}
\usepackage{geometry}
\geometry{margin=1in}
\usepackage[makestderr]{pythontex}

\title{Linear Regression: OLS, Gradient Descent, and Regularization}
\author{Machine Learning Foundations}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This document presents a comprehensive analysis of linear regression methods including ordinary least squares (OLS), gradient descent optimization, and regularization techniques (Ridge and Lasso). We examine model diagnostics, multicollinearity detection, and cross-validation for hyperparameter tuning.
\end{abstract}

\section{Introduction}
Linear regression models the relationship between features $X$ and target $y$:
\begin{equation}
y = X\beta + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2 I)
\end{equation}

\textbf{OLS Solution:}
\begin{equation}
\hat{\beta}_{OLS} = (X^T X)^{-1} X^T y
\end{equation}

\textbf{Ridge Regression:}
\begin{equation}
\hat{\beta}_{Ridge} = (X^T X + \lambda I)^{-1} X^T y
\end{equation}

\textbf{Lasso Regression:}
\begin{equation}
\hat{\beta}_{Lasso} = \arg\min_\beta \|y - X\beta\|_2^2 + \lambda \|\beta\|_1
\end{equation}

\section{Computational Environment}
\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

plt.rc('text', usetex=True)
plt.rc('font', family='serif')
np.random.seed(42)

def save_plot(filename, caption):
    plt.savefig(filename, bbox_inches='tight', dpi=150)
    print(r'\begin{figure}[H]')
    print(r'\centering')
    print(r'\includegraphics[width=0.85\textwidth]{' + filename + '}')
    print(r'\caption{' + caption + '}')
    print(r'\end{figure}')
    plt.close()
\end{pycode}

\section{Data Generation}
\begin{pycode}
# Generate regression dataset with known coefficients
n_samples = 200
n_features = 5
true_coefs = np.array([3.0, -2.0, 1.5, 0, 0])  # Last two are irrelevant

# Generate features with some correlation
X = np.random.randn(n_samples, n_features)
X[:, 1] = X[:, 0] * 0.5 + X[:, 1] * 0.5  # Introduce correlation

# Generate target with noise
y = X @ true_coefs + np.random.normal(0, 1, n_samples)

# Split data
n_train = 150
X_train, X_test = X[:n_train], X[n_train:]
y_train, y_test = y[:n_train], y[n_train:]

# Standardize features
X_mean = X_train.mean(axis=0)
X_std = X_train.std(axis=0)
X_train_scaled = (X_train - X_mean) / X_std
X_test_scaled = (X_test - X_mean) / X_std

# Plot data overview
fig, axes = plt.subplots(1, 3, figsize=(12, 4))

# Feature distributions
for i in range(n_features):
    axes[0].hist(X_train[:, i], bins=20, alpha=0.5, label=f'X{i+1}')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Frequency')
axes[0].set_title('Feature Distributions')
axes[0].legend()

# Correlation matrix
corr = np.corrcoef(X_train.T)
im = axes[1].imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
axes[1].set_xticks(range(n_features))
axes[1].set_yticks(range(n_features))
axes[1].set_xticklabels([f'X{i+1}' for i in range(n_features)])
axes[1].set_yticklabels([f'X{i+1}' for i in range(n_features)])
axes[1].set_title('Feature Correlation')
plt.colorbar(im, ax=axes[1])

# Target distribution
axes[2].hist(y_train, bins=20, color='steelblue', edgecolor='black', alpha=0.7)
axes[2].set_xlabel('y')
axes[2].set_ylabel('Frequency')
axes[2].set_title('Target Distribution')

plt.tight_layout()
save_plot('lr_data.pdf', 'Dataset overview showing feature distributions and correlations.')
\end{pycode}

Dataset: $n = \py{n_samples}$ samples, $p = \py{n_features}$ features.

\section{Ordinary Least Squares}
\begin{pycode}
class LinearRegression:
    def __init__(self):
        self.coef_ = None
        self.intercept_ = None

    def fit(self, X, y):
        # Add intercept
        X_b = np.column_stack([np.ones(len(X)), X])
        # OLS solution
        self.coef_ = np.linalg.lstsq(X_b, y, rcond=None)[0]
        self.intercept_ = self.coef_[0]
        self.coef_ = self.coef_[1:]
        return self

    def predict(self, X):
        return X @ self.coef_ + self.intercept_

    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

# Fit OLS model
ols = LinearRegression()
ols.fit(X_train_scaled, y_train)

y_train_pred = ols.predict(X_train_scaled)
y_test_pred = ols.predict(X_test_scaled)

train_r2 = ols.score(X_train_scaled, y_train)
test_r2 = ols.score(X_test_scaled, y_test)
train_mse = np.mean((y_train - y_train_pred)**2)
test_mse = np.mean((y_test - y_test_pred)**2)

# Residual analysis
residuals = y_train - y_train_pred

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# Predicted vs Actual
axes[0, 0].scatter(y_train, y_train_pred, alpha=0.5, s=30)
lims = [min(y_train.min(), y_train_pred.min()), max(y_train.max(), y_train_pred.max())]
axes[0, 0].plot(lims, lims, 'r--', linewidth=2)
axes[0, 0].set_xlabel('Actual')
axes[0, 0].set_ylabel('Predicted')
axes[0, 0].set_title(f'Predicted vs Actual (R2 = {train_r2:.3f})')
axes[0, 0].grid(True, alpha=0.3)

# Residuals vs Fitted
axes[0, 1].scatter(y_train_pred, residuals, alpha=0.5, s=30)
axes[0, 1].axhline(0, color='r', linestyle='--')
axes[0, 1].set_xlabel('Fitted Values')
axes[0, 1].set_ylabel('Residuals')
axes[0, 1].set_title('Residuals vs Fitted')
axes[0, 1].grid(True, alpha=0.3)

# Q-Q plot
stats.probplot(residuals, dist="norm", plot=axes[1, 0])
axes[1, 0].set_title('Normal Q-Q Plot')
axes[1, 0].grid(True, alpha=0.3)

# Residual histogram
axes[1, 1].hist(residuals, bins=20, density=True, alpha=0.7, color='steelblue', edgecolor='black')
x_norm = np.linspace(residuals.min(), residuals.max(), 100)
axes[1, 1].plot(x_norm, stats.norm.pdf(x_norm, 0, np.std(residuals)), 'r-', linewidth=2)
axes[1, 1].set_xlabel('Residual')
axes[1, 1].set_ylabel('Density')
axes[1, 1].set_title('Residual Distribution')

plt.tight_layout()
save_plot('lr_diagnostics.pdf', 'OLS regression diagnostics showing residual analysis.')
\end{pycode}

OLS Performance: Train $R^2 = \py{f"{train_r2:.3f}"}$, Test $R^2 = \py{f"{test_r2:.3f}"}$.

\section{Gradient Descent}
\begin{pycode}
class GradientDescentRegression:
    def __init__(self, learning_rate=0.01, n_iterations=1000):
        self.lr = learning_rate
        self.n_iter = n_iterations
        self.coef_ = None
        self.intercept_ = None
        self.loss_history = []

    def fit(self, X, y):
        n_samples, n_features = X.shape
        self.coef_ = np.zeros(n_features)
        self.intercept_ = 0
        self.loss_history = []

        for _ in range(self.n_iter):
            y_pred = X @ self.coef_ + self.intercept_
            error = y_pred - y
            loss = np.mean(error**2)
            self.loss_history.append(loss)

            # Gradients
            grad_coef = (2/n_samples) * X.T @ error
            grad_intercept = (2/n_samples) * np.sum(error)

            # Update
            self.coef_ -= self.lr * grad_coef
            self.intercept_ -= self.lr * grad_intercept

        return self

    def predict(self, X):
        return X @ self.coef_ + self.intercept_

# Fit with gradient descent
gd = GradientDescentRegression(learning_rate=0.1, n_iterations=500)
gd.fit(X_train_scaled, y_train)

# Compare learning rates
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

learning_rates = [0.001, 0.01, 0.1, 0.5]
for lr in learning_rates:
    gd_temp = GradientDescentRegression(learning_rate=lr, n_iterations=200)
    gd_temp.fit(X_train_scaled, y_train)
    axes[0].plot(gd_temp.loss_history, label=f'lr={lr}', linewidth=1.5)

axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('MSE Loss')
axes[0].set_title('Gradient Descent Convergence')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_yscale('log')

# Coefficient comparison
x_pos = np.arange(n_features)
width = 0.35
axes[1].bar(x_pos - width/2, ols.coef_, width, label='OLS', alpha=0.8)
axes[1].bar(x_pos + width/2, gd.coef_, width, label='GD', alpha=0.8)
axes[1].set_xticks(x_pos)
axes[1].set_xticklabels([f'X{i+1}' for i in range(n_features)])
axes[1].set_xlabel('Feature')
axes[1].set_ylabel('Coefficient')
axes[1].set_title('OLS vs Gradient Descent Coefficients')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
save_plot('lr_gradient_descent.pdf', 'Gradient descent convergence and coefficient comparison with OLS.')
\end{pycode}

\section{Ridge Regression}
\begin{pycode}
class RidgeRegression:
    def __init__(self, alpha=1.0):
        self.alpha = alpha
        self.coef_ = None
        self.intercept_ = None

    def fit(self, X, y):
        n_features = X.shape[1]
        X_b = np.column_stack([np.ones(len(X)), X])
        I = np.eye(n_features + 1)
        I[0, 0] = 0  # Don't regularize intercept
        coef = np.linalg.solve(X_b.T @ X_b + self.alpha * I, X_b.T @ y)
        self.intercept_ = coef[0]
        self.coef_ = coef[1:]
        return self

    def predict(self, X):
        return X @ self.coef_ + self.intercept_

    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

# Ridge path
alphas = np.logspace(-3, 3, 50)
ridge_coefs = []
ridge_scores = []

for alpha in alphas:
    ridge = RidgeRegression(alpha=alpha)
    ridge.fit(X_train_scaled, y_train)
    ridge_coefs.append(ridge.coef_)
    ridge_scores.append(ridge.score(X_test_scaled, y_test))

ridge_coefs = np.array(ridge_coefs)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Coefficient path
for i in range(n_features):
    axes[0].plot(alphas, ridge_coefs[:, i], label=f'X{i+1}', linewidth=1.5)
axes[0].set_xscale('log')
axes[0].set_xlabel('Alpha (Regularization)')
axes[0].set_ylabel('Coefficient Value')
axes[0].set_title('Ridge Regression Coefficient Path')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Validation curve
axes[1].plot(alphas, ridge_scores, 'b-o', linewidth=1.5, markersize=4)
best_alpha_idx = np.argmax(ridge_scores)
best_alpha = alphas[best_alpha_idx]
axes[1].axvline(best_alpha, color='r', linestyle='--', label=f'Best alpha={best_alpha:.3f}')
axes[1].set_xscale('log')
axes[1].set_xlabel('Alpha')
axes[1].set_ylabel('Test R2 Score')
axes[1].set_title('Ridge Regression Validation Curve')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
save_plot('lr_ridge.pdf', 'Ridge regression coefficient path and validation curve.')

best_ridge = RidgeRegression(alpha=best_alpha)
best_ridge.fit(X_train_scaled, y_train)
ridge_test_r2 = best_ridge.score(X_test_scaled, y_test)
\end{pycode}

Best Ridge alpha: \py{f"{best_alpha:.3f}"} with test $R^2 = \py{f"{ridge_test_r2:.3f}"}$.

\section{Lasso Regression}
\begin{pycode}
class LassoRegression:
    def __init__(self, alpha=1.0, n_iterations=1000):
        self.alpha = alpha
        self.n_iter = n_iterations
        self.coef_ = None
        self.intercept_ = None

    def _soft_threshold(self, x, threshold):
        return np.sign(x) * np.maximum(np.abs(x) - threshold, 0)

    def fit(self, X, y):
        n_samples, n_features = X.shape
        self.coef_ = np.zeros(n_features)
        self.intercept_ = np.mean(y)

        # Coordinate descent
        for _ in range(self.n_iter):
            for j in range(n_features):
                residual = y - self.intercept_ - X @ self.coef_ + X[:, j] * self.coef_[j]
                rho = X[:, j] @ residual
                z = np.sum(X[:, j]**2)
                self.coef_[j] = self._soft_threshold(rho, self.alpha * n_samples / 2) / z

            self.intercept_ = np.mean(y - X @ self.coef_)

        return self

    def predict(self, X):
        return X @ self.coef_ + self.intercept_

    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

# Lasso path
lasso_alphas = np.logspace(-4, 0, 50)
lasso_coefs = []
lasso_scores = []

for alpha in lasso_alphas:
    lasso = LassoRegression(alpha=alpha, n_iterations=500)
    lasso.fit(X_train_scaled, y_train)
    lasso_coefs.append(lasso.coef_)
    lasso_scores.append(lasso.score(X_test_scaled, y_test))

lasso_coefs = np.array(lasso_coefs)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Coefficient path
for i in range(n_features):
    axes[0].plot(lasso_alphas, lasso_coefs[:, i], label=f'X{i+1}', linewidth=1.5)
axes[0].set_xscale('log')
axes[0].set_xlabel('Alpha (Regularization)')
axes[0].set_ylabel('Coefficient Value')
axes[0].set_title('Lasso Regression Coefficient Path')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Validation curve
axes[1].plot(lasso_alphas, lasso_scores, 'g-s', linewidth=1.5, markersize=4)
best_lasso_idx = np.argmax(lasso_scores)
best_lasso_alpha = lasso_alphas[best_lasso_idx]
axes[1].axvline(best_lasso_alpha, color='r', linestyle='--', label=f'Best alpha={best_lasso_alpha:.4f}')
axes[1].set_xscale('log')
axes[1].set_xlabel('Alpha')
axes[1].set_ylabel('Test R2 Score')
axes[1].set_title('Lasso Regression Validation Curve')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
save_plot('lr_lasso.pdf', 'Lasso regression coefficient path showing feature selection.')

best_lasso = LassoRegression(alpha=best_lasso_alpha, n_iterations=500)
best_lasso.fit(X_train_scaled, y_train)
lasso_test_r2 = best_lasso.score(X_test_scaled, y_test)
n_selected = np.sum(np.abs(best_lasso.coef_) > 0.01)
\end{pycode}

Best Lasso alpha: \py{f"{best_lasso_alpha:.4f}"} with test $R^2 = \py{f"{lasso_test_r2:.3f}"}$. Selected \py{n_selected} features.

\section{Model Comparison}
\begin{pycode}
# Compare all methods
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Coefficient comparison
methods = ['True', 'OLS', 'Ridge', 'Lasso']
coefs = [true_coefs, ols.coef_, best_ridge.coef_, best_lasso.coef_]

x = np.arange(n_features)
width = 0.2

for i, (method, coef) in enumerate(zip(methods, coefs)):
    axes[0].bar(x + i*width, coef, width, label=method, alpha=0.8)

axes[0].set_xticks(x + 1.5*width)
axes[0].set_xticklabels([f'X{i+1}' for i in range(n_features)])
axes[0].set_xlabel('Feature')
axes[0].set_ylabel('Coefficient')
axes[0].set_title('Coefficient Comparison Across Methods')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Performance comparison
methods_names = ['OLS', 'Ridge', 'Lasso']
test_scores = [test_r2, ridge_test_r2, lasso_test_r2]
train_scores = [train_r2, best_ridge.score(X_train_scaled, y_train),
                best_lasso.score(X_train_scaled, y_train)]

x = np.arange(len(methods_names))
width = 0.35
axes[1].bar(x - width/2, train_scores, width, label='Train', alpha=0.8)
axes[1].bar(x + width/2, test_scores, width, label='Test', alpha=0.8)
axes[1].set_xticks(x)
axes[1].set_xticklabels(methods_names)
axes[1].set_ylabel('R2 Score')
axes[1].set_title('Model Performance Comparison')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
save_plot('lr_comparison.pdf', 'Comparison of OLS, Ridge, and Lasso regression methods.')
\end{pycode}

\section{Results Summary}
\begin{table}[H]
\centering
\caption{Regression Model Performance}
\begin{tabular}{lccc}
\toprule
Model & Train $R^2$ & Test $R^2$ & Hyperparameter \\
\midrule
OLS & \py{f"{train_r2:.3f}"} & \py{f"{test_r2:.3f}"} & - \\
Ridge & \py{f"{best_ridge.score(X_train_scaled, y_train):.3f}"} & \py{f"{ridge_test_r2:.3f}"} & $\alpha = \py{f"{best_alpha:.3f}"}$ \\
Lasso & \py{f"{best_lasso.score(X_train_scaled, y_train):.3f}"} & \py{f"{lasso_test_r2:.3f}"} & $\alpha = \py{f"{best_lasso_alpha:.4f}"}$ \\
\bottomrule
\end{tabular}
\end{table}

\begin{pycode}
print(r'\begin{table}[H]')
print(r'\centering')
print(r'\caption{Coefficient Estimates}')
print(r'\begin{tabular}{ccccc}')
print(r'\toprule')
print(r'Feature & True & OLS & Ridge & Lasso \\')
print(r'\midrule')
for i in range(n_features):
    print(f'X{i+1} & {true_coefs[i]:.2f} & {ols.coef_[i]:.2f} & {best_ridge.coef_[i]:.2f} & {best_lasso.coef_[i]:.2f} \\\\')
print(r'\bottomrule')
print(r'\end{tabular}')
print(r'\end{table}')
\end{pycode}

\section{Conclusion}
This analysis demonstrated:
\begin{itemize}
    \item OLS as the baseline with closed-form solution
    \item Gradient descent as iterative optimization method
    \item Ridge regression for handling multicollinearity (L2 penalty)
    \item Lasso regression for feature selection (L1 penalty)
    \item Hyperparameter tuning via validation curves
\end{itemize}

The Lasso method successfully identified the relevant features and set irrelevant coefficients close to zero, demonstrating its utility for sparse solutions.

\end{document}