Linear & Logistic Regression

Before fitting any model, we need a precise vocabulary for the problem.

We are given a training set of $n$ examples: ${(x^{(i)}, y^{(i)})}_{i=1}^{n}$, where $x^{(i)} \in \mathbb{R}^d$ is a feature vector (a $d$-dimensional description of the $i$-th example) and $y^{(i)} \in \mathbb{R}$ is the label or target we want to predict.

We want to find a function $h: \mathbb{R}^d \to \mathbb{R}$ — called the hypothesis — that maps new feature vectors to predictions. For linear regression, we restrict $h$ to the family of linear (technically, affine) functions:

$$h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \cdots + \theta_d x_d$$

The $\theta_j$ are the parameters or weights of the model. We add a bias (intercept) term $\theta_0$ by prepending a constant $x_0 = 1$ to every feature vector, so the whole thing collapses to:

$$h_\theta(x) = \theta^\top x$$

where now $x \in \mathbb{R}^{d+1}$ and $\theta \in \mathbb{R}^{d+1}$.

The learning problem: find $\theta$ that makes $h_\theta(x^{(i)})$ close to $y^{(i)}$ for all training examples — and, crucially, also generalizes to unseen data.

Analogy: think of your features as dials on a mixing board and $\theta$ as the dial settings. You are tuning the board so the output signal matches a target recording.

import numpy as np

# n=4 examples, d=2 features (house size, bedrooms)
# We prepend a column of 1s for the bias term
X_raw = np.array([[2104, 3],
                   [1600, 3],
                   [2400, 4],
                   [1416, 2]])  # shape (4, 2)

n = X_raw.shape[0]
X = np.hstack([np.ones((n, 1)), X_raw])  # shape (4, 3): [1, x1, x2]
y = np.array([400, 330, 369, 232])        # prices in $k

print(X.shape)  # (4, 3)
print(y.shape)  # (4,)

We need a way to measure how well our current $\theta$ fits the data. The standard choice is the mean squared error (MSE) cost:

$$J(\theta) = \frac{1}{2n} \sum_{i=1}^{n} \left(h_\theta(x^{(i)}) - y^{(i)}\right)^2$$

The factor of $\frac{1}{2}$ is purely for algebraic convenience (it cancels the 2 from the derivative). Minimizing $J(\theta)$ is called least squares regression — we are minimizing the sum of squared residuals.

Why squared and not absolute value? Several reasons:

1. Squaring is differentiable everywhere, making calculus-based optimization clean.
2. It penalizes large errors more heavily than small ones — a feature if you believe large errors are disproportionately bad.
3. It connects to maximum likelihood estimation under Gaussian noise (derived later).

Geometric intuition: Each residual $h_\theta(x^{(i)}) - y^{(i)}$ is a vertical distance from the prediction line to the actual data point. We are minimizing the sum of the squares of these vertical distances. Imagine pushing the regression line around until the total area of the squares drawn on those residuals is minimized.

In matrix notation, with design matrix $X \in \mathbb{R}^{n \times (d+1)}$ and label vector $y \in \mathbb{R}^n$:

$$J(\theta) = \frac{1}{2n}|X\theta - y|^2 = \frac{1}{2n}(X\theta - y)^\top(X\theta - y)$$

def cost(X, y, theta):
    """MSE cost for linear regression."""
    n = len(y)
    residuals = X @ theta - y      # shape (n,)
    return (1 / (2 * n)) * np.dot(residuals, residuals)

# Random initial theta: very bad cost
theta_init = np.zeros(3)
print(cost(X, y, theta_init))  # large number

Because $J(\theta)$ is a convex quadratic function of $\theta$, it has a unique global minimum. We can find it analytically by setting the gradient to zero.

Expand the cost: $$J(\theta) = \frac{1}{2n}(X\theta - y)^\top(X\theta - y) = \frac{1}{2n}(\theta^\top X^\top X \theta - 2y^\top X\theta + y^\top y)$$

Take the gradient with respect to $\theta$ (using the matrix calculus identities $\nabla_\theta(\theta^\top A\theta) = 2A\theta$ for symmetric $A$, and $\nabla_\theta(b^\top\theta) = b$):

$$\nabla_\theta J = \frac{1}{n}(X^\top X \theta - X^\top y)$$

Set to zero: $$X^\top X \theta^* = X^\top y$$

If $X^\top X$ is invertible: $$\boxed{\theta^* = (X^\top X)^{-1} X^\top y}$$

