Currently this package eagerly constructs and uses a dense kernel matrix. However, iterative GP methods are very useful when the data is large enough that the kernel matrix can't even fit into memory. It would be better to construct and use a lazy representation of the kernel matrix.
I've started work on adding this functionality to the KernelFunctions API: JuliaGaussianProcesses/KernelFunctions.jl#515