Bayesian optimization on a polynomial regression as the surrogate model My understanding of Bayesian optimization is that it is generally used in conjunction with Gaussian process (GP) as the surrogate model. Because GP inherently produces an uncertainty of estimate, the acquisition function of expected improvement (EI) can take that into account the uncertainity and strike a balance between exploration and exploitation.
If we were to use a (2nd order) polynomial regression model as the surrogate, as what is normally used in response surface methodology (RSM), is it possible to still use expected improvement as the acquisition function? If so, what is the method to obtain an uncertainty of estimate at any point from a 2nd order polynomial model in order to calculate the EI?

g g already answered this in their comment, but since this led to a further question in the comments, I will try to clarify this further.
A Gaussian process is equivalent to Bayesian linear regression with the inputs mapped into some high-dimensional space such that the dot product of two mapped inputs corresponds to the kernel function applied to the original inputs. See Rasmussen & Williams for a detailed exploration of how the Bayesian linear regression 'weight-space view' and the 'function-space view' are equivalent: https://gaussianprocess.org/gpml/chapters/RW.pdf
What this means for your purposes here is what you are trying to do is equivalent to fitting a Gaussian process with a polynomial kernel. You could therefore fit a Gaussian process with a degree-2 polynomial kernel to generate your uncertainty estimates, and that would be equivalent. Alternatively, you could use Bayesian linear regression. To make this more concrete, let's say you have two input variables, $x_1$ and $x_2$, so that each input datapoint $x_i$ has the associated vector $x_i = [x_{i_1},x_{i_2}]$. If using a 2nd degree polynomial, you can now map each input datapoint to the following $z_i$ (I would standardize your data first...):
$z_i = [x^2_{i_1}, x^2_{i_2},  x_{i_1}x_{i_2}, x_{i_1}, x_{i_2}]$
and fit the mapped input data to a linear model using Bayesian linear regression (e.g. in scikit-learn -- the scikit-learn implementation of Bayesian linear regression will calculate the uncertainty on each prediction for you), and use this to generate your uncertainty values. Polynomial regression is just linear regression with polynomial basis functions.
The two approaches (GP with a polynomial kernel, or Bayesian linear regression with polynomial basis functions) are equivalent. Which one to use depends on whether you have a small number of datapoints in a high dimensional input space (in this case a GP with a polynomial kernel will be faster), or whether you have a large number of datapoints in a low dimensional input space (in this case, performing a matrix decomposition on the kernel matrix will be expensive, so use Bayesian linear regression).
