I'm trying to use polynomial regression to fit the curve
$$X = [0, 1]$$ $$Y = \sin(2 \pi X) + \epsilon$$
where $\epsilon$ is normally distributed with the same $\sigma$ for all $X$
For every value of $x$ I'm creating the vector $[x^0, x^1, \ldots x^n]$ (that is, consisting of the value $x$ raised to the powers $[0, n]$) and trying to fit it using linear regression with ordinary least squares. However, as the size of $n$ increases (empirically, beyond $10$), the MSE gets significantly worse (the result of also obvious when plotting)
How does that happen? Isn't ordinary least squares guaranteed to find the optimal solution? Since for $n=2$ the solution is much better MSE than for $n=10$, the estimator $[b_0, b_1, b_2, 0, 0, 0, \ldots]$ would produce better MSE with $n=10$ than the one that ordinary least squares finds. What I'm I misunderstanding?
Below is sample python code (if anyone is interested in running it)
import numpy as np from matplotlib import pyplot as plt def inputs(X, n): return np.column_stack([X**i for i in range(n+1)]) def targets(X): T0 = np.sin(2* np.pi * X) + 5 T = T0 + 0.1 * np.random.randn(*X.shape) return T, T0 class Regressor: def fit(self, X, T): T = T.reshape(-1, 1) self.B = np.linalg.inv(X.T @ X) @ X.T @ T def predict(self, X): Y = X @ self.B return Y.ravel() num_samples = 100 num_dimensions = 20 regressor = Regressor() X = np.linspace(0, 1, num_samples) T, T0 = targets(X) X_ = inputs(X, num_dimensions) regressor.fit(X_, T) Y = regressor.predict(X_) MSE = ((Y - T)**2).sum() plt.plot(X, T0, 'C0') plt.plot(X, T, 'oC0') plt.plot(X, Y, 'C1') plt.show()