ML2.4 · LESSON 4 OF 6

Closed-form fitting

Fit a linear predictor by hand, implement its loss and gradient, then compare with a library solution.

THE IDEA

Skip the walk: solve for the bottom directly

Gradient descent walks down the bowl. For linear regression the bottom can be written down. Set the gradient to zero and the equations are linear in w and b, so a linear solver finds the answer in one go.

Step 1 of 4

Gradient equal to zero. At the bottom of the bowl the gradient is 0. Writing the model as X·w with a column of x values and a column of ones, that condition is Xᵀ(Xw − y) = 0.

points used = 4 every x the same (1.5)

X^T X = [[ 33.14, -0.60 ], X^T y = [ 67.08 ] [ -0.60, 4 ]] [ 3.90 ] (first column x, second column ones) det = 4 * 33.14 - (-0.60)^2 = 132.200

w = (n*Sxy - Sx*Sy) / det = 2.0474 b = (Sxx*Sy - Sx*Sxy) / det = 1.2821

gradient descent iterations = 0

	w	b
closed form	2.0474	1.2821
gradient descent after 0	0.0000	0.0000
gap	-2.0474	-1.2821

Gradient descent uses learning rate 0.0603 here (half of the largest stable rate), starting from w = 0, b = 0. Green solid: closed form. Orange dashed: gradient descent.

Choose how many points to use and read the live matrices. The green line is the exact closed-form answer. The orange dashed line is gradient descent after the iterations you choose.

Closed form versus gradient descent

The closed form gives the exact least-squares line in one solve. Gradient descent approaches it, getting closer with each iteration and never quite arriving. With 12 points and a safe learning rate, about 30 iterations get both w and b within 0.01 of the exact answer. With two points the fitted line passes through both: two points determine a line, so the loss is zero. For models with millions of parameters the closed form is too expensive and most models have no closed form at all, which is why gradient descent matters. For small linear problems, the solve wins.

Predict first. Two features are exact copies of each other. Is XᵀX still invertible?

In code: solve, do not invert

The textbook formula and what reliable code does

# textbook:  w = (X^T X)^-1 X^T y        # forms an inverse: slower, less accurate
w = np.linalg.solve(X.T @ X, X.T @ y)   # solves the system directly
w, *_ = np.linalg.lstsq(X, y, rcond=None)   # QR/SVD based, also handles rank-deficient X

Chapter context

Linear regression fits a weighted sum of features. The model is simple enough to inspect numerically, making it a useful bridge between geometry, optimization and software implementation.

TRY IT YOURSELF

Part 1. Do it with the visual

The points are added in a fixed, well-spread order. The gradient-descent line uses a safe fixed learning rate, so only the iteration count changes.

points used = 12 every x the same (1.5)

X^T X = [[ 76.56, -0.60 ], X^T y = [ 155.07 ] [ -0.60, 12 ]] [ 9.80 ] (first column x, second column ones) det = 12 * 76.56 - (-0.60)^2 = 918.360

w = (n*Sxy - Sx*Sy) / det = 2.0327 b = (Sxx*Sy - Sx*Sxy) / det = 0.9183

gradient descent iterations = 0

	w	b
closed form	2.0327	0.9183
gradient descent after 0	0.0000	0.0000
gap	-2.0327	-0.9183

Gradient descent uses learning rate 0.0784 here (half of the largest stable rate), starting from w = 0, b = 0. Green solid: closed form. Orange dashed: gradient descent.

To do: On all 12 points, run gradient descent for enough iterations that both w and b are within 0.01 of the closed-form answer
To do: Use only 2 points: the closed-form line must pass exactly through both, so its MSE is 0
To do: Break the method: make every house the same size so that X^T X is singular and there is no single best line

Part 2. Do it in code

design(x) returns the n by 2 matrix whose first column is x and whose second column is all ones.
fit_normal(X, y) solves the normal equations X^T X w = X^T y with np.linalg.solve (not an inverse). With dependent columns it has no usable answer (an error or nan).
fit_lstsq(X, y) uses np.linalg.lstsq and returns just the weights. fit(X, y) uses fit_normal when X has full column rank (np.linalg.matrix_rank) and falls back to fit_lstsq otherwise.

import numpy as np

def design(x):
    # columns: x, then ones
    pass

def fit_normal(X, y):
    # solve X^T X w = X^T y
    pass

def fit_lstsq(X, y):
    # same fit through np.linalg.lstsq (return only the weights)
    pass

def fit(X, y):
    # normal equations if possible, least squares otherwise
    pass

print(design(np.array([1.0, 2.0, 3.0])))

Ctrl/⌘ + Enter runs it

Stuck, or something not working? You can come back to it any time.

CHECK YOUR UNDERSTANDING

Does a duplicated feature guarantee an invertible X-transpose-X?

ANoBYes

Chapter reference: math, worked example & code

This section belongs to the whole chapter, so it is the same on every lesson in it and may cover a different lesson's topic. The code for the experiment you just ran is inside the experiment itself.

Math and notation

\hat y_i=x_i^Tw+b,\quad L=\frac1n\sum_i(\hat y_i-y_i)^2.

\nabla_w L=\frac2nX^T(Xw+b\mathbf1-y),\quad \partial_bL=\frac2n\sum_i(\hat y_i-y_i).

When invertible, a no-intercept least-squares solution is $(X^TX)^{-1}X^Ty$ . Use a pseudoinverse or suitable solver otherwise.

Worked example

Use x=[1,2], y=[2,4], w=1 and b fixed at zero. Predictions [1,2] give residuals [-1,-2], so MSE=(1+4)/2=2.5. The w gradient is 2/2×(−1×1−2×2)=−5. With eta=0.1, w becomes 1.5 and MSE falls to (0.25+1)/2=0.625.

Reference code

The following example uses the Python standard library and does not call a remote service.

x, y, w = [1,2], [2,4], 1.0
mse = lambda weight: sum((weight*a-b)**2 for a,b in zip(x,y))/len(x)
grad = 2/len(x)*sum((w*a-b)*a for a,b in zip(x,y))
assert mse(w) == 2.5 and grad == -5
w -= .1*grad
assert w == 1.5 and mse(w) == .625

Python examples are provided for study and local execution; this page does not execute Python.

Common mistakes

A low residual does not show that a weight is causal. The closed form is not always a literal matrix inverse in reliable code. Include or exclude bias consistently in derivations.