# Ordinary least squares does not optimize error

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()

• There are many, many things that can go wrong, which is why you really need to display your results after reducing your problem to a minimal reproducible example. I suspect (without any proof yet) that you are running into floating point problems due to the severe near-collinearity of all those powers of x. If that's the case, the solution is to use an orthogonal polynomial basis. (A far better solution is never to do polynomial regression, but presumably this is not intended as a practical example.) – whuber Dec 8 '18 at 21:25
• @whuber: this code is a minimai reproducible example. It's hard to list the results in any other form. Theoretically though, do I correctly understand that increasing variables should produce at least equally small MSE on the training set, by definition? – blue_note Dec 9 '18 at 11:18
• You haven't given us an example yet: in addition to the code, exactly what does it output? – whuber Dec 9 '18 at 18:38