This is the normal equation. The matrix $(X^\top X)^{-1}X^\top$ is called the Moore-Penrose pseudoinverse of $X$.

Geometric interpretation: $X\theta^$ is the orthogonal projection of $y$ onto the column space of $X$ — the closest point in that subspace to $y$. The residual vector $y - X\theta^$ is perpendicular to every column of $X$, which is exactly what $X^\top(y - X\theta^*) = 0$ says.

When does $X^\top X$ fail to be invertible? When features are linearly dependent (e.g., one feature is a linear combination of others). In practice we use numpy.linalg.lstsq or regularization (ridge regression adds $\lambda I$ to $X^\top X$, always making it invertible).

Complexity: Computing $(X^\top X)^{-1}$ costs $O(d^3)$ and the matrix multiply $X^\top X$ costs $O(nd^2)$. For $d \ll n$, this is fine. For $d = 10^6$ (e.g., NLP), iterative methods like gradient descent are far preferable.

def normal_equation(X, y):
    """Closed-form OLS solution: theta = (X^T X)^{-1} X^T y"""
    return np.linalg.lstsq(X, y, rcond=None)[0]

theta_opt = normal_equation(X, y)
print(theta_opt)   # [bias, w_sqft, w_bedrooms]

# Verify: cost should be near zero on training data
print(cost(X, y, theta_opt))

The normal equation is elegant but requires $O(d^3)$ matrix inversion. For large-scale problems, gradient descent is the workhorse: iteratively walk in the direction of steepest descent of $J$.

Update rule (for all parameters simultaneously): $$\theta_j \leftarrow \theta_j - \alpha \frac{\partial J}{\partial \theta_j}$$

For linear regression, the partial derivative is: $$\frac{\partial J}{\partial \theta_j} = \frac{1}{n}\sum_{i=1}^{n}\left(h_\theta(x^{(i)}) - y^{(i)}\right)x_j^{(i)}$$

In matrix form: $\nabla_\theta J = \frac{1}{n}X^\top(X\theta - y)$, so the update is: $$\theta \leftarrow \theta - \frac{\alpha}{n} X^\top(X\theta - y)$$

Three flavors:

• Batch GD: use all $n$ examples per step. Stable but slow per epoch when $n$ is huge.
• Stochastic GD (SGD): use one random example per step. Noisy but very fast updates; often escapes saddle points.
• Mini-batch GD: use a batch of $B \approx 32$–$256$ examples. Balances stability and speed; the standard in deep learning.

Learning rate $\alpha$: the critical hyperparameter. Too large: diverges. Too small: converges slowly. Strategies: fixed $\alpha$, exponential decay, Adam optimizer (adaptive per-parameter rates).

Feature normalization is essential: if $x_1 \in [0, 2000]$ (sqft) and $x_2 \in [1, 5]$ (bedrooms), the cost surface is a narrow elongated ellipse and gradient descent oscillates. After normalizing (zero mean, unit variance), it is nearly circular and converges in far fewer steps.

def gradient_descent(X, y, alpha=0.01, n_iters=1000):
    n, d = X.shape
    theta = np.zeros(d)
    history = []
    for _ in range(n_iters):
        grad = (1 / n) * X.T @ (X @ theta - y)  # shape (d,)
        theta -= alpha * grad
        history.append(cost(X, y, theta))
    return theta, history

# Always normalize features before gradient descent!
X_norm = X.copy().astype(float)
X_norm[:, 1:] = (X_norm[:, 1:] - X_norm[:, 1:].mean(0)) / X_norm[:, 1:].std(0)

theta_gd, hist = gradient_descent(X_norm, y, alpha=0.1, n_iters=1000)
print(f"Final cost: {hist[-1]:.4f}")

Why is the squared-error loss the "right" loss for regression? The answer comes from maximum likelihood estimation (MLE).

Assume the target is generated by our linear model plus independent Gaussian noise: $$y^{(i)} = \theta^\top x^{(i)} + \epsilon^{(i)}, \quad \epsilon^{(i)} \sim \mathcal{N}(0, \sigma^2)$$

This means the likelihood of observing $y^{(i)}$ given $x^{(i)}$ and $\theta$ is: $$p(y^{(i)} \mid x^{(i)}; \theta) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y^{(i)} - \theta^\top x^{(i)})^2}{2\sigma^2}\right)$$

The joint likelihood over all $n$ independent examples is: $$L(\theta) = \prod_{i=1}^{n} p(y^{(i)} \mid x^{(i)}; \theta)$$

Maximizing $L(\theta)$ is equivalent to minimizing its negative log-likelihood: $$-\log L(\theta) = \frac{n}{2}\log(2\pi\sigma^2) + \frac{1}{2\sigma^2}\sum_{i=1}^n (y^{(i)} - \theta^\top x^{(i)})^2$$

Dropping constants, we recover exactly the squared-error cost! Least squares is MLE under Gaussian noise. This justifies the MSE loss not just as a computational convenience but as the principled probabilistic choice when you believe noise is Gaussian and i.i.d.

This framing also gives us the Gauss-Markov theorem: among all linear unbiased estimators, OLS has minimum variance — it is BLUE (Best Linear Unbiased Estimator).

Now shift from predicting continuous values to predicting class membership. The label $y^{(i)} \in {0, 1}$ (binary classification: spam/not-spam, malignant/benign, etc.).

Naive approach: apply linear regression and threshold at 0.5. This fails badly:

1. $h_\theta(x)$ can output values like $-3$ or $7$ — not interpretable as probabilities.
2. The squared-error loss treats a prediction of $3$ for $y=1$ as a huge error (residual $= -2$), discouraging confident correct predictions.
3. A single outlier far from the boundary can drag the regression line and destroy the threshold.

What we really want is a model that outputs $p(y=1 \mid x; \theta)$ — a probability between 0 and 1. We achieve this by passing the linear combination through the sigmoid function (also called the logistic function):

$$\sigma(z) = \frac{1}{1 + e^{-z}}$$

The hypothesis becomes: $$h_\theta(x) = \sigma(\theta^\top x) = \frac{1}{1 + e^{-\theta^\top x}}$$

Properties of the sigmoid:

• $\sigma(z) \in (0, 1)$ for all $z \in \mathbb{R}$.
• $\sigma(0) = 0.5$.
• $\sigma(-z) = 1 - \sigma(z)$ (antisymmetric around 0.5).
• $\sigma'(z) = \sigma(z)(1 - \sigma(z))$ (elegant self-referential derivative).

Think of it as a soft step function: it smoothly transitions from 0 to 1 as the linear score crosses 0.

import numpy as np
import matplotlib.pyplot as plt

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

z = np.linspace(-8, 8, 200)
print(sigmoid(0))    # 0.5
print(sigmoid(2))    # ~0.88
print(sigmoid(-2))   # ~0.12

# The derivative identity: sigma'(z) = sigma(z) * (1 - sigma(z))
s = sigmoid(z)
dsigma = s * (1 - s)  # peaks at z=0 with value 0.25

We interpret the logistic regression output as a conditional probability:

$$P(y=1 \mid x; \theta) = h_\theta(x) = \sigma(\theta^\top x)$$ $$P(y=0 \mid x; \theta) = 1 - h_\theta(x)$$

We predict class 1 when $h_\theta(x) \geq 0.5$, which is equivalent to $\theta^\top x \geq 0$.

The decision boundary is the set of $x$ where $\theta^\top x = 0$ — a hyperplane (a line in 2D, a plane in 3D). This is determined entirely by $\theta$, not by the sigmoid; the sigmoid only converts the score to a probability.

Non-linear decision boundaries: Logistic regression can represent non-linear boundaries by adding polynomial features. If we engineer $\phi(x) = [1, x_1, x_2, x_1^2, x_1 x_2, x_2^2]$, the decision boundary $\theta^\top \phi(x) = 0$ is a conic section (circle, ellipse, hyperbola) in the original feature space. This is feature engineering, not a change to the model.

The model is linear in the parameters, not necessarily in the features. This distinction is fundamental — it separates linear models from truly nonlinear ones (neural networks, kernel SVMs).

Combined with a probabilistic output, this makes logistic regression a generalized linear model (GLM) — a linear predictor passed through a link function (the sigmoid) to produce a distribution parameter.

# Decision boundary visualization (2D toy example)
# theta = [bias, w1, w2]; boundary is: bias + w1*x1 + w2*x2 = 0
# => x2 = -(bias + w1*x1) / w2

def decision_boundary_line(theta, x1_range):
    """Returns x2 values along the decision boundary."""
    bias, w1, w2 = theta
    return -(bias + w1 * x1_range) / w2

theta_ex = np.array([-3.0, 1.0, 1.0])  # boundary: x1 + x2 = 3
x1 = np.array([0, 1, 2, 3])
x2_boundary = decision_boundary_line(theta_ex, x1)
print(x2_boundary)  # [3. 2. 1. 0.]  — the line x1+x2=3

