Implementing cross-validation to tune the hyperparameters of an unsupervised model I have extensively researched the application of cross validation for unsupervised learning  (as it is a requirement by my project manager) but it seems that there is no clear consensus as to how to put it in action.
I am currently working on a Clustering project and I have selected 2 algorithms I want to tune :

*

*K-means Clustering

*HDBSCAN.

The business need is the followings : separating different customer segments and identifying each segment characteristics.
My performance metric is DBCV (Density Based Clustering Validation), which I implemented through the HDBSCAN validity_index.
The first problem is that GridSearchCV, RandomSearchCV and BayesSearchCV are not compatible with scoring metrics that do not have ground truth (y_true), so I have to implemement cross validation and grid search manually.
Here is what I have for Kfold cross validation of my Kmeans algorithm (k as the only hyperparameter):
from sklearn.model_selection import KFold
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from hdbscan.validity import validity_index

scores = []
kf = KFold(n_splits=3)
for train_index, test_index in kf.split(X_train_tr):
    X_train = X_train_tr[train_index]
    X_test = X_train_tr[test_index]
    for k in range (2,16):
        print(k)
        km = KMeans(n_clusters=k, random_state=10, n_init = 10)
        km.fit(X_train)
        cluster_labels = km.predict(X_test)
        ssd = km.inertia_
        sil = silhouette_score(X_test, cluster_labels)
        val = validity_index(X_test, cluster_labels)
    
        scores.append({'K': k, 'SSD': ssd, 'Silhouette': sil, 'Validity': val})

Manual Gridsearch with HDBSCAN :
import hdbscan
from sklearn.model_selection import ParameterGrid

hdb = hdbscan.HDBSCAN(gen_min_span_tree=True)

# specify parameters and distributions to sample from
param_grid = list(ParameterGrid(
    {
    'min_samples': [10,100,300,500,1000],
    'min_cluster_size': [100,1000,3000,5000,10000,20000],
    'cluster_selection_method': ['eom','leaf'],
    'metric': ['euclidean','manhattan']
}))

scores = []
i = 0
kf = KFold(n_splits=3)
for train_index, test_index in kf.split(X_train_tr):
    X_train = X_train_tr[train_index]
    X_test = X_train_tr[test_index]
    #Performing manual gridsearch
    for params in param_grid:
        print(i)
        i+=1
        hdb = hdbscan.HDBSCAN(prediction_data=True,
                              gen_min_span_tree=True,
                              min_samples=params['min_samples'],
                              min_cluster_size=params['min_cluster_size'],
                              cluster_selection_method=params['cluster_selection_method'],
                              metric=params['metric']).fit(X_train)
        
        cluster_labels = hdbscan.approximate_predict(hdb,X_test)[0]
        #sil = silhouette_score(X_test, cluster_labels)
        val = validity_index(X_test, cluster_labels)
        val_train = hdb.relative_validity_
    
        scores.append({'min_samples': params['min_samples'],
                       'min_cluster_size': params['min_cluster_size'],
                       'cluster_selection_method': params['cluster_selection_method'],
                       'metric': params['metric'],
                      # 'Silhouette': sil,
                       'Validity Test': val,
                      'Validity Train': val_train})

scores_df = pd.DataFrame(scores).groupby(
    ["min_samples",'min_cluster_size','cluster_selection_method','metric']).mean().reset_index()

This is working, in the sense that it gives sensible results, but something in the back of my mind tells me that this is not right. The point of the model is to figure out the best clusters based on the whole dataset, not to predict the best clusters on new data points.
That being said, one of the business requirement was to have stable clusters, which explains why I think cross validation could help selecting the best model and hyperparameters.
So here are the 3 choices I'm pondering :

*

*Should I cross-validate based on the predicted clusters on my test set? (as implemented with the code above)

*Should I cross-validate only based on the performance of the clustering on the train set? (using the different splits to make sure my model is consistent)

*Should I use another way of cross validating my results?


So I've reached a tentative answer to my problem, which I don't claim is better than @PatrickBormann 's edited answer but allowed me to reach sensible conclusions, which I'm sharing now for other people.
The best way in my opinion to perform the hyperparameter tuning of an unsupervised machine learning algorithm (when stability is a concern) is the method outlined here How to perform Validation on Unsupervised learning? which is described in the following figure (credits to @Vincentdebakker) :

I have used the implementation of Predictive Strength (PS) described in the 100 page ML book Github (https://github.com/aburkov/theMLbook/blob/master/prediction_strength.py).
The problem for my dataset is that Prediction Strength as it is implemented in Python takes a very long time to compute. My (pretty performant) computer was not able to predict a single prediction strength in 2 hours, so for my dataset of ~100k samples and 33 features it would not be feasible to use that metric during cross validation.
I have resorted to using the method outlined here : https://amueller.github.io/aml/04-model-evaluation/17-cluster-evaluation.html which gave me sensible results and allowed me to optimize the DBSCAN relatively quickly.
Here is my implementation with a manual gridsearch for DBSCAN :
from sklearn.model_selection import ParameterGrid
from sklearn.model_selection import KFold
from sklearn.cluster import DBSCAN
from sklearn.metrics import silhouette_score, adjusted_rand_score
import gc

#Specify parameters to sample from
param_grid = list(ParameterGrid(
    {
    'min_samples': [200,300,400,500,700,900],
    'eps' : [0.25,0.35,0.45,0.55,0.65],
    'metric': ['euclidean','manhattan']
}))

scores = []
kf = KFold(n_splits=3)
SIL = []
i = 0
#Performing manual gridsearch
for params in param_grid:
    
    labels = []
    indices = []
    scores_stab = []
    SIL = []
    
    for train_index, test_index in kf.split(X_train_tr):
        X_train = X_train_tr[train_index]
        X_test = X_train_tr[test_index]

        dbs = DBSCAN(n_jobs=-1,
                      min_samples=params['min_samples'],
                      eps=params['eps'],
                      metric=params['metric']).fit(X_train)
        clusters = dbs.labels_
        sil = silhouette_score(X_train, clusters)
        SIL.append(sil)
        
        #Calculating stability on the train set
        indices.append(train_index)
        relabel = -np.ones(X_train_tr.shape[0], dtype=np.int32)
        relabel[train_index] = clusters
        labels.append(relabel)
        
        min_samples_pct = params['min_samples'] / len(X_train)
        
        gc.enable()
        del X_train, X_test
        gc.collect()
    
    for l, i in zip(labels, indices):
        for m, j in zip(labels, indices):
        # we also compute the diagonal
            in_both = np.intersect1d(i, j)
            scores_stab.append(adjusted_rand_score(l[in_both], m[in_both]))

    print("min_samples: {}, eps: {}, metric: {} ==>  Silhouette: {}, Stability {}".format(
    params['min_samples'], params['eps'], params['metric'], np.mean(SIL), np.mean(scores_stab)))        
    scores.append({'min_samples': params['min_samples'], 'min_samples_pct': min_samples_pct,'eps': params['eps'],
                   'metric': params['metric'],'Silhouette': np.mean(SIL), 'Stability': np.mean(scores_stab)})
    
scores_dbs = pd.DataFrame(scores)

scores_dbs

It should be underlined that the stability calculation as defined in this line of code is not as rigorous as the cross validation based on 10 iterations defined in the link above, but I think this is still a good approximation that can be merged within the cross validation. I may be wrong though because the clusters in my dataset are quite stable so maybe this implementation would fail with less cleanly separated clusters.
Important note the min_samples hyperparameter defined here is relative to the number of samples it is calculated on (here 2/3 of my initial X_train (3 folds) ,and I kept 10% of my initial dataset as a test set). That's why it's important to calculate the min_samples to dataset length ratio that will then be multiplied by the length of the whole dataset when applying the clustering algorithm to my initial dataset.
