Week 4: Least Squares and Conditioning

Least Squares as Optimization#

Given a design matrix $X \in \mathbb{R}^{n \times d}$ and labels $y \in \mathbb{R}^n$, ordinary least squares (OLS) minimises:

$f(w) = \frac{1}{2}\|Xw - y\|_2^2 = \frac{1}{2} \sum_{i=1}^n (w^T x_i - y_i)^2$

Expanding (the $\frac{1}{2}$ is a convenient factor that cancels with the derivative of the square):

$f(w) = \frac{1}{2} w^T (X^T X) w - y^T X w + \frac{1}{2} y^T y$

This is a convex quadratic in $w$ — the Hessian is $X^T X$, which is always PSD. The gradient is:

$\nabla f(w) = X^T X w - X^T y = X^T (X w - y)$

GD update: $w_{k+1} = w_k - \eta X^T (X w_k - y)$. Since the Hessian $X^T X$ is constant, $f$ is $L$-smooth with $L = \lambda_{\max}(X^T X)$, the largest eigenvalue of the Gram matrix.

Worked example — $n=3, d=2$:

$X = \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 0 & 2 \end{pmatrix}, \quad y = \begin{pmatrix} 2 \\ 3 \\ 6 \end{pmatrix}$

$X^T X = \begin{pmatrix} 2 & 1 \\ 1 & 5 \end{pmatrix}, \quad X^T y = \begin{pmatrix} 5 \\ 15 \end{pmatrix}$

Eigenvalues of $X^T X$: the characteristic polynomial is $\lambda^2 - 7\lambda + 9$, so $\lambda \in \{\frac{7 \pm \sqrt{13}}{2}\} \approx \{5.30, 1.70\}$. So $L \approx 5.30$, $\mu \approx 1.70$, $\kappa \approx 3.13$ — well-conditioned.

Closed Form vs Iterative#

The OLS solution solves the normal equations:

$X^T X w^* = X^T y \quad \implies \quad w^* = (X^T X)^{-1} X^T y$

When $X^T X$ is invertible (i.e. $X$ has full column rank), this gives the exact minimiser in one linear solve.

So why ever use gradient descent?

1. $n$ is huge: Computing $X^T X$ costs $\mathcal{O}(n d^2)$. If $n = 10^9$, this is prohibitive — but computing $\nabla f_i(w) = (w^T x_i - y_i) x_i$ for one sample is $\mathcal{O}(d)$, which is why SGD (Week 6) shines.
2. $d$ is huge: Inverting $X^T X$ costs $\mathcal{O}(d^3)$. If $d = 10^6$, the closed form is impossible — but a GD step is $\mathcal{O}(nd)$.
3. Regularisation changes the game: Adding terms beyond $\ell_2$ (like $\ell_1$ in Week 10) breaks the closed form entirely.

The tradeoff: the closed form gives exact answers in one shot for small-to-moderate $d$. Iterative methods trade exactness for scalability and flexibility.

Ridge Regression#

Adding $\ell_2$ regularisation gives ridge regression:

$f(w) = \frac{1}{2}\|Xw - y\|_2^2 + \frac{\lambda}{2}\|w\|_2^2$

The gradient becomes $\nabla f(w) = X^T(Xw - y) + \lambda w$, and the Hessian is $X^T X + \lambda I$.

The crucial effect: even if $X^T X$ is nearly singular (some eigenvalues near zero), adding $\lambda I$ shifts all eigenvalues up by $\lambda$:

$\lambda_i(X^T X + \lambda I) = \lambda_i(X^T X) + \lambda$

This guarantees strong convexity with $\mu \geq \lambda > 0$ and improves the condition number:

$\kappa_{\text{ridge}} = \frac{\lambda_{\max} + \lambda}{\lambda_{\min} + \lambda} \leq \kappa_{\text{OLS}}$

Regularisation makes GD converge faster — a happy synergy between statistical regularisation and optimisation efficiency.

The closed form is $w^*_{\text{ridge}} = (X^T X + \lambda I)^{-1} X^T y$.

Condition Number and Convergence#

Recall from Week 3: for strongly convex $f$, GD's convergence factor is $\frac{\kappa - 1}{\kappa + 1}$.

Practical intuition — iterations needed to reach tolerance $\varepsilon$:

When $\kappa$ is large, we look to momentum (reduces effective $\kappa$ to $\sqrt{\kappa}$, Week 7), Newton/quasi-Newton methods (directly use Hessian to undo ill-conditioning), or adaptive methods (per-parameter step sizes, Week 8).

Numerical Example — $\kappa = 1$ vs $\kappa = 100$#

Consider two least squares problems in $\mathbb{R}^2$, both with optimum $w^* = (0, 0)$:

• Well-conditioned: $X^T X = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$, $\kappa = 1$
• Ill-conditioned: $X^T X = \begin{pmatrix} 100 & 0 \\ 0 & 1 \end{pmatrix}$, $\kappa = 100$

GD with $\eta = 1/L$ on the ill-conditioned problem will shrink $w_1$ quickly (large gradient component) but $w_2$ very slowly (small gradient component). The iterates trace a zigzag down the narrow valley.

Knowledge Check#

Browser lab: conditioning and convergence speed#

Run GD on LS with condition numbers $\kappa = 1$, $10$, and $100$. Plot $\|w_k - w^*\|$ vs iteration on a log scale and compare convergence speed.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
d = 2
n = 200

kappas = [1, 10, 100]
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))

for ax, kappa in zip(axes, kappas):
    diag = np.array([np.sqrt(kappa), 1.0])
    w_true = rng.normal(size=d)
    X = rng.normal(size=(n, d)) @ np.diag(diag)
    y = X @ w_true

    XTX = X.T @ X
    L = np.linalg.eigvalsh(XTX).max()
    eta = 1.0 / L

    w = np.zeros(d)
    K = 200
    dists = [np.linalg.norm(w - w_true)]
    for _ in range(K):
        g = X.T @ (X @ w - y)
        w = w - eta * g
        dists.append(np.linalg.norm(w - w_true))

    ax.semilogy(dists, linewidth=1.5)
    ax.set_xlabel("Iteration")
    ax.set_ylabel(r"$\|w_k - w^*\|$")
    ax.set_title(f"$\\kappa \\approx {kappa}$")
    ax.grid(True, alpha=0.3)
    print(f"kappa={kappa:3d}: final distance to w* = {dists[-1]:.3e}")

plt.tight_layout()
plt.show()

print("Notice: larger κ means slower convergence — GD crawls along the shallow direction.")

Further Reading#

• Stephen Boyd & Lieven Vandenberghe, Introduction to Applied Linear Algebra (Cambridge University Press, 2018). Chapters 12–13 derive least squares and its data-fitting uses. Freely available at vmls-book.stanford.edu.
• Lloyd N. Trefethen & David Bau, Numerical Linear Algebra (SIAM, 1997). Lectures 11–12 cover least squares and the conditioning of the normal equations.
• Jorge Nocedal & Stephen Wright, Numerical Optimization (2nd ed., Springer, 2006). Chapter 10 is the reference treatment of linear least-squares problems.