Part 1: Dataset and Problem Setup

• The synthetic dataset has $n=40$ examples and $d=3$ features, yielding an overdetermined system ($n > d$).
• The design matrix $X \in \mathbb{R}^{40 \times 3}$ has full column rank (with high probability for random data), so the least-squares solution is unique.
• The target vector $y \in \mathbb{R}^{40}$ is generated from true coefficients plus noise, ensuring a noisy regression problem (residuals won’t be zero).
• The goal is to find $\hat{w} \in \mathbb{R}^3$ that minimizes $\|y - Xw\|_2^2 = \sum_{i=1}^{40} (y_i - \mathbf{x}_i^\top w)^2$ (sum of squared prediction errors).
• Both methods (lstsq and normal equations) should find the same $\hat{w}$ for well-conditioned $X$, demonstrating that algorithm choice doesn’t affect the solution when stability is not an issue.

Part 2: Two Paths to the Same Solution

• Method 1: SVD-based lstsq
  Internally computes $X = U \Sigma V^\top$ and solves $\hat{w}_1 = V \Sigma^{-1} U^\top y$. The SVD factors $X$ directly (no matrix-matrix multiplication), avoiding condition number squaring. The rcond=None parameter sets the rank cutoff to machine precision: singular values below rcond * max(σ) are treated as zero. If $X$ is rank-deficient, lstsq returns the minimum-norm solution $\hat{w} = V \Sigma^+ U^\top y$ (pseudoinverse via SVD). Complexity: $O(nd^2)$ for thin SVD (assuming $n > d$), though constant factors are larger than Cholesky/LU.

• Method 2: Normal equations via direct solve
  Forms Gram matrix $X^\top X \in \mathbb{R}^{3 \times 3}$ explicitly (two matrix multiplications: $O(nd^2)$ flops), then solves $(X^\top X) \hat{w}_2 = X^\top y$ via LU or Cholesky decomposition. The solve function assumes $X^\top X$ is invertible; if nearly singular, it may return garbage weights without warning. Complexity: $O(nd^2 + d^3)$ (form Gram matrix + solve $d \times d$ system). For small $d$ (< 100), this is faster than SVD.

• Why both work here: For well-conditioned $X$ with full rank, both methods converge to the same least-squares solution (the unique minimizer of $\|y - Xw\|_2^2$). The normal equations are mathematically equivalent to the SVD solution, and lstsq implements the same minimization via stable factorization.

Part 3: Verifying Solution Agreement

• Weight difference $\|\hat{w}_1 - \hat{w}_2\|_2$: Measures Euclidean distance between the two weight vectors. For well-conditioned problems, this should be near machine precision ($\sim 10^{-15}$), confirming both methods found the same solution. If this exceeds $10^{-10}$, $X$ is ill-conditioned and numerical errors have diverged the solutions.

• Residual norms $\|r_1\|_2$ and $\|r_2\|_2$: Measure the magnitude of prediction errors. Since both $\hat{w}_1$ and $\hat{w}_2$ minimize the same objective $\|y - Xw\|_2^2$, their residual norms should be identical (or differ only by numerical error). The residual norm is $\|r\|_2 = \sqrt{\sum_{i=1}^{40} (y_i - \mathbf{x}_i^\top \hat{w})^2}$, the square root of the minimized loss.

• Typical output:

  ||w1-w2||: 1.234e-15
  ||r1||, ||r2||: 2.456 2.456

  The weight difference is at machine precision (effectively zero), and residual norms are identical, confirming both methods found the same least-squares solution.

• Orthogonality condition: Both solutions should satisfy $X^\top r \approx 0$ (residuals orthogonal to column space). Verify via np.linalg.norm(X.T @ r1) and np.linalg.norm(X.T @ r2), which should be near zero (typically $\sim 10^{-14}$). This is the first-order optimality condition for least-squares.

Part 4: When Normal Equations Fail

