I have a regression problem with 5-6k variables. I divide my data into 3 non-overlapping sets: training, validation, and testing. I train using only the training set, and generate a lot of different linear regression models by choosing a different set of 200 variables for each model (I try about 100k such subsets). I score a model as $\min(R^2_{\text{training data}}, R^2_{\text{validation data}})$. Using this criterion, I end up choosing a model. It turns out that the model chosen has very similar $R^2$ on the training and the validation data. However, when I try this model on the testing data, it has much lower $R^2$. So it seems I am somehow overfitting on both the training and the validation data. Any ideas on how can I get a more robust model?
I tried increasing the training data size, but that didn't help. I am thinking of perhaps shrinking the size of each subset.
I have tried using regularization. However, the models I obtain using the lasso or the elastic net have much lower $R^2$ on the training set as well as the validation set, as compared to the model I obtain by doing the subset selection approach. Therefore, I don't consider these models, since I assume that if Model A performs better than Model B on both the training set as well as the validation set, Model A is clearly better than Model B. I would be very curious if you disagree with this.
On a related note, do you think $R^2$ is a bad criteria for choosing my models?