If choice of learning algorithm is an unstable hyperparameter in (nested) CV, is the estimate of generalization error still valid?

Question

Suppose I do a grid search over different algorithms (eg, LASSO, RF, MLP), fitting the choice of algorithm as a hyperparameter in nested CV using Approach 2 from this answer by @cbeleites (Approach 2: "average performance of all winning models regardless of their individual parameter sets"---unlike Approach 1, this does not require stable optimization results).

The model (I refer to model as parameters+hyperparameters (one of which is the learning algorithm)) is unstable so that there are different winning learning algorithms per outer fold. (Given the flexibility of learning algorithms and also with a randomized search over each model's hyperparameters, and also because it seems that some models are both very good even with different parameters (local minima?), I expect that this kind of instability happens even more frequently when searching over both the model-specific and model choice hyperparameters.)

As an example, I run 5 outer folds

outer fold 1 LASSO wins with C=0.03 
outer fold 2 LASSO wins with C=0.02
outer fold 3 RF wins with max depth=4
outer fold 4 LASSO wins with C=0.03
outer fold 5 LASSO wins with C=0.03

Now, we have 5 estimates of the generalization error, one of which is from the RF. Average these 5 estimates and call it the estimate of the generalization error. Repeat an identical modeling process on the full dataset (again, treating choice of algorithm as a hyperparameter). Even when the hyperparameters were not stable across folds (where hyperparameters include intrinsic learning algorithm hyperparameters and also the choice of learning algorithm itself), does the generalization error estimate hold? As described, in Approach 2 the error estimate does hold even when the hyperparameters are different over folds. But, the question concerned a a single type of learning alogrithm---is this the same when the learning algorithm is a also a hyperparameter? I think it should be, given @cbeleites description and I think Cawley and Talbot, but I wanted to confirm. If so, this would open the door to some very interesting optimization methods for me, so thanks for your help!

Update:

I didn't have a particular dataset when I asked this question, but I designed an experiment in response to the comments. I am using the breast cancer dataset from sklearn. It has 569 observations and 30 features with a binary response. I am performing 2 hide/see splits (I hide 284 obs and "see" 285" obs). I perform nested CV on the 284 "seen" cases. For nested CV, I make 10 outer splits so that they are each roughly 255 train and 29 test. I make 5 inner splits to optimize hyperparameters, which in this case include the models themselves and their respective hyperparameters:

ess = [
    lr(penalty='l1'),
    lr(penalty='l2'),
    rf(),
    mlp(),
    svc(kernel='rbf',probability=True)
  ]

And their respective parameters:

paramss = [
        {
        'C': uniform(1e-10,10)
        },
        {
        'C': uniform(1e-10,10)
        },
        {
          "max_features": randint(2, 29),
          "min_samples_split": randint(2, 1000),
          "min_samples_leaf": randint(2, 1000),
          "bootstrap": [True, False],
          "criterion": ["gini", "entropy"],
          "max_depth": randint(2, 100)
        },
        {
        'alpha':uniform(1e-10,10)
        },
        {
        'C': uniform(1e-10,10),
        'gamma':np.logspace(-9, 3, 13)
        }
        ]

I am doing a randomized CV search for each inner fold, but the random seed is set so that each fold sees the same set of hyperparameters. For each, I do 10 iterations of random search (so 10 combinations of parameters).

Results:

pr :
     Est 0.983813310117+/-0.0322076978758
     True 0.99458419491+/-0.00188667067285
ll :
     Est 0.155164194443+/-0.16402795992
     True 0.121746396698+/-0.0457845523029
roc :
     Est 0.983693206449+/-0.0267251644912
     True 0.991010575306+/-0.00386771816282
bri :
     Est 0.0353970233858+/-0.0265623516895
     True 0.0361436496984+/-0.0160694593715

The following shows the selected parameters that was evaluated on each outer fold collected over folds, so if a particular learning algorithm never won, it has a blank list for its parameters. You seen that Ridge won 12, LASSO won 7, and SVM won 1 (we have a total of 20 outer folds because we repeat the experiment twice with 10 outer folds in each trial).

'RF': {'bootstrap': [], 'min_samples_leaf': [], 'max_features': [], 'criterion': [], 'min_samples_split': [], 'max_depth': []}, 

'SVM': {'C': [7.9402135046538351], 'gamma': [0.0001]}, 

'Ridge': {'C': [2.3754122004491229, 2.3754122004491229, 2.3754122004491229, 6.4161334476906919, 2.3754122004491229, 2.3754122004491229, 2.3754122004491229, 4.5344924742731223, 2.3754122004491229, 2.3754122004491229, 7.2201822952694714, 2.3754122004491229]}, 

'MLP': {'alpha': []}, 

'LASSO': {'C': [2.3754122004491229, 2.3754122004491229, 2.3754122004491229, 6.0904246277127791, 9.7260111391489339, 9.657491980529997, 9.7260111391489339]}

Further note:

np.std(param_res['LASSO']['C'])
3.3922064823749056
np.std(param_res['Ridge']['C'])
1.6905531536433458

In general, the estimated error is pretty close to the real error, although the variance of the log loss is huge--which seems to always be the case. I am using log loss to optimize these parameters as well, so it's odd.

Repeated with a synthetic nonlinear dataset:

from sklearn.datasets import make_gaussian_quantiles
X, y = make_gaussian_quantiles(n_features=2, n_classes=2, n_samples=600)


pr :
     Est 0.996811863475+/-0.00551945904917
     True 0.989877985512+/-0.009249290943
ll :
     Est 0.0755535031225+/-0.038931334513
     True 0.0969572501675+/-0.051841178538
roc :
     Est 0.997278446635+/-0.00464250919492
     True 0.992784285778+/-0.00636942675159
bri :
     Est 0.0224036276652+/-0.0151780407018
     True 0.0179532075878+/-0.0068763012141


    defaultdict(<type 'dict'>, {
'RF': {'bootstrap': [], 'min_samples_leaf': [], 'max_features': [], 'criterion': [], 'min_samples_split': [], 'max_depth': []}, 
'SVM': {'C': [7.2201822952694714, 7.2201822952694714, 7.7770241058382013, 7.5858400355869096, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.5858400355869096, 7.5858400355869096, 7.5858400355869096, 7.2201822952694714, 7.5858400355869096, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.5858400355869096, 7.5858400355869096, 7.2201822952694714, 7.7770241058382013], 'gamma': [1.0, 1.0, 10.0, 10.0, 1.0, 1.0, 1.0, 1.0, 10.0, 10.0, 10.0, 1.0, 10.0, 1.0, 1.0, 1.0, 10.0, 10.0, 1.0, 10.0]}, 
'Ridge': {'C': []}, 
'MLP': {'alpha': []}, 
'LASSO': {'C': []}})

Still, the estimate of log loss and brier score are quite variable.

If we increase sample size from 600 to 5000, everything's better:

X, y = make_gaussian_quantiles(n_features=2, n_classes=2, n_samples=5000)


pr :
     Est 0.999776032821+/-0.000273469116399
     True 0.999877941448+/-2.88452390002e-05
ll :
     Est 0.0227828679184+/-0.00713550673936
     True 0.0221161827013+/-0.00279061197794
roc :
     Est 0.999775036095+/-0.000274563194329
     True 0.999877689549+/-2.36936475411e-05
bri :
     Est 0.0062857004316+/-0.00280176776558
     True 0.00600250426244+/-0.000664927935621



defaultdict(<type 'dict'>, {
'RF': {'bootstrap': [], 'min_samples_leaf': [], 'max_features': [], 'criterion': [], 'min_samples_split': [], 'max_depth': []}, 
'SVM': {'C': [7.2201822952694714, 7.7770241058382013, 7.2201822952694714, 7.5858400355869096, 7.2201822952694714, 7.7770241058382013, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714, 7.2201822952694714], 'gamma': [1.0, 10.0, 1.0, 10.0, 1.0, 10.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}, 
'Ridge': {'C': []}, 
'MLP': {'alpha': []}, 
'LASSO': {'C': []}})

Even in this case, however

(Pdb) np.std(param_res['SVM']['C'])
0.17950329298698336
(Pdb) np.std(param_res['SVM']['gamma'])
3.2136427928442828

The C and gamma tend to covary, however. So perhaps when there are two hyperparameters that could each have similar effects on model fitting, it's impossible to find one optimization that gives best performance, exactly like when there is multicollinearity, it's impossible to find a single value for model parameters, since you can always subtract from one whatever you add to the other?

Answer 1

Yes, I think the estimate holds for "you can achieve this performance* by training a model that auto-tunes hyperparameters this way".

Still, there's always the question: if optimiziation is not stable, did it work at all?
So the important point is that you need to think what this instability means.

You may have a situation where there truly is no single best set of hyperparameters. While you could reasonably expect that for a single hyperparameter like the penalty of a LASSO a (global) optimum exists. But when comparing different algorithms you may have very similar performance achieved by different algorithms (think e.g. LDA, logistic regression and linear SVM).

• Best case: easy problem, different hyperparameter sets achieve perfect classification

• a not so good variant of this is: splitting into the CV sets does not yield independent sets of cases (e.g. because of hierarchical/clustered data structure) -> there may be heavy overfitting as a consequence -> if the inner CV sees only "perfect" classification, it doesn't know how to choose.
  => check inner performance estimates and use proper scoring rules for optimization

• not perfect, but e.g. the 3 linear models mentioned above all show up in the winner parameters -> linear model achieves xxx performance, and it doesn't matter which one you choose.

* IMHO you need to put all interpretations into relation of what you know about "this performance": measurement of generalization error has variance uncertainty (and bias) like any other measurement. Bias is going to be slightly pessimistic if the splitting is done properly in your case, but variance uncertainty crucially depends on the number of cases available.

So the questions are:

• Is the uncertainty on the final (outer) performance estimate small enough to distinguish whether the model answers your application problem well or not so well or not at all?
• Is the uncertainty on the inner performance estimates small enough to make selection of an apparently best model a sensible approach?

One question I'd always ask as well: compared to a simple "baseline" model where you fix all hyperparameters in advance by your knowledge about the problem, did the optimization yield any improvement (keeping in mind the uncertainty)?
(In my field I encounter problems where - due to the very small sample size that is available - I can say after checking such a baseline model performance: I don't need to start an optimization because I won't be able to show that the "optimal" model is better than the one I already have unless I get far more cases.)