We cannot use squared error for logistic regression. If we did, $J(\theta)$ would be non-convex (many local minima), and gradient descent would fail. Instead, we derive the correct loss from maximum likelihood.

The likelihood of one example, combining the two cases elegantly: $$p(y \mid x; \theta) = h_\theta(x)^y \cdot (1 - h_\theta(x))^{1-y}$$

Verify: when $y=1$ this gives $h_\theta(x)$; when $y=0$ it gives $1 - h_\theta(x)$.

The joint log-likelihood over all $n$ examples is: $$\log L(\theta) = \sum_{i=1}^n \left[ y^{(i)} \log h_\theta(x^{(i)}) + (1 - y^{(i)}) \log(1 - h_\theta(x^{(i)})) \right]$$

We maximize this (or equivalently minimize the binary cross-entropy loss): $$J(\theta) = -\frac{1}{n}\sum_{i=1}^n \left[ y^{(i)} \log h_\theta(x^{(i)}) + (1 - y^{(i)}) \log(1 - h_\theta(x^{(i)})) \right]$$

Intuition: When $y=1$, we want $h_\theta(x) \to 1$; the term $-\log(h)$ penalizes predictions far from 1 (the cost explodes as $h \to 0$). When $y=0$, the term $-\log(1-h)$ explodes as $h \to 1$. This asymmetric, logarithmic penalty is exactly what creates a convex, well-behaved loss surface.

A remarkable fact: the gradient of the cross-entropy loss for logistic regression has the same form as linear regression: $$\frac{\partial J}{\partial \theta_j} = \frac{1}{n}\sum_{i=1}^{n}(h_\theta(x^{(i)}) - y^{(i)})x_j^{(i)}$$

The only difference is that $h$ is now the sigmoid rather than a linear function. This is no coincidence — it follows from the exponential family structure of the Bernoulli distribution.

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

def logistic_cost(X, y, theta):
    n = len(y)
    h = sigmoid(X @ theta)           # shape (n,)
    # clip to avoid log(0)
    h = np.clip(h, 1e-12, 1 - 1e-12)
    return -(1/n) * (y @ np.log(h) + (1 - y) @ np.log(1 - h))

def logistic_gradient(X, y, theta):
    n = len(y)
    h = sigmoid(X @ theta)
    return (1/n) * X.T @ (h - y)     # same form as linear regression!

def logistic_gradient_descent(X, y, alpha=0.1, n_iters=1000):
    theta = np.zeros(X.shape[1])
    for _ in range(n_iters):
        theta -= alpha * logistic_gradient(X, y, theta)
    return theta

Both linear and logistic regression overfit when features are many relative to examples, or when features are correlated. The solution is regularization — adding a penalty term to the cost that discourages large weights.

Ridge regression (L2): adds $\frac{\lambda}{2n}|\theta|2^2 = \frac{\lambda}{2n}\sum{j=1}^d \theta_j^2$ (exclude $\theta_0$).

$$J_{\text{ridge}}(\theta) = J(\theta) + \frac{\lambda}{2n}\sum_{j=1}^d \theta_j^2$$

