Validity of regression diagnostics for deterministic computer experiments The question
Is it possible to justify the use of regression diagnostics such as F-tests or the AIC for the analysis of deterministic computer experiments?
Background
One way to analyse computer experiments is by building surrogate models. The unknown function $f$ is approximated by a nice and explicitly known function $\hat f$, determined from some training data $(x_i)$ and the according responses $(f(x_i))$.
To make the question more concrete it suffices to focus on multivariate polynomials as surrogate functions. They are fit using standard regression tools (e.g. "lm" in R or "sklearn.linear_model.LinearRegression" in Python). Furthermore, since design matrices can get large very quick, methods relying on variable selection such as stepwise regression are often used.
Notice that in most cases neither of the two main regression assumptions are fulfilled: More often than not the unknown function $f$ will be ill-specified, i.e. not be a polynomial of any degree, and there will be no noisy observations, i.e. independent errors with a normal distribution.
My issue
The minimisation of residual sum of squares makes perfect sense for an approximation task. Hence, I do understand the use of regression tools in general and diagnostics such as $R^2$ even though the regression assumptions are violated.
But it seems to me that diagnostics such as F-tests and the AIC do rely crucially on the likelihood or other statistical assumptions. While tools calculate them readily (but quite formally it seems), I do not see any reason why stepwise regression or variable selection based on such diagnostics should work in the deterministic and ill-specified case.
That said, I have seen this approach used in practice and it is also mentioned in textbooks (e.g. Chapter 7.2.2 of [1]). So maybe there is something to it?
In other words: Can stepwise regression or variable selection using AIC (or similar statistics based on the likelihood or distributional assumptions of the regression model) be justified also in the case of deterministic computer experiments?
[1] "The Design and Analysis of Computer Experiments" by Thomas J. Santner, Brian J. Williams, William I. Notz

This question reminds me a bit of the question Could a mismatch between loss functions used for fitting vs. tuning parameter selection be justified? , where I answered that it can be useful to use a different loss function than the loss function which needs to be optimized. The reason is that statistical fluctuations might be best minimized by using a loss function that relates to the statistical distribution of the errors. For example: when we are estimating the location parameter of a Laplace distribution, then finding an estimator that minimizes the sum of absolute residuals will be better in optimizing the expected square error than an estimator that minimizes the sum of squared residuals.
In your example there is also in some way statistical variation. You have a deterministic function (the population), but you only take a sample $x_i$ where you evaluate $f(x_i)$, so you have a sub-selection of the population*, and that resembles the randomness of sampling.

I do not see any reason why stepwise regression or variable selection based on such diagnostics should work in the deterministic and ill-specified case.

The error distribution of this 'randomness' is not very clear, and to optimize AIC (based on an arbitrary likelihood) or F statistic will only make sense as a heuristic diagnostic. Yet they do work and that is because they introduce some way to penalize the number of terms. You can fit a high order polynomial and reach an $R^2 = 1$ but that is possibly introducing a larger prediction error.

Instead, to reduce this error, it might be better to use, for instance, leave one out cross validation, which relates to the mean squared error (Proof of LOOCV formula), and is asymptotically similar to optimizing a likelihood (Equivalence of AIC and LOOCV under mismatched loss functions)

*(If you would know the entire function $f$ then you could just optimize the entire function straight away, and if you have no limits to the number of terms then you get something like a perfectly fitting Taylor series or Fourier series).
