Can a linear regression model without polynomial features overfit?
It sure can!
Throw in a bunch of predictors that have minimal or no predictive ability, and you’ll get parameter estimates that make those work. However, when you try it out of sample, your predictions will be awful.
set.seed(2020)

# Define sample size
#
N &lt;- 1000

# Define number of parameters
#
p &lt;- 750

# Simulate data
#
X &lt;- matrix(rnorm(N*p), N, p)

# Define the parameter vector to be 1, 0, 0, ..., 0, 0
#
B &lt;- rep(0, p)#c(1, rep(0, p-1))

# Simulate the error term
# 
epsilon &lt;- rnorm(N, 0, 10)

# Define the response variable as XB + epsilon 
#
y &lt;- X %*% B + epsilon

# Fit to 80% of the data
#
L &lt;- lm(y[1:800]~., data=data.frame(X[1:800,]))

# Predict on the remaining 20%
#
preds &lt;- predict.lm(L, data.frame(X[801:1000, ]))

# Show the tiny in-sample MSE and the gigantic out-of-sample MSE
#
sum((predict(L) - y[1:800])^2)/800
sum((preds - y[801:1000,])^2)/200

I get an in-sample MSE of $7.410227$ and an out-of-sample MSE of $1912.764$.
It is possible to simulate this hundreds of times to show that this wasn't just a fluke.
set.seed(2020)

# Define sample size
#
N &lt;- 1000

# Define number of parameters
#
p &lt;- 750

# Define number of simulations to do
#
R &lt;- 250

# Simulate data
#
X &lt;- matrix(rnorm(N*p), N, p)

# Define the parameter vector to be 1, 0, 0, ..., 0, 0
#
B &lt;- c(1, rep(0, p-1))

in_sample &lt;- out_of_sample &lt;- rep(NA, R)

for (i in 1:R){
  
  if (i %% 50 == 0){print(paste(i/R*100, &quot;% done&quot;))}

  # Simulate the error term
  # 
  epsilon &lt;- rnorm(N, 0, 10)
  
  # Define the response variable as XB + epsilon
  #
  y &lt;- X %*% B + epsilon
  
  # Fit to 80% of the data
  #
  L &lt;- lm(y[1:800]~., data=data.frame(X[1:800,]))
  
  # Predict on the remaining 20%
  #
  preds &lt;- predict.lm(L, data.frame(X[801:1000, ]))

  # Calculate the tiny in-sample MSE and the gigantic out-of-sample MSE
  #
  in_sample[i] &lt;- sum((predict(L) - y[1:800])^2)/800
  out_of_sample[i] &lt;- sum((preds - y[801:1000,])^2)/200
}

# Summarize results
#
boxplot(in_sample, out_of_sample, names=c(&quot;in-sample&quot;, &quot;out-of-sample&quot;), main=&quot;MSE&quot;)
summary(in_sample)
summary(out_of_sample)
summary(out_of_sample/in_sample)

The model has overfit badly every time.
In-sample MSE summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.039   5.184   6.069   6.081   7.029   9.800 

Out-of-sample MSE summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  947.8  1291.6  1511.6  1567.0  1790.0  3161.6 

Paired Ratio Summary (always (!) much larget than 1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  109.8   207.9   260.2   270.3   319.6   566.9