Hyperparameter Optimization using chemml.optimization.GeneticAlgorithm

We use a sample dataset from ChemML library which has the SMILES codes and Dragon molecular descriptors for 500 small organic molecules with their densities in \(kg/m^3\).

For more information on Genetic Algorithm, please refer to our paper

from chemml.datasets import load_organic_density
smiles,density,features = load_organic_density()

print(density.shape, features.shape)
density, features = density.values, features.values

from sklearn.preprocessing import StandardScaler
scalerx = StandardScaler()
features = scalerx.fit_transform(features)
density = scalerx.fit_transform(density)

(500, 1) (500, 200)

Defining hyperparameter space

Lets consider kernel ridge regression from scikit-learn for training. The hyperparameters of interest are: alpha, kernel and degree.

The space variable is a tuple of dictionaries for each hyperparameter. The dictionary is specified as:

{'name' : {'type' : <range>}}

An additional mutation key, with its value as: (mean, standard deviation) of a Gaussian distribution, is also required for the ‘uniform’ hyperparameter type.

from sklearn.kernel_ridge import KernelRidge
space = (
        {'alpha'   :   {'uniform' : (0.1, 10), 'mutation': (0,1)}},
        {'kernels' :   {'choice'  : ['rbf', 'sigmoid', 'polynomial', 'linear']}},
        {'degree'  :   {'int'     : (1,5)}} )

Defining objective function

The objective function is defined as a function that receives one ‘individual’ of the genetic algorithm’s population that is an ordered list of the hyperparameters defined in the space variable. Within the objective function, the user does all the required calculations and returns the metric (as a tuple) that is supposed to be optimized. If multiple metrics are returned, all the metrics are optimized according to the fitness defined in the initialization of the Genetic Algorithm class.

from sklearn.metrics import mean_absolute_error
from chemml.utils import regression_metrics
def obj(individual):
    krr = KernelRidge(alpha=individual[0], kernel=individual[1], degree=individual[2])
    krr.fit(features[:400], density[:400])
    pred = krr.predict(features[400:])
    mae = regression_metrics(density[400:],pred)['MAE'].values[0]
    return mae

Optimize the model

from chemml.optimization import GeneticAlgorithm
import warnings
warnings.filterwarnings('ignore')

ga = GeneticAlgorithm(evaluate=obj, space=space, fitness=("min", ),
                    pop_size = 8, crossover_size=6, mutation_size=2, algorithm=3)
fitness_df, final_best_hyperparameters = ga.search(n_generations=5)

ga.search returns:

• a dataframe with the best individuals of each generation along with their fitness values and the time taken to evaluate the model

• a dictionary containing the best individual

fitness_df

	Best_individual	Fitness_values	Time (hours)
0	(5.757142857142857, linear, 2)	0.107350	0.000031
1	(2.846733237062951, linear, 3)	0.102321	0.000026
2	(1.3355690373961644, linear, 3)	0.096700	0.000038
3	(1.3355690373961644, linear, 3)	0.096700	0.000028
4	(1.3355690373961644, linear, 3)	0.096700	0.000039

print(final_best_hyperparameters)

{'alpha': np.float64(1.3355690373961644), 'kernels': 'linear', 'degree': 3}

Resume optimization

The Genetic Algorithm can resume the search for a combination of the best hyperparameters from the last checkpoint. This feature can be useful when the objective function is computationally expensive.

fitness_df_resume, final_best_hyperparameters_resume = ga.search(n_generations=5)

fitness_df_resume

	Best_individual	Fitness_values	Time (hours)
0	(1.3355690373961644, linear, 3)	0.0967	0.000031
1	(1.3355690373961644, linear, 3)	0.0967	0.000035
2	(1.3355690373961644, linear, 3)	0.0967	0.000029
3	(1.3355690373961644, linear, 3)	0.0967	0.000030
4	(1.3355690373961644, linear, 3)	0.0967	0.000032

Hyperparameter optimization for multi-output models

Multi-output models can be broadly divided into two classes for the sake of optimization; scikit-learn like and neural network. This subsection gives a brief pipeline for optimizing both types.

from rdkit.Chem import Descriptors, MolFromSmiles
import pandas as pd
import numpy as np

mols = [MolFromSmiles(s) for s in smiles['smiles'].tolist()]
descs = [(Descriptors.MolWt(m), Descriptors.NumHDonors(m)) for m in mols]
targets_mop = np.array(descs).reshape(-1,2)

from chemml.chem import RDKitFingerprint

