Why is r squared lowered when adding polynomial features? I am trying to find a best fit line f(x) = ? for a random set of x,y coordinates.

Linear Regression with polynomial features works well for around 10 different polynomials but beyond 10 the r squared actually starts to drop!

If the new features are not useful to the Linear Regression I would assume that they would be given a coefficient of 0 and therefore adding features should not hurt the overall r squared.

I reproduced this problem when housing price predictions when creating a large amount of interaction features.

I have my python code below:

Create Random Data

import numpy as np
import matplotlib.pyplot as plt


def pol(x):
    return x * np.cos(x)

x = np.linspace(0, 12, 100)
rng = np.random.RandomState(1234)
rng.shuffle(x)
x = np.sort(x[:25])
y = pol(x) + np.random.randn(25)*2


plt.scatter(x, y, color='green', s=50, marker='.')

plt.show()


Regress and Check Each R Squared

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

for p in range(1,30):
    plot_range = [i/10 for i in range(0,120)]
    poly = PolynomialFeatures(p)
    X_fin = poly.fit_transform([[samp] for samp in x])
    X_fin_plot = poly.fit_transform([[samp] for samp in plot_range])
    reg = LinearRegression().fit(X_fin, y)

    from sklearn.metrics import mean_squared_error, r2_score
    print(p,r2_score(y, reg.predict(X_fin)))


Display Last Regression Line

plt.scatter(x, y, color='green', s=50, marker='.')
plt.plot(plot_range,reg.predict(X_fin_plot))
plt.show()


I also have two plots to compare. The first is with 10 polynomial features and the second is with 40. Notice how the second misses the majority of the first points.
My original answer was not correct, so here is a corrected answer:

When you use PolynomialFeatures(), you don't get the intended polynomials. Instead you get polynomials plus an interaction term:

from sklearn.preprocessing import PolynomialFeatures import numpy as 
np    z = np.array([[0, 1],
                    [2, 3],
                    [4, 5]]) 
poly = PolynomialFeatures(2)
print(poly.fit_transform(z))


Output is:

[[ 1.  0.  1.  0.  0.  1.]
 [ 1.  2.  3.  4.  6.  9.]
 [ 1.  4.  5. 16. 20. 25.]]


A raw polynomial should look like:

new_z = np.hstack((z**(i+1) for i in range(2)))
print(new_z)


Output is:

[[ 0  1  0  1]
 [ 2  3  4  9]
 [ 4  5 16 25]]


Here is a quick R implementation of your problem with raw polynomials:

x = c(0.12121212, 1.09090909, 3.27272727, 3.51515152, 4, 4.24242424,
  4.72727273, 4.84848485, 5.09090909, 6.18181818, 6.78787879, 7.15151515,
  7.39393939, 7.63636364, 8.24242424, 8.60606061, 9.09090909, 9.81818182,
  9.93939394, 10.3030303, 10.54545455, 10.66666667, 11.39393939, 11.63636364,
  11.87878788)

y = c(-2.87011136,1.77132943,-1.23698978,-3.09768628,-2.11919042,-4.11234626,
  -1.1684339, 1.34601699, -2.37623758,4.20290438, 6.16349341, 3.60661197,
  2.58898819, 3.80785471, -2.96359566, -5.672873, -9.71694313, -7.62778351,
  -8.95730409, -8.04664475, -5.18464423, -6.54562138, 3.45527603, 6.11936457,
  9.30106747)

regdata = data.frame(x,y)
colnames(regdata) &lt;- c("x","y")

r2list = list()
r2adjlist = list()
plist = list()

for (p in seq(1:29)){
  reg = lm(y~poly(x,p, raw=T), data=regdata)
  print(paste0("Poly: ", p))
  print(paste0("  R2      ", summary(reg)$r.squared))
  print(paste0("  R2_adj. ", summary(reg)$adj.r.squared))
  r2list[[p]] &lt;-  summary(reg)$r.squared
  r2adjlist[[p]] &lt;- summary(reg)$adj.r.squared
  plist[[p]] &lt;- p
}

plot(plist, r2list,xlab="Polynomial", ylab="R2")
lines(plist, r2list)


The R2 contingent on the degree of the polynomial is shown below:


So your initial intuition was (of course) correct, but your treatment of data was not correct.