• If $X$ has condition number $\kappa(X) = 10^6$ (moderately ill-conditioned), then $X^\top X$ has condition number $\kappa(X^\top X) = 10^{12}$ (severely ill-conditioned).
• Numerical errors in forming $X^\top X$ (matrix multiplication introduces $O(\epsilon_{\text{machine}} \cdot \kappa(X)^2)$ relative error) and solving the system (LU decomposition amplifies errors by $\kappa(X^\top X)$) can produce wildly incorrect weights.
• The lstsq approach avoids this by working directly with $X$ (condition number $10^6$, not $10^{12}$). SVD computes singular values accurately, and the solution $\hat{w} = V \Sigma^{-1} U^\top y$ has relative error bounded by $O(\epsilon_{\text{machine}} \cdot \kappa(X))$, not $\kappa(X)^2$.
• Silent failure: If $X^\top X$ is nearly singular, np.linalg.solve may return garbage weights with no warning (only raises LinAlgError if exactly singular). The lstsq approach issues warnings for rank-deficient matrices and returns the minimum-norm solution.

Why Gram Matrix Structure Matters in ML

• Gram matrix $X^\top X$ is always symmetric: $(X^\top X)^\top = X^\top X$. This enables faster solvers (Cholesky is 2x faster than LU for symmetric PD matrices).
• Gram matrix is PSD: For any $v \in \mathbb{R}^d$, $v^\top (X^\top X) v = \|Xv\|_2^2 \ge 0$. It’s PD (strictly positive) if and only if $X$ has full column rank.
• Eigenvalues of $X^\top X$ are squared singular values: If $X = U \Sigma V^\top$, then $X^\top X = V \Sigma^2 V^\top$. Eigenvalues $\lambda_i = \sigma_i^2$ are nonnegative.
• Ridge regression stabilizes normal equations: Adding $\lambda I$ to $X^\top X$ ensures positive definiteness: $(X^\top X + \lambda I) \hat{w}_\lambda = X^\top y$. Even small $\lambda$ (e.g., $10^{-6}$) dramatically improves conditioning for nearly singular $X^\top X$.

ML Examples and Patterns

• Ridge regression comparison: Ridge solves $(X^\top X + \lambda I) \hat{w}_\lambda = X^\top y$. Normal equations approach:

  w_ridge = np.linalg.solve(X.T @ X + lambda_reg * np.eye(d), X.T @ y)

  But lstsq doesn’t directly support ridge. Use augmented least-squares instead:

  X_aug = np.vstack([X, np.sqrt(lambda_reg) * np.eye(d)])
  y_aug = np.hstack([y, np.zeros(d)])
  w_ridge = np.linalg.lstsq(X_aug, y_aug, rcond=None)[0]

  This avoids forming $X^\top X$ explicitly, maintaining numerical stability.

• Polynomial regression with high degree: Fitting $y = w_0 + w_1 x + w_2 x^2 + \ldots + w_p x^p$ creates design matrix with columns $[1, x, x^2, \ldots, x^p]$. These become highly collinear for large $p$ (Vandermonde matrix has exponentially large condition number). Normal equations fail catastrophically ($p > 10$ often breaks); lstsq handles $p \sim 20$ gracefully.

• Cholesky vs. LU decomposition: If $X$ has full rank, $X^\top X$ is PD. Cholesky factorization $X^\top X = LL^\top$ is faster and more stable than LU:

  L = np.linalg.cholesky(X.T @ X)  # L L^T = X^T X
  z = np.linalg.solve(L, X.T @ y)  # Forward substitution
  w_chol = np.linalg.solve(L.T, z)  # Backward substitution

  Roughly 2x faster than solve for symmetric systems, but requires verifying positive definiteness (raises error if $X^\top X$ is singular or indefinite).

• Large-scale regression ($n \gg d$): For many examples and few features, normal equations are efficient: form $X^\top X$ in $O(nd^2)$ time, solve $d \times d$ system in $O(d^3)$. For $d < 100$, this is faster than SVD. For $d \gg n$ (more features than examples), use dual formulation or kernel methods.

• Iterative refinement for ill-conditioned problems: Solve normal equations in low precision, then refine via lstsq:

  w_approx = np.linalg.solve(X.T @ X, X.T @ y)  # Initial solve
  for _ in range(5):
      r = y - X @ w_approx  # Residual in double precision
      delta_w = np.linalg.lstsq(X, r, rcond=None)[0]  # Refine via lstsq
      w_approx += delta_w

  Modern libraries (LAPACK, Intel MKL) implement this automatically for lstsq.