rdkfp = RDKitFingerprint(fingerprint_type='morgan', radius=2, n_bits=1024)
features_mop = rdkfp.represent(smiles['smiles'].tolist())

Generating Morgan Fingerprints: 100%|██████████| 500/500 [00:00<00:00, 16302.11molecule/s]

Multioutput method 1: Scikit-learn MultiOutputRegressor wrapper

from sklearn.multioutput import MultiOutputRegressor

def multi_obj(individual):
    mop = MultiOutputRegressor(KernelRidge(alpha=individual[0], kernel=individual[1], degree=individual[2]))
    mop.fit(features_mop[:400], targets_mop[:400])
    pred = mop.predict(features_mop[400:])
    mae_molwt = regression_metrics(targets_mop[400:,0],pred[:,0])['MAE'].values[0]
    mae_hbd = regression_metrics(targets_mop[400:,1],pred[:,1])['MAE'].values[0]
    return np.mean([mae_molwt, mae_hbd])

ga = GeneticAlgorithm(evaluate=multi_obj, space=space, fitness=("min", ),
                    pop_size = 8, crossover_size=6, mutation_size=2, algorithm=3)
fitness_df, final_best_hyperparameters = ga.search(n_generations=5)

fitness_df

	Best_individual	Fitness_values	Time (hours)
0	(4.3428571428571425, linear, 2)	9.810322	0.000120
1	(3.7428511684765304, linear, 1)	9.792120	0.000098
2	(3.7428511684765304, linear, 1)	9.792120	0.000106
3	(3.7428511684765304, linear, 1)	9.792120	0.000097
4	(3.7428511684765304, linear, 1)	9.792120	0.000097

final_best_hyperparameters

{'alpha': np.float64(3.7428511684765304), 'kernels': 'linear', 'degree': 1}

Multioutput method 2: Neural networks (ChemML MLP)

from chemml.models import MLP

mlp_space = (
    {'nneurons_0'           : {'choice' : np.arange(10, 100, 10)}},
    {'nneurons_1'           : {'choice' : np.arange(10, 100, 10)}},
    {'nneurons_2'           : {'choice' : np.arange(10, 100, 10)}},
    {'activation_0'          : {'choice' : ['ReLU', 'Tanh']}},
    {'activation_1'          : {'choice' : ['ReLU', 'Tanh']}},
    {'activation_2'          : {'choice' : ['ReLU', 'Tanh']}},
    {'alpha'                : {'uniform': (0.0001, 0.1), 'mutation': (0,1)}}
)

def multi_obj_mlp(individual):

    layers = [individual[i] for i in range(0,3) if individual[i] != 0]
    activations = [individual[i] for i in range(3,6) if individual[i-3] != 0]

    mlp = MLP(nneurons=layers, activations=activations, alpha=individual[6],
              engine='pytorch', nfeatures=features_mop.shape[1], noutputs=targets_mop.shape[1])
    mlp.fit(features_mop[:400], targets_mop[:400])
    pred = mlp.predict(features_mop[400:])
    mae_molwt = regression_metrics(targets_mop[400:,0],pred[:,0])['MAE'].values[0]
    mae_hbd = regression_metrics(targets_mop[400:,1],pred[:,1])['MAE'].values[0]
    return np.mean([mae_molwt, mae_hbd])

ga = GeneticAlgorithm(evaluate=multi_obj_mlp, space=mlp_space, fitness=("min", ),
                    pop_size = 8, crossover_size=6, mutation_size=2, algorithm=3)
fitness_df, final_best_hyperparameters = ga.search(n_generations=5)

fitness_df

	Best_individual	Fitness_values	Time (hours)
0	(40, 10, 30, ReLU, Tanh, Tanh, 0.0143714285714...	20.587449	0.000793
1	(40, 10, 30, ReLU, Tanh, Tanh, 0.0143714285714...	20.587449	0.000745
2	(60, 40, 20, ReLU, Tanh, Tanh, 0.0018723872531...	20.571680	0.000809
3	(50, 70, 30, ReLU, Tanh, Tanh, 0.0144638850443...	20.563786	0.000718
4	(50, 70, 30, ReLU, Tanh, Tanh, 0.0144638850443...	20.563786	0.000741

final_best_hyperparameters

{'nneurons_0': np.int64(50),
 'nneurons_1': np.int64(70),
 'nneurons_2': np.int64(30),
 'activation_0': 'ReLU',
 'activation_1': 'Tanh',
 'activation_2': 'Tanh',
 'alpha': np.float64(0.014463885044310704)}