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

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?

• Can you post sample size and the error estimates for the different outer folds? – cbeleites May 29 '17 at 7:29
• +1 to cbeleites comment as well as post. If your data is large enough potentially using a hold-out set to report performance statistics might be the easiest option. – usεr11852 May 29 '17 at 9:45
• I updated with a full experiment – user0 May 31 '17 at 17:53
• @usεr11852 how large does it have to be for that to work? – user0 May 31 '17 at 18:41
• I would not suggest any particular hard limit but the hold-out set would have to be at least in the tens of thousands. – usεr11852 May 31 '17 at 23:24

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.)