Connection to Linear Algebra Theory

• Normal equations derivation: Minimize $\mathcal{L}(w) = \|y - Xw\|_2^2 = (y - Xw)^\top (y - Xw)$. Gradient: \[\frac{\partial \mathcal{L}}{\partial w} = -2X^\top (y - Xw) = 0 \quad \Rightarrow \quad X^\top X w = X^\top y.\] If $X$ has full column rank, $X^\top X$ is invertible: $\hat{w} = (X^\top X)^{-1} X^\top y$.

• SVD solution: For $X = U \Sigma V^\top$, the least-squares solution is $\hat{w} = V \Sigma^{-1} U^\top y$. If $X$ is rank-deficient, use the pseudoinverse $\Sigma^+ = \text{diag}(\sigma_1^{-1}, \ldots, \sigma_r^{-1}, 0, \ldots, 0)$ (invert only nonzero singular values).

• Orthogonality condition: Both solutions satisfy $X^\top r = 0$ (residuals orthogonal to column space). This is the first-order optimality condition for least-squares. Geometrically, $r$ is perpendicular to all columns of $X$.

• Projection matrix: The least-squares solution projects $y$ onto $\text{col}(X)$: \[\hat{y} = X\hat{w} = X(X^\top X)^{-1} X^\top y = P y,\] where $P = X(X^\top X)^{-1} X^\top$ is the orthogonal projection matrix. Residual $r = y - \hat{y} = (I - P)y$ is orthogonal to $\text{col}(X)$.

• Condition number bounds: For perturbation $\Delta y$ in $y$, the relative error in $\hat{w}$ is bounded by: \[\frac{\|\Delta \hat{w}\|}{\|\hat{w}\|} \le \kappa(X) \frac{\|\Delta y\|}{\|y\|}.\] For normal equations with perturbation $\Delta (X^\top y)$, the bound uses $\kappa(X^\top X) = [\kappa(X)]^2$: \[\frac{\|\Delta \hat{w}\|}{\|\hat{w}\|} \le \kappa(X)^2 \frac{\|\Delta (X^\top y)\|}{\|X^\top y\|}.\] This explains the condition number squaring effect.

Numerical and Implementation Notes

• Always check condition number before using normal equations: Compute $\kappa(X)$ via np.linalg.cond(X). If $\kappa(X) > 10^8$, normal equations are unreliable; use lstsq instead.

• Verify symmetry exploitation: NumPy’s solve doesn’t automatically detect symmetry. For large $d$, use scipy.linalg.solve(..., assume_a='pos') to enable Cholesky (faster for PD matrices).

• rcond parameter tuning: rcond=None uses machine precision. For noisy data, a larger rcond (e.g., 1e-10) can improve stability by treating small singular values as zero (truncated SVD).

• QR factorization as middle ground: np.linalg.qr followed by triangular solve is faster than SVD (roughly 2x) and more stable than normal equations. Not exposed in lstsq, but available via:

  Q, R = np.linalg.qr(X)
  w_qr = np.linalg.solve(R, Q.T @ y)

  Condition number of $R$ equals $\kappa(X)$ (no squaring), making QR a good compromise.

Numerical and Shape Notes

• $X \in \mathbb{R}^{40 \times 3}$: design matrix (40 examples, 3 features)
• $y \in \mathbb{R}^{40}$: target vector
• $\hat{w}_1, \hat{w}_2 \in \mathbb{R}^3$: weight vectors from both methods
• $X^\top X \in \mathbb{R}^{3 \times 3}$: Gram matrix (symmetric, square)
• $X^\top y \in \mathbb{R}^3$: projected targets
• $r_1, r_2 \in \mathbb{R}^{40}$: residual vectors
• X.T @ X computes $(3 \times 40) \cdot (40 \times 3) = (3 \times 3)$ Gram matrix
• X.T @ y computes $(3 \times 40) \cdot (40,) = (3,)$ projected targets
• X @ w1 computes $(40 \times 3) \cdot (3,) = (40,)$ predictions

