Linear Regression by Normal Equation

Solve linear regression in one shot instead of iterating.

Setting the gradient of the squared error to zero and solving gives the normal equation:

theta = (X_b^T X_b)^-1 X_b^T y

where X_b is X with a leading column of ones so the intercept is fitted as just another weight.

Return (w, b) with w of shape (d,) and b a float, so the result is interchangeable with the gradient descent version.

X_b^T X_b is singular whenever features are collinear or you have fewer rows than columns, and inverting it then blows up. Use np.linalg.pinv, the pseudo-inverse, which returns the minimum norm solution instead of failing.

This is exact and needs no learning rate, but it costs O(d^3), which is why anything with many features is trained by gradient descent instead.

Function signature

def fit_normal_equation(X, y):

Example

Input

X = np.array([[1.0], [2.0], [3.0]])
y = np.array([3.0, 5.0, 7.0])
fit_normal_equation(X, y)

Output

(array([2.]), 0.9999999999999885)

Reasoning

The points lie exactly on y = 2x + 1. With a column of ones added, the normal equation solves for theta = [1, 2], so b = 1 and w = [2], with zero error.

Constraints and notes

• Add a leading column of ones for the intercept.
• Use np.linalg.pinv, so collinear features do not fail.
• Return w with shape (d,) and b as a float.