Population risk (what we want to minimize):

Cannot compute without knowing \(p(x,y)\)

Empirical risk (what we can compute):

Empirical Risk Minimization (ERM):

Convergence guarantee (fixed \(\mathcal{F}\), not too rich):

Risk of the sample minimizer approaches the best risk in \(\mathcal{F}\)

Training error: \(\hat{R}_n(\hat{f}_n)\), the empirical risk of the fit on the data that produced it

Test error: \(\hat{R}(\hat{f}_n)\) on held-out data, an unbiased estimate of \(R(\hat{f}_n)\)

The order never changes (\(f^*\) the best in \(\mathcal{F}\)):

• Left: \(\hat{f}_n\) was chosen to make \(\hat{R}_n\) small
• Right: \(f^*\) is the best on the population, \(\hat{f}_n\) is not

Least squares, \(p\) weights, \(n\) samples, noise variance \(\sigma^2\):

• \(\mathbb{E}[\text{training MSE}] = \sigma^2(1 - p/n)\)
• \(\mathbb{E}[\text{MSE on fresh noise at the same inputs}] = \sigma^2(1 + p/n)\)
• \(p = n\): training error 0, the fit interpolates the noise

Overfitting: the gap \(2\sigma^2 p/n\) grows with every added weight

Assume data comes from a process:

where \(\epsilon\) is random noise (\(\mathbf{x}\) includes a constant 1, so the bias is in \(\mathbf{w}\))

What we observe:

• Pairs \((\mathbf{x}_i, y_i)\)
• Never \(\mathbf{w}^T \mathbf{x}_i\) and \(\epsilon_i\) separately

Questions:

1. What distribution for \(\epsilon\)?
2. How does this affect estimation?

Common assumption: \(\epsilon \sim \mathcal{N}(0, \sigma^2)\)

Why Gaussian? Often reasonable in practice:

• Central limit theorem: sum of many small effects
• Maximum entropy for fixed variance
• Mathematical tractability

Model: \(y = \mathbf{w}^T \mathbf{x} + \epsilon\), where \(\epsilon \sim \mathcal{N}(0, \sigma^2)\)

This implies \(y|\mathbf{x}\) is Gaussian:

Density of \(y\) given \(\mathbf{x}\):

For entire dataset \(\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^n\):

Assuming independent samples:

Maximum likelihood principle: Find \(\mathbf{w}\) that makes observed data most likely:

Equivalently, minimize \(-\frac{1}{n}\sum_i \log p(y_i|\mathbf{x}_i; \mathbf{w})\): ERM with loss \(\ell = -\log p\)

Log-likelihood:

Substitute Gaussian pdf:

Simplify:

To maximize \(\log L(\mathbf{w})\):

• First term: constant w.r.t. \(\mathbf{w}\)
• Second term: minimize \(\sum_i(y_i - \mathbf{w}^T \mathbf{x}_i)^2\)

Result: MLE = Least squares

Gaussian noise assumption leads to squared loss

Real data has outliers:

• Sensor errors
• Data entry mistakes
• Rare events

Laplace noise model:

Heavier tails than Gaussian: large residuals are not unusual

Likelihood with Laplace noise:

Log-likelihood:

Maximizing \(\log L(\mathbf{w})\) equivalent to minimizing:

L1 loss emerges directly from the Laplace noise assumption.

This is least absolute deviations (LAD): the fit tracks the conditional median, not the mean

• Not Lasso: Lasso keeps squared loss and adds an L1 penalty on \(\mathbf{w}\)

Notation: Weight vector \(\mathbf{w}\), bias \(b\) (replacing LMMSE parameters \(\mathbf{a}, b\))

Model structure:

where:

• \(\mathbf{x} \in \mathbb{R}^p\): feature vector
• \(\mathbf{w} \in \mathbb{R}^p\): weight vector
• \(b \in \mathbb{R}\): bias/intercept
• \(\hat{y} \in \mathbb{R}\): prediction

Vectorized over dataset:

• \(\mathbf{X} \in \mathbb{R}^{n \times p}\): data matrix
• \(\mathbf{y} \in \mathbb{R}^n\): target vector
• \(\mathbf{1} \in \mathbb{R}^n\): vector of ones

Augmented notation (absorb the bias):

Then: \(\hat{y} = \tilde{\mathbf{w}}^T \tilde{\mathbf{x}}\) and \(\hat{\mathbf{y}} = \tilde{\mathbf{X}}\tilde{\mathbf{w}}\). From here on \(\mathbf{X}\) includes the column of ones.

Empirical risk (squared loss):

Matrix form:

Expanded:

This is a quadratic in \(\mathbf{w}\):

where:

• \(\mathbf{A} = \mathbf{X}^T\mathbf{X}\) (always positive semidefinite)
• \(\mathbf{c} = \mathbf{X}^T\mathbf{y}\)

Convexity: \(\nabla^2 J = \frac{1}{n}\mathbf{X}^T\mathbf{X} \succeq 0\)

• \(\succeq 0\): positive semidefinite, \(\mathbf{v}^T (\nabla^2 J) \mathbf{v} \geq 0\) for every \(\mathbf{v}\)
• Curvature never negative in any direction, so every local minimum is global

Minimize: \(J(\mathbf{w}) = \frac{1}{2n}\|\mathbf{y} - \mathbf{X}\mathbf{w}\|^2\)

Take gradient:

Set to zero:

Normal equations:

Observations:

1. \(\mathbf{X}^T\mathbf{X} \in \mathbb{R}^{p \times p}\): Gram matrix
2. \(\mathbf{X}^T\mathbf{y} \in \mathbb{R}^p\): cross-correlation
3. System is \(p\) equations in \(p\) unknowns
4. Unique solution if \(\mathbf{X}^T\mathbf{X}\) invertible

Hence name: residual \(\mathbf{y} - \mathbf{X}\hat{\mathbf{w}}\) is orthogonal (normal) to every column of \(\mathbf{X}\) - the sample form of the orthogonality principle

Normal equations: \(\mathbf{X}^T\mathbf{X}\mathbf{w} = \mathbf{X}^T\mathbf{y}\)

Divide by \(n\):

With the column of ones: \(\frac{1}{n}\mathbf{X}^T\mathbf{X}\) holds second moments, not covariances

Center first (\(\mathbf{X}_c\): each column minus its mean, \(\mathbf{y}_c = \mathbf{y} - \bar{y}\)): the ones column drops out and the equations become

• \(\hat{\mathbf{K}}_{XX} = \frac{1}{n}\mathbf{X}_c^T\mathbf{X}_c\): sample covariance matrix
• \(\hat{\mathbf{k}}_{XY} = \frac{1}{n}\mathbf{X}_c^T\mathbf{y}_c\): sample cross-covariance
• \(b = \bar{y} - \mathbf{w}^T\bar{\mathbf{x}}\) recovers the intercept

Solution:

This is empirical LMMSE.

Population: \(\mathbf{w}^* = \mathbf{K}_{XX}^{-1} \mathbf{k}_{XY}\) Sample: \(\hat{\mathbf{w}} = \hat{\mathbf{K}}_{XX}^{-1} \hat{\mathbf{k}}_{XY}\)

Consistency: \(\hat{\mathbf{w}} \to \mathbf{w}^*\) as \(n \to \infty\)

Solution (when \(\mathbf{X}^T\mathbf{X}\) invertible):

Moore-Penrose pseudoinverse (full column rank):

So: \(\mathbf{w} = \mathbf{X}^+\mathbf{y}\)

Properties of \(\mathbf{X}^+\):

1. \(\mathbf{X}^+\mathbf{X} = \mathbf{I}_p\) (left inverse)
2. \(\mathbf{X}\mathbf{X}^+ = \mathbf{P}_X\) (projection onto col(\(\mathbf{X}\)))

When is \(\mathbf{X}^T\mathbf{X}\) invertible?

• Need \(\text{rank}(\mathbf{X}) = p\) (full column rank)
• Requires \(n \geq p\) (more data than features)
• Columns of \(\mathbf{X}\) linearly independent

Rank deficient case: infinitely many solutions; the SVD pseudoinverse picks the one with minimum \(\|\mathbf{w}\|\)

Method 1: Normal equations

w = np.linalg.solve(X.T @ X, X.T @ y)

Complexity:

• Form \(\mathbf{X}^T\mathbf{X}\): \(O(np^2)\)
• Form \(\mathbf{X}^T\mathbf{y}\): \(O(np)\)
• Solve system: \(O(p^3)\)
• Total: \(O(np^2 + p^3)\), about \(np^2 + p^3/3\) flops

Condition number: \(\kappa(\mathbf{A}) = \sigma_{\max}(\mathbf{A}) / \sigma_{\min}(\mathbf{A})\), the ratio of largest to smallest singular value

• The largest factor by which solving with \(\mathbf{A}\) amplifies a relative error in the data
• For \(\mathbf{X}^T\mathbf{X}\): \(\lambda_{\max} / \lambda_{\min}\) of its eigenvalues (lecture 03’s \(\mathbf{K}_{XX}\))

Squaring:

• Squaring doubles the digits lost
• float64 carries ~16 digits: at \(\kappa(\mathbf{X}) = 10^8\), none survive

Singular Value Decomposition:

where:

• \(\mathbf{U} \in \mathbb{R}^{n \times p}\): orthonormal columns (left singular vectors)
• \(\boldsymbol{\Sigma} \in \mathbb{R}^{p \times p}\): diagonal, \(\sigma_1 \geq \cdots \geq \sigma_p \geq 0\)
• \(\mathbf{V} \in \mathbb{R}^{p \times p}\): orthogonal (right singular vectors)

Solution via SVD:

where \(\boldsymbol{\Sigma}^+\) = pseudoinverse of \(\boldsymbol{\Sigma}\):

• If \(\sigma_i > \delta\): \(\sigma_i^+ = 1/\sigma_i\)
• If \(\sigma_i \leq \delta\): \(\sigma_i^+ = 0\)

Rank deficient: smallest \(\|\mathbf{w}\|\) among all minimizers

Component form: one term per singular direction

• Small \(\sigma_i\) divides: that direction carries the noise
• The filter factor \(f_i\) keeps or drops each direction whole

Cost:

• About \(2np^2 + 11p^3\) flops, the most of the three
• np.linalg.lstsq uses the SVD

Truncated SVD (regularization): \(f_i = 1\) for the \(k\) largest \(\sigma_i\), else \(0\)

• Ridge replaces the 0/1 switch by a smooth factor between them

First-order Taylor expansion:

• For a fixed step length, the decrease is largest at \(\Delta \propto -\nabla J(\mathbf{w})\)

Gradient points uphill:

Direction of maximum increase of \(J\)

Update rule:

where \(\eta\) = step size / learning rate

Convergence criterion:

At optimum: \(\nabla J(\mathbf{w}^*) = \mathbf{0}\) (critical point)

NumPy (\(\mathbf{X}\) includes the ones column):

def gradient_descent(X, y, eta=0.01, steps=1000):
    n, p = X.shape
    w = np.zeros(p)
    for _ in range(steps):
        grad = X.T @ (X @ w - y) / n
        w = w - eta * grad
    return w

Definition: \(f\) is convex if for all \(\lambda \in [0,1]\):

Squared loss is convex:

Proof via Hessian:

Since \(\mathbf{X}^T\mathbf{X} \succeq 0\) (positive semidefinite):

Therefore \(J\) is convex.

1. Any local minimum is global minimum
2. A minimizer always exists; unique if \(\mathbf{X}^T\mathbf{X} \succ 0\)
3. Gradient descent converges to global minimum for a small enough step size
4. No need for random restarts

Strong convexity if \(\mathbf{X}\) full column rank:

where \(\mu = \lambda_{\min}(\mathbf{X}^T\mathbf{X})/n > 0\)

Error recursion (quadratic \(J\), Hessian \(\mathbf{H} = \frac{1}{n}\mathbf{X}^T\mathbf{X}\)):

• Along eigenvector \(i\) of \(\mathbf{H}\), the error scales by \(1 - \eta\lambda_i\) each step

Too small: Slow convergence

Too large: Overshoot, divergence

Theoretical bounds: need \(|1 - \eta\lambda_i| < 1\) for every \(i\):

Optimal step size:

Convergence rate (\(\kappa = \lambda_{\max}/\lambda_{\min}\)):

Scale mismatch inflates \(\kappa(\mathbf{H})\) (figure):

• Raw features: \(\kappa \approx 10^5\)
• Standardized: \(\kappa \approx 1\)

Rescaling:

• Changes gradient descent speed and the ridge penalty
• Leaves least-squares predictions unchanged

Standardization (z-score): \(x' = (x - \mu)/\sigma\)

• Zero mean, unit variance (the default choice)
• Preserves outliers

Min-Max scaling: \(x' = (x - x_{min})/(x_{max} - x_{min})\)

• Maps to [0, 1]
• Sensitive to outliers
• Bounded output
• Does not center: \(\kappa \approx 260\) in the figure

Robust scaling: \(x' = (x - \text{median})/\text{IQR}\)

• Uses median and interquartile range
• Robust to outliers

Use one sample at a time:

where \(i_t\) randomly selected

For squared loss (\(\ell_i = \frac{1}{2}(y_i - \mathbf{w}^T\mathbf{x}_i)^2\)):

LMS algorithm (Widrow-Hoff) is SGD for squared loss: \(\mathbf{w}_{t+1} = \mathbf{w}_t + \eta e_t \mathbf{x}_t\)

Unbiased gradient estimate:

Noise in gradient:

Cost and use:

• Cheap iteration: \(O(p)\) vs \(O(np)\)
• Online learning capability

Fixed step size: Converges to neighborhood

where \(\sigma_g^2\) = gradient noise variance, \(\mu = \lambda_{\min}(\mathbf{H})\)

Decreasing step size: Converges to exact solution

Required conditions (Robbins-Monro):

Example: \(\eta_t = \eta_0/t\)

Lowering the floor:

• Mini-batching: gradient variance \(\sigma_g^2/m\)
• Decreasing \(\eta_t\)
• Momentum, Adam, RMSprop change the step, not the gradient variance

Convergence rates:

Mini-batch gradient:

where \(|\mathcal{B}| = m\) = batch size

Variance reduction:

Consequences:

1. Gradient noise (std) falls by \(\sqrt{m}\)
2. The noise floor \(O(\eta\sigma_g^2/(m\mu))\) allows a larger \(\eta\) at the same floor
3. The \(m\) gradients compute in parallel (one pass on a GPU)

Epoch: One pass through dataset

• Batch: 1 update per epoch
• SGD: \(n\) updates per epoch
• Mini-batch: \(n/m\) updates per epoch

Equal compute:

• One step with \(m\) samples: variance \(\sigma_g^2/m\)
• \(m\) single-sample steps: same samples, \(m\) updates
• The gain from batching is parallel hardware, not fewer samples

From \(\hat{\mathbf{w}} - \mathbf{w}^* = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\epsilon\):

• Large eigenvalues of \(\mathbf{X}^T \mathbf{X}\): precise
• Small eigenvalues: uncertain along those directions

Fisher information (curvature of the log-likelihood):

• \(\text{Cov}(\hat{\mathbf{w}}) = \mathcal{I}^{-1}\): sharper peak, smaller variance

Confidence region: the ellipsoid \((\hat{\mathbf{w}} - \mathbf{w}^*)^T \mathcal{I} (\hat{\mathbf{w}} - \mathbf{w}^*) \leq c\)

• Axes along the eigenvectors of \(\mathbf{X}^T\mathbf{X}\), longest where its eigenvalue is smallest

Cramér-Rao bound: any unbiased estimator \(\tilde{\mathbf{w}}\) has \(\text{Cov}(\tilde{\mathbf{w}}) \succeq \mathcal{I}^{-1}\)

• Least squares attains it for every \(n\)
• Lower variance requires bias

Setup: \(y = f(x) + \epsilon\), \(\;\mathbb{E}[\epsilon] = 0\), \(\;\text{Var}(\epsilon) = \sigma^2\)

Expand at point \(x_0\), with \(\hat{f}\) random (trained on random data):

Independence of \(\epsilon\) and \(\hat{f}\):

Add and subtract \(\mathbb{E}[\hat{f}(x_0)]\):

For linear models, averaged over the training inputs (what the figure measures):

• Variance: \(\frac{\sigma^2}{n}\text{tr}\left(\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\right) = p\sigma^2/n\)
• Bias decreases as \(p\) increases (nested models)

Measured (polynomial degree, \(n = 20\), \(\sigma = 0.1\), the figure’s setup):

Overfitting:

• Variance dominates bias
• Training error \(\ll\) test error