Gotchas:

1. Silent failure of normal equations: If $X^\top X$ is singular or nearly singular, solve may return garbage weights with no warning. Always verify: np.allclose(X.T @ X @ w2, X.T @ y).

2. Condition number not checked: The code doesn’t compute $\kappa(X)$ explicitly. For production use:

   cond_X = np.linalg.cond(X)
   if cond_X > 1e10:
       print(f"Warning: X is ill-conditioned (cond={cond_X:.2e})")

3. Symmetry of $X^\top X$: Always symmetric: $(X^\top X)^\top = X^\top X$. Exploit this for faster solvers (Cholesky).

4. Residual norms equality: Only holds if both methods found the global minimizer. If normal equations fail, $\|r_2\|$ can be arbitrarily large.

Verification checks:

# Verify condition number squaring
cond_X = np.linalg.cond(X)
cond_XTX = np.linalg.cond(X.T @ X)
print(f"cond(X): {cond_X:.2e}, cond(X^T X): {cond_XTX:.2e}")
print(f"Ratio (should be ~cond(X)^2): {cond_XTX / cond_X**2:.4f}")

# Verify normal equations optimality condition
print(f"||X^T r1||: {np.linalg.norm(X.T @ r1):.2e}")  # Should be ~0
print(f"||X^T r2||: {np.linalg.norm(X.T @ r2):.2e}")  # Should be ~0

# Check if solutions are identical (within tolerance)
assert np.allclose(w1, w2, atol=1e-10), "Solutions differ significantly!"

# Verify both minimize the same objective
loss1 = np.linalg.norm(r1)**2
loss2 = np.linalg.norm(r2)**2
print(f"Loss (lstsq): {loss1:.10f}")
print(f"Loss (normal eqs): {loss2:.10f}")
assert abs(loss1 - loss2) < 1e-10, "Losses differ!"

# Check rank (should be 3 for full column rank)
rank_X = np.linalg.matrix_rank(X)
print(f"Rank(X): {rank_X}/3")

Pedagogical Significance

This example is the canonical comparison of least-squares solvers, demonstrating that algorithm choice matters even when mathematical equivalence holds. Two paths to the same solution have radically different numerical properties.

Key takeaways: 1. lstsq and normal equations are mathematically equivalent but numerically different: Both minimize $\|y - Xw\|_2^2$, but lstsq avoids condition number squaring. 2. Condition number squaring is unavoidable when forming $X^\top X$: SVD avoids this by working with $X$ directly. 3. Normal equations can fail silently: Garbage weights with no warning for ill-conditioned $X$; lstsq issues rank warnings. 4. Solver hierarchy reflects stability-speed trade-off: SVD is slowest but most robust; normal equations are fastest but fragile. 5. Residual orthogonality $X^\top r = 0$ characterizes optimality: Unifying condition for both methods.

Common misconceptions addressed: - “Normal equations are always fine”: False—they fail silently for ill-conditioned $X$, producing garbage weights. - “lstsq is overkill for simple regression”: False—it’s the recommended default precisely because it handles edge cases gracefully. - “Faster means better”: False—normal equations are faster but fragile; SVD is slightly slower but robust. Stability favors lstsq for ML. - “Both methods give the same answer”: True for well-conditioned problems, false for ill-conditioned—normal equations diverge, lstsq remains accurate.

Connection to other chapters: - Chapter 10 (SVD): SVD directly yields least-squares solution; singular values reveal conditioning. - Chapter 12 (Least-squares): This is the canonical demonstration of solver choice. - Chapter 14 (Conditioning): Condition number squaring explains why normal equations are risky. - Chapter 15 (Cholesky/LU): Cholesky is faster for normal equations but assumes $X^\top X$ is PD.

Why this snippet is pedagogically powerful: It proves that two mathematically equivalent formulations have radically different numerical properties. Understanding this distinction is fundamental to reliable ML engineering. The code is minimal (10 lines), but the lesson is profound: always prefer lstsq unless you’ve explicitly verified that normal equations are safe for your data.