How to optimize Gaussian-process parameters for multiple tasks with GPML?

You can write a likelihood for the data in a usual way as far as I can understand.

Let we have data $D = {X, Y} = {mathbf{x}_i, mathbf{y}_i}$_{i = 1}^{n} with $mathbf{x}_i in mathbb{R}^p$, and outputs $mathbf{y}_i in mathbb{R}^m$ is multidimensional. Denote by $mathbf{y}^j = {y_{1j}, y_{2j}, ldots, y_{nj}}$ - A vector which corresponds to $j$-th curve.

You suppose that each curve is a realization of a Gaussian process, but all Gaussian processes have zero mean value and the same covariance function $k(mathbf{x}_i, mathbf{x}_j)$, so you get for each curve the loglikelihood of the the form: $$ L(X, mathbf{y}^j) = -frac12 ( (mathbf{y}^j)^T K^{-1} mathbf{y}^j + |K| + n log (2 pi)), $$ here $K = {k(mathbf{x}_i, mathbf{x}_j)}_{i, j = 1}^{n}$ is a sample covariance matrix. And the joint likelihood is a product of separate likelihoods, so we can use common techniques to optimize it. For example, for most covariance function we can calculate derivatives w.r.t covariance function parameters and use a derivative-based optimization.