Demo: Linear Models and Least Squares
View ordinary least squares as projection onto a one-dimensional column space, then compare it with the geometry induced by weighted least squares.
Mathematical setup
The toy model is $y=h\theta+e$ with
\[h= \begin{bmatrix} 1\\ \text{slope} \end{bmatrix}.\]Ordinary least squares solves
\[\hat\theta_{\mathrm{OLS}}=(h^Th)^{-1}h^Ty, \qquad h^T(y-h\hat\theta_{\mathrm{OLS}})=0.\]Weighted least squares with $W=\operatorname{diag}(1,w_2)$ solves
\[\hat\theta_{\mathrm{WLS}}=(h^TWh)^{-1}h^TWy, \qquad h^TW(y-h\hat\theta_{\mathrm{WLS}})=0.\]The weighted residual is orthogonal in the $W$-inner product, so it need not be Euclidean-orthogonal to the column space.
What to try
- Set $w_2=1$. OLS and WLS coincide because the two geometries are the same.
- Increase $w_2$. The fitted point moves to respect the second coordinate more strongly.
- Change the slope. The column space rotates, and both projections change even when the observed $y$ is fixed.
For weighted least squares, and for GLS with $W=\Sigma^{-1}$, the normal equation becomes $H^T W(y-H\hat\theta)=0$. The unweighted residual dot product need not be zero for the weighted fit.
Try it in Python
This cell computes the OLS and weighted projections, then checks the two normal equations directly.
import numpy as np
import matplotlib.pyplot as plt
y = np.array([2.0, 1.0])
slope = 0.6
w2 = 1.0
h = np.array([1.0, slope])
W = np.diag([1.0, w2])
theta_ols = h @ y / (h @ h)
fit_ols = h * theta_ols
resid_ols = y - fit_ols
theta_wls = (h @ W @ y) / (h @ W @ h)
fit_wls = h * theta_wls
resid_wls = y - fit_wls
line_t = np.linspace(-3, 3, 100)
line = np.outer(line_t, h)
plt.plot(line[:, 0], line[:, 1], color="gray", label="column space")
plt.scatter(*y, color="black", label="observed y")
plt.scatter(*fit_ols, label="OLS fit")
plt.scatter(*fit_wls, label="WLS fit")
plt.plot([y[0], fit_ols[0]], [y[1], fit_ols[1]], "--", label="OLS residual")
plt.plot([y[0], fit_wls[0]], [y[1], fit_wls[1]], ":", label="WLS residual")
plt.axis("equal")
plt.xlabel("coordinate 1")
plt.ylabel("coordinate 2")
plt.legend()
plt.show()
print(f"theta_OLS = {theta_ols:.3f}")
print(f"theta_WLS = {theta_wls:.3f}")
print(f"h dot OLS residual = {h @ resid_ols:.3e}")
print(f"h^T W WLS residual = {h @ W @ resid_wls:.3e}")
print(f"h dot WLS residual = {h @ resid_wls:.3e}")