Can Gaussian Process be fit incrementally? I am using Gaussian Process Regressor to fit data for a Bayesian Optimiser. This is a relevant part of my Python code.

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import (RBF, Matern, RationalQuadratic,
                                          ExpSineSquared, DotProduct,
                                          ConstantKernel)
kernel = Matern(length_scale=lenSc,nu=2.5)
gp = GaussianProcessRegressor(kernel=kernel)
l2=10
Xi=np.concatenate((nprand.uniform(10.0,30.0,[l2,6]),nprand.uniform(-0.3,0.3,[l2,42]),nprand.uniform(0.0,30.0,[l2,1]),nprand.uniform(-30.0,30.0,[l2,1]),nprand.uniform(0.0,0.3,[l2,6]),nprand.uniform(5.0,20.0,[l2,6])),axis=1)
yi=map(lambda x: f1(x), Xi)
gp.fit(Xi, yi)
while(1):
   #Some code to get the new values- cc and nv
   Xi=numpy.row_stack((Xi,cc))
   yi=np.concatenate((yi,[nv]))
   gp.fit(Xi, yi)


Each iteration of the Bayesian Optimiser, I add a new element to $Xi$ and $yi$ and fit $gp$, the Gaussian Process, again. As the size of $Xi$ and $yi$ reach over 200, the time taken to fit the data becomes noticeable, this decreases the overall efficiency of my Bayesian Optimiser. 

The Bayesian Optimiser has to make around 1000 iterations, each iteration needs to fit the data all over again. So, I was wondering if there's a way to incrementally fit the data? If that's not already being done by sklearn's GaussianProcessRegressor.

Thank you!
I think with GPR, the computational cost increases on an order of N^3, so the amount of time with more data increases drastically.

I don't think the current version of sklearn's GPR implements the kind of incremental fit you are expecting. With larger data and iterations, I suggest implementing your own version on a GPU or somethings since GPR implementation is not that difficult.