The regularized normal equation becomes $\theta^* = (X^\top X + \lambda I')^{-1} X^\top y$, where $I'$ has 0 in the bias position. Adding $\lambda I'$ makes the matrix always invertible — an added bonus.

Bayesian interpretation: ridge regression is MAP estimation with a Gaussian prior $\theta_j \sim \mathcal{N}(0, \sigma^2/\lambda)$. The regularizer encodes a prior belief that weights are small.

Lasso regression (L1): adds $\frac{\lambda}{n}|\theta|1 = \frac{\lambda}{n}\sum{j=1}^d |\theta_j|$.

Lasso has the remarkable property of producing sparse solutions — many $\theta_j$ are driven exactly to zero — making it a built-in feature selector. This occurs because the L1 ball has corners at the axes, and the optimal solution often hits a corner.

Choosing $\lambda$: via cross-validation. Split training data into $k$ folds; for each candidate $\lambda$, train on $k-1$ folds and validate on the remaining fold; pick $\lambda$ minimizing average validation loss.

The bias-variance tradeoff: as $\lambda$ increases, model variance decreases (less overfitting) but bias increases (underfitting). The optimal $\lambda$ balances these.

def ridge_normal_equation(X, y, lam=1.0):
    """Ridge: (X^T X + lambda * I')^{-1} X^T y"""
    n, d = X.shape
    I_prime = np.eye(d)
    I_prime[0, 0] = 0  # do not regularize bias
    A = X.T @ X + lam * I_prime
    return np.linalg.solve(A, X.T @ y)

# from sklearn for comparison
from sklearn.linear_model import Ridge, Lasso
ridge = Ridge(alpha=1.0).fit(X[:, 1:], y)  # sklearn handles bias internally
print(ridge.coef_, ridge.intercept_)

Logistic regression handles binary labels ($K=2$ classes). For $K > 2$ classes, the natural generalization is softmax regression (also called multinomial logistic regression).

We have one weight vector $\theta^{(k)} \in \mathbb{R}^d$ per class. The probability of class $k$ is:

$$P(y=k \mid x; \Theta) = \text{softmax}(\theta^{(1)\top}x, \ldots, \theta^{(K)\top}x)k = \frac{e^{\theta^{(k)\top}x}}{\sum{j=1}^{K}e^{\theta^{(j)\top}x}}$$

The softmax function:

• Maps $K$ real-valued scores (logits) to a valid probability distribution summing to 1.
• Amplifies the largest score (temperature controlled by scaling).
• Reduces to sigmoid when $K=2$.

The loss is categorical cross-entropy (negative log-likelihood): $$J(\Theta) = -\frac{1}{n}\sum_{i=1}^n \sum_{k=1}^K \mathbf{1}{y^{(i)}=k}\log P(y^{(i)}=k \mid x^{(i)}; \Theta)$$

This is the exact same loss used in the output layer of classification neural networks — logistic/softmax regression is literally the last layer of a deep net with no hidden layers. Understanding it here means you already understand the output layer of ResNet-152.

def softmax(logits):
    """Numerically stable softmax."""
    logits = logits - logits.max(axis=-1, keepdims=True)  # subtract max
    exp_l = np.exp(logits)
    return exp_l / exp_l.sum(axis=-1, keepdims=True)

# K=3 classes, batch of 4 examples
logits = np.array([[2.0, 1.0, 0.1],
                   [1.0, 3.0, 0.2],
                   [0.1, 0.2, 4.0],
                   [1.5, 1.5, 1.5]])
probs = softmax(logits)
print(probs.sum(axis=1))  # [1. 1. 1. 1.] — valid distributions

Training loss is not the only measure of model quality. For regression and classification, different metrics answer different questions.

For regression:

• MSE / RMSE: same units as $y$; penalizes large errors heavily.
• MAE: mean absolute error; more robust to outliers.
• $R^2$ (coefficient of determination): fraction of variance in $y$ explained by the model. $R^2 = 1 - \frac{\text{SS}{\text{res}}}{\text{SS}{\text{tot}}}$ where $\text{SS}{\text{res}} = \sum(y - \hat{y})^2$ and $\text{SS}{\text{tot}} = \sum(y - \bar{y})^2$. $R^2 = 1$ is perfect; $R^2 = 0$ means the model is no better than predicting the mean; $R^2 < 0$ is possible on test data.

For classification:

• Accuracy: fraction correct. Misleading on imbalanced data (99% negative class → 99% accuracy by predicting all negative).
• Precision / Recall: precision = $\frac{TP}{TP+FP}$ (of positives predicted, how many were correct?); recall = $\frac{TP}{TP+FN}$ (of true positives, how many did we catch?).
• F1 score: harmonic mean of precision and recall.
• ROC-AUC: area under the receiver-operating characteristic curve; measures ranking quality across all thresholds. AUC = 1 is perfect; AUC = 0.5 is random.
• Log-loss (cross-entropy): measures probability calibration — are predicted probabilities actually accurate?

Always separate your training error from test/validation error. Overfitting: low training error, high test error. Underfitting: high training error (and high test error).

from sklearn.metrics import r2_score, roc_auc_score, log_loss

# Regression example
y_true = np.array([3.0, 2.5, 4.0, 5.0])
y_pred = np.array([2.8, 2.7, 3.9, 5.2])
print(f"R2: {r2_score(y_true, y_pred):.4f}")

# Classification example
y_cls = np.array([0, 1, 1, 0, 1])
y_prob = np.array([0.1, 0.9, 0.8, 0.3, 0.7])
print(f"AUC: {roc_auc_score(y_cls, y_prob):.4f}")
print(f"Log-loss: {log_loss(y_cls, y_prob):.4f}")