I've read many threads on this website try to understand why we need to break the data on hand into 3 parts, the training, validation and test data set. I am still thinking it is enough just to break the data set on into into 2 pieces, i.e., the training and validation, then I can fit different models on the training set and apply them on the validation set, then eventually I will choose the model that has best performance on the validation set and I should expect this model will also beat other models I've fitted on a completely independent real test data set provided that the training, validation and the independent test data set all representative of the population. Therefore, I claim that the validation score is enough for model selection (or hyper parameter tuning).
First, I've experienced empirically that my above claim is wrong here. However, I've come up with the following theoretical argument that I cannot refute myself. I hope one could point out where I made mistakes in the following argument.
A Machine Learning Model has 2 types of parameters: (1) Hyper Parameter, e.g., the regularization parameter $\lambda$ in LASSO, denoted by $\boldsymbol{\lambda}$, which is also an input stuff that you have to specify (together with the data) before you run the model; (2) Weight Parameter, e.g., the linear coefficients in LASSO, denoted by $\boldsymbol{w}$, which is automatically generated by the model after you specified $\boldsymbol{\lambda}$ and the data. Note that in tree-based models, $\boldsymbol{w}$ does not have explicit expressions, in that case, $\boldsymbol{w}$ represents the specific branches of the tree.
For a given class of machine learning algorithm, e.g., LASSO, for different value of the hyper parameter, it is essentially a different model.
Let the data set you have on hand that have both the value of the predictors and the value of the response variable be $D$ and $f(\cdot |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D)$ be the model you use to submit your prediction over the independent test data set $T$ which you only observe the value of the predictors but not the value of the response variables. Note that in $f(\cdot |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D)$, the value of $\boldsymbol{\lambda}$ is specified by you, the value of $\boldsymbol{w_\lambda}(D)$ is determined both by $\boldsymbol{\lambda}$ and the training data set $D$. With this function, if you plug in a value of the predictors, it will output a prediction about the corresponding value of the response variable. Let $L(y,f(x|\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D))$ be the loss function of the true response value versus the prediction.
According to Elements of Statistical Learning (ESL), we are looking at two types of error:
Conditional Generalization Error: \begin{equation} Err_D=\mathbb{E}[L(Y,f(X |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D))|D] \end{equation} where $(Y,X)$ are drawn from the joint distribution $F_{X,Y}$ and $D$ is fixed. ESL says that this conditional generalization error is our goal.
Comment: this makes sense, since we only have data set $D$ to develope our model, and we want to know exactly, based on $D$ and the model we developed, what is its performance over unseen data that come from $F_{X,Y}$. Theoretically, we would like to know the distribution of $L(Y,f(X |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D))|D$. And theoretically, this can be done by draw i.i.d. samples $\{(y_i,x_i)\}_{i=1}^m$ from $F_{X,Y}$ and plug $x_i$ into $f(\cdot |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D)$ and compute the $L(y_i,f(x_i |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D))|D$ and finally, if you do this for enough number of times, you can recover the distribution $L(Y,f(X |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D))|D$. Practically, you just can't do that.
But practically, you can do the following: leave-one-out cross-validation. Each time you take one row $(y_j,x_i)$ from $D$ and fit the model with the remaining of $D$, denoted the remaining of $D$ by $D_{-j}$ indicating that we've removed one row of $D$. Then you can get $f(x_j |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^{D_{-j}})$. I don't know how to show it theoretically, but intuitively, we should have \begin{equation} f(\cdot |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^{D_{-j}})\approx f(\cdot |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^D)~,\forall~j, \end{equation} as we only removed one row of $D$ and if $D$ is large (i.e., it has many rows), then the above approximation should be really good. Let's put the computational concern aside at this moment. Each value of $f(x_j |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^{D_{-j}})$ should be an unbiased estimator of $Err_D$. So I think that if you $D$ is large and you have enough computational power, then after you do this leave-one-out method and take the average, you should be able to get a very close number to the true $Err_D$. Then you should be able to expect that after you submit your prediction over $T$ (assume $T$ is also large, i.e., has many rows), the error rate on $T$ that should be very close to $Err_D$.
My Question is: ESL on page 242 second last paragraph also mentioned this intuition, however, it just says that it turns out this leave-one-out method does not effectively estimate the $Err_D$ and in the end of this chapter, it also says there is no good way to estimate $Err_D$. It does not give proof or theoretical argument to illustrate the idea of why leave-one-out cannot recover $Err_D$ (it does give empirical evidence).
I believe the leave-one-out method will generate a quantity almost equal to $Err_D$, then I am questioning, why we need to break the data set on hand, i.e., $D$ into 3 pieces: training ,validation and test? Specifically, if I want to find the value of hyper parameter $\boldsymbol{\lambda}^*$ that gives me the lowest $Err_D$, I could just do a grid search over the space of $\boldsymbol{\lambda}$ and for each value of $\boldsymbol{\lambda}$, I just run the leave-one-out cross-validation and it will give me the (almost) TRUE value of $Err_D$ for the specific $\boldsymbol{\lambda}$. Then I could just choose the one that generates the lowest $Err_D$.
In general, if I just randomly split my $D$ into $D_F$ and $D_T$ and I use $D_F$ to fit my $M$ different models (e.g., $M$ different values of $\boldsymbol{\lambda}$, the hyper parameter), then I apply these models on $D_T$. Then for each model and each $(x_i,y_i)\in D_T$, the quantity $L(y_i,f(x_i |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^{D_F}))$ should be an unbiased estimator of $Err_{D_F}$ for model with specific $\boldsymbol{\lambda}$. Then if I take average of $L(y_i,f(x_i |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^{D_F}))$ for $(x_i,y_i)\in D_T$, the resulting number should be pretty close to $Err_{D_F}(\boldsymbol{\lambda})$. Then if I bring another independent test data set (which should be essentially not very different from $D_T$) and apply $f(\cdot |\boldsymbol{\lambda},\boldsymbol{w_\lambda}^{D_F})$ on that new test data set, I should expect to see that the average loss on it is roughly the same as what I've got from $D_T$, isn't it? So then I think this procedure did model selection and performance assessment at the same time. Because the quantity I got for each model on $D_T$ is an unbiased (and could be also with little variance if $D_T$, $D_F$ are all large) estimator of $Err_{D_F}(\boldsymbol{\lambda})$ and our goal is to find a model that has the least average loss when applying to a new test data set. Hence, I could just choose the model with the least $Err_{D_F}(\boldsymbol{\lambda})$, in this case, I not only pin down the best model among all of my candidate models, but also found a good estimate of the average loss rate in predicting new instances. So why we need to split $D$ into 3 pieces?
The below answer is very insightful. When we randomly divide $D$ into $F$ and $V$, then fit the model on $F$ and apply on $V$, the estimates (i.e., average error on $V$) is the realization of the unbiased estimator of $f(\cdot|\lambda,w^{D_F}_\lambda)$. But as $V$ is not infinite, then this unbiased estimator will have non-zero variance and if you pick the value of the hyper parameters that gives you the least realization of this unbiased estimator, then certainly it is upward bias.
Also there is another thing worth mentioning:
Say we have 2 values of the hyper parameter, $\lambda_1$ and $\lambda_2$. Now, let's say we are in the setting of $D=F+V$. You get 2 models from $F$, which is denoted by $f(\cdot|\lambda_1,w^{F}_{\lambda_1})$ and $f(\cdot|\lambda_2,w^{F}_{\lambda_2})$. Take model $f(\cdot|\lambda_1,w^{F}_{\lambda_1})$ as an example. I want to know the generalization error of this model, which can be compute theoretically by: generating infinitely amount of data based on the distribution of $D$, then apply this model to those data and get the average error rate, which is a constant, not a random variable. Note that $F$ and $V$ are both come from the distribution of $D$, as we randomly split $D$ into $F$ and $V$. Now, of course, we cannot generate infinite data. Note that if we only generate one data and apply the model to that data, then the error (rate) is a random variable, but it is still an unbiased estimator (which is a random variable) of the true error rate. Now we apply the model to $V$, which is finite, and we can get an error rate, which is a realization of an unbiased estimator of the true error rate. It has variance, in the sense that if you have another $V'$ with same lengths as $V$, then you will get a different number. Same thing happens for $f(\cdot|\lambda_2,w^{F}_{\lambda_2})$. So for this specific $V$, using $f(\cdot|\lambda_1,w^{F}_{\lambda_1})$, you get a number $e_1$ and using $f(\cdot|\lambda_2,w^{F}_{\lambda_2})$, you get a number $e_2$. You have no idea about how close $e_1$ to $\mathbb{E}[e_1]$ and how close $e_2$ to $\mathbb{E}[e_2]$. So if you pick model 2 because of $e_2<e_1$, this does not guarantee that if you get infinite copies of $V$ (same distribution of $V$ but with different realization), then you have $\mathbb{E}[e_2]<\mathbb{E}[e_1]$.
The second thing is that for using $\lambda_1$, you can do a cross validation say 10 fold, and note that in each fold, the $F_i$ is different and $V_i$ is also different. You will have $e^{(i)}_1$ for $i=1,2,...,10$ and theoretically speaking, each $e_i$ is the realization of a different thing, i.e., model $f(\cdot|\lambda_1,w^{F_i}_{\lambda_1})$ trained on $F_i$ and the error on $V_i$. And you will get a mean of $\{e^{(i)}_1\}_{i=1}^{10}$ and std. It could also be that this std for $\lambda_1$ is greater than std for $\lambda_2$, even if the mean is the opposite.
