Multi-output Gaussian Process in GPFlow

Question

I need to implement a Multi-Output GP that works of batch inputs. I have time series where:

• X has shape (6,100,1), so 6 batches each of shape (100,1), simply values between 0 and 1
• Y has shape (6,100,21) So i need to work with 21 outputs.
  I followed this code of the GPFlow tutorial:

import numpy as np
import gpflow

P = L = 21
M = 10

# Define a kernel suitable for batched input
kern_list = [gpf.kernels.SquaredExponential() for _ in range(L)]
kernel = gpf.kernels.LinearCoregionalization(kern_list, W=np.random.randn(P, L))

# Define inducing points
Zinit = np.linspace(0,1,M)[:,None]  # Assuming you have Zinit defined
Z = Zinit.copy()
iv = gpf.inducing_variables.SharedIndependentInducingVariables(
    gpf.inducing_variables.InducingPoints(Z)
)

# LMC initialization
q_mu = np.zeros((M, L))
# initialize \sqrt(Σ) of variational posterior to be of shape LxMxM
q_sqrt = np.repeat(np.eye(M)[None, ...], L, axis=0) * 1.0

# Create an SVGP model
model = gpflow.models.SVGP(kernel=kernel, 
                           likelihood=gpflow.likelihoods.Gaussian(),
                           inducing_variable=iv,
                           q_mu=q_mu, q_sqrt=q_sqrt, num_latent_gps=21)

# Optimize the model
def optimize_model_with_scipy(model):
    optimizer = gpf.optimizers.Scipy()
    optimizer.minimize(
        model.training_loss_closure((X,Y)),
        variables=model.trainable_variables,
        method="l-bfgs-b",
        options={"disp": True, "maxiter": 1000},
    )

optimize_model_with_scipy(model)

The code runs but it's not fitting the data... I'm wondering if it's a problem of initialization