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

I have a lot of test curves and I want to optimize the length and scale parameters simultaneously for all curves. Is this possible?

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.