Estimating hyper-parameters of a Gaussian Process

Hello all,

I am not sure if this is even doable, so I thought of asking the community to start with. I've been working with MATLAB's fitrgp command to fit the data using a GP. I have multiple GPR objects (total of M), each with their own observations, and the hyper-parameters for each GP are estimated independently from each other. My understanding is that the algorithm searches for the set of hyper-parameters that maximize the log-marginal likelihood log p(yi|Xi,θi), where Xi,yi is the training data set for the i- th GP.

Ideally, I want to train these GPs jointly. Namely, I want to find the set of hyper-parametrs that minimize the sum of log-marginal likelihoods, θ = arg max θ \sum_(i=1)^(M) [log p(yi|Xi,θ)]. Is it possible to actually change this in MATLAB's fitrgp functions?

Any help would be much appreciated!