QM9 kernel models using molSimplify, RACs and R: Part 2

Tuesday, February 20, 2018
Jon Paul Janet

In the second part of our tutorial, we will demonstrate how to use R to conduct kernel based prediction of atomization energies based on RACs. You’ll need QM9_descriptor_file.csv, which we prepared using molSimplify in the previous tutorial and also provide here. We’ll use R to conduct a simple KRR model using a radial basis function kernel – we recommend the excellent, free RStudio IDE but a similar procedure could be followed in your favorite language. In particular, we’ll use the CVST package to train our model and DRR to conduct accelerated cross validation. Both packages are available from CRAN and depend on kernlab.

We’ll use 4000 points for training and the remaining ~131k molecules as a test set and construct a model the following general form:

We use the same, slightly uncommon form of the kernel (sigma being the inverse width) as in the kernlab implementation. The sum here runs over the n=4000 training examples. Effectively, we are performing a non-linear transformation from our space of descriptors into to a high-(in this case, infinite-)dimensional space of transformed features. The kernel function returns the inner product in this space, and as can be seen from the above, the influence of each training points is controlled by the Euclidean distance between the training point and the query point. The rate of influence decay with distance is controlled by the hyperparameter σ. The weights α are learned during training, where we seek to minimize a loss function along with a regularization term controlled by λ:

While the intention of this guide is not to review KRR methods, we'll make a brief comment that the model complexity is controlled by n, the number of training points used, and solving for the model coefficients typically involves the (approximate/blockwise) inversion of the nxn kernel matrix, and so in general shows poor scaling with large numbers of training examples.

Again, we have provided a full script that will allow you conduct the calculation (simpleKRR.R), and we’ll skip some of the boilerplate and focus on the key lines. The first step is to read in the data and partition into test and training sets:

## load the file prepared in part 1
data <- read.csv("QM9_descriptor_file.csv",header=TRUE,row.names = 1)
## choose number of training points
training_examples <- 4000
## partition the data
train_ind = sample(nrow(data),training_examples)
df_test = data[-train_ind,]
df_train=data[train_ind,]

The next step is to apply normalization to the data. We’ll follow best practice and normalize (to mean 0 and variance 1) using the training data only:

## normalize data 
df_train <- scale(df_train,center=TRUE,scale=TRUE)
train_center <- attr(df_train, 'scaled:center')
train_scale <- attr(df_train, 'scaled:scale')
df_train <- as.data.frame(df_train)
# note: best practice, only use training scales!
df_test <- scale(df_test,center = train_center,scale=train_scale)

The next step is to convert the data to a format CVST understands, and then build the test and training objects:

## Load into CVST
x_train <- df_train[,2:ncol(df_train)]
y_train <-  df_train[,'u0']
x_test <- df_test[,2:ncol(df_test)]
y_test <-  df_test[,'u0']
train_data_CVST <-constructData(x=as.matrix(x_train), y=y_train)
test_data_CVST <-constructData(x=as.matrix(x_test), y=y_test)

For our first try, we won’t worry about optimizing hyperparameters (σ, λ) and instead simply use some reasonable values:


## set up learner
krr_learner = constructFastKRRLearner()   ## Build the base learner
## first time around, don't do CV but use a reasonable kernel:
params_krr = list(kernel="rbfdot", sigma=1E-6, lambda=1E-11,nblocks=4) 

The constructFastKRRLearner() comes from the DRR package and implements fast, divide-and-conquer matrix inversion that will be useful for doing cross validation a little later and is critical if we want to use larger kernel matrices. Now, we can train and test the model. We will convert back to original units (Ha) and then further convert units to kcal/mol, and print the results:


## train model
krr_trained = krr_learner$learn(train_data_CVST, params_krr)
## use model for prediction
test_krr = krr_learner$predict(krr_trained, test_data_CVST)
train_krr = krr_learner$predict(krr_trained, train_data_CVST)
pred_test_y_krr =  as.vector(test_krr*train_scale["u0"] +  train_center["u0"])
pred_train_y_krr = as.vector(train_krr* train_scale["u0"] +  train_center["u0"])
pred_test_y_krr = as.vector(test_krr*train_scale["u0"] + train_center["u0"])
pred_train_y_krr = as.vector(train_krr*train_scale["u0"] + train_center["u0"])
test_target = as.vector(test_data_CVST$y*train_scale["u0"]  +  train_center["u0"])
train_target = as.vector(train_data_CVST$y*train_scale["u0"]  +  train_center["u0"])
## evaluate
train_mse_krr = mse(pred_train_y_krr,train_target)*(HA_to_kcal_mol^2)
test_mse_krr = mse(pred_test_y_krr,test_target)*(HA_to_kcal_mol^2)
train_e_krr =  (pred_train_y_krr - train_target)*HA_to_kcal_mol
test_e_krr =  (pred_test_y_krr - test_target)*HA_to_kcal_mol

Running the code produces the followings output:


krr train rmse = 10.36 kcal/mol       
krr train MAE = 7.37 kcal/mol       
krr test rmse = 10.98 kcal/mol       
krr test MAE =  7.79 kcal/mol 

Note that this corresponds to one random sample of 4000 training points and the results might differ a little bit each time the code is run. It would be best practice to run this code multiple times and average the results.

In order to improve the result using cross-validation, we will do a quick 5x5 grid search in (σ, λ) around our initial guess, using the following code:

## set parameters for CV
params_CVST_krr = constructParams(kernel = 'rbfdot',sigma=10^(seq(-8,-4,length=5)),lambda=10^(seq(-12,-8,length=5)),nblocks=4)
## conduct 10-fold CV
CV_ret_krr <-  CV(data= train_data_CVST, learner = krr_learner,params =  params_CVST_krr,fold = 10, verbose = TRUE)
## capture the optimal parameter choices
params_krr <- CV_ret_krr[[1]]

These lines define our grid of hyperparameter values and conduct 10-fold cross-validation to choose the optimal combination. In general a much larger field of hyperparameter values should be considered and the grid should be iteratively refined, but we have preemptively narrowed down the range in order to reduce run time. Using the new params_krr in the code above (the provided script does both) yields the following results for a typical run:

krr train rmse = 9.53 kcal/mol       
krr train MAE = 6.68 kcal/mol       
krr test rmse = 10.42 kcal/mol       
krr test MAE =  7.20 kcal/mol 

We will warn that conducting detailed cross-validation is expensive, and adding the above lines of codes decreases the run time from seconds to minutes on a standard workstation.

Scripts: 

About Us

The Kulik group focuses on the development and application of new electronic structure methods and atomistic simulations tools in the broad area of catalysis.

Our Interests

We are interested in transition metal chemistry, with applications from biological systems (i.e. enzymes) to nonbiological applications in surface science and molecular catalysis.

Our Focus

A key focus of our group is to understand mechanistic features of complex catalysts and to facilitate and develop tools for computationally driven design.

Contact Us

Questions or comments? Let us know! Contact Dr. Kulik: