Need help in modeling 5D input, multi-output GP

Questions version agnostic
1

Hi,

Let me start by saying I am very new to GP modeling.

I have adapted the codes from the Multi-output Gaussian Processes Notebook to my own system. The ICM model in the notebook worked on my data for 1 batch of the data set. However, I need help in modifying the model when I want to train it on 6 batches of data. I’m guessing the modification has to do with the shapes of the parameters.

Each batch has an output of y_species(t), for 11 time points and for 6 species.

Here is what my input array looks like for all 6 Batches:

# First column
time = np.tile([1, 15, 30, 45, 60, 120, 180, 300, 420, 720, 1080], 36)

# Second, third, and fourth columns. Features related to each batch.
X1 = np.repeat([0.427, 0.427, 0.639, 0.249, 0.860, 0.427], 66)
X2 = np.repeat([0.389, 0.389, 0.389, 0.227, 0.782, 0.389], 66)
X3 = np.repeat([0.467, 0.542, 0.467, 0.273, 0.935, 0.427], 66)

# Fifth column
speciesIDX = np.tile(np.repeat([0, 1, 2, 3, 4, 5], 11), 6)

# Combine into a single array
input_array = np.column_stack((time, X1, X2, X3, speciesIDX))

The following model works when it is trained on the Batch 1 data (first 66 elements):

def get_icm(input_dim, kernel, W=None, kappa=None, B=None, active_dims=None):
    """
    This function generates an ICM kernel from an input kernel and a Coregion kernel.
    """
    coreg = pm.gp.cov.Coregion(input_dim=input_dim, W=W, kappa=kappa, B=B, active_dims=active_dims)
    icm_cov = kernel * coreg  # Use Hadamard Product for separate inputs
    return icm_cov

X = input_array[:66,:]
Y = output[:66]

with pm.Model() as model:
    # Priors
    ell = pm.Gamma("ell", alpha=5, beta=2)
    eta = pm.Gamma("eta", alpha=0.5, beta=0.1)
    kernel = eta**2 * pm.gp.cov.ExpQuad(input_dim=5, ls=ell, active_dims=[0])
    sigma = pm.HalfNormal("sigma", sigma=0.06)

    # Get the ICM kernel
    W = pm.Normal("W", mu=0, sigma=1, shape=(6, 2), initval=np.random.randn(6, 2))
    kappa = pm.Gamma("kappa", alpha=3.0, beta=2.0, shape=6)
    B = pm.Deterministic("B", pt.dot(W, W.T) + pt.diag(kappa))
    cov_icm = get_icm(input_dim=5, kernel=kernel, B=B, active_dims=[4])

    # Define a Multi-output GP
    mogp = pm.gp.Marginal(cov_func=cov_icm)
    y_ = mogp.marginal_likelihood("f", X, Y, sigma=sigma)

How do I modify this model when I want to train it on 6 batches of data?

Many thanks,

2

What do you mean by batch exactly? Is each batch independent, or is there some relation between them?

3

Hi,

I have 6 batches of data. Each data is the concentration at a particular time: Ci(t) for i = 1, 2, …, 6 species, with 11 time points. I have successfully made a multi-output GP for 1 batch of data just by using the following input:

time = np.tile([1, 15, 30, 45, 60, 120, 180, 300, 420, 720, 1080], 6) # for 6 species
speciesIDX = np.repeat([0, 1, 2, 3, 4, 5], 11) #species index
input_array = np.column_stack((time, speciesIDX))

Now, I have X1, X2, and X3 which has unique values per batch. You can think of them as additional features that describe the batch data. How do I combine this multi-ouput GP with a multi dimension input? I have been reading about the Kronecker product. Do you think I can use that here?

My ultimate goal is to have some prediction of the Ci(t) for i = 1, 2, …, 6 species, given X1, X2, X3.

4

Update: I created a covariance per dimension and set the overall covariance as the overall product.

with pm.Model() as model:
    # Priors
    ell_t = pm.Gamma("ell", alpha=5, beta=2)
    ls1 = pm.Gamma("ls1", alpha=0.2, beta=0.2)
    ls2 = pm.Gamma("ls2", alpha=0.2, beta=0.2)
    ls3 = pm.Gamma("ls3", alpha=0.2, beta=0.2)
    
    eta = pm.Gamma("eta", alpha=0.5, beta=0.1)
    kernel = eta**2 * pm.gp.cov.ExpQuad(input_dim=5, ls=ell_t, active_dims=[0])
    
    # Get the ICM kernel
    W = pm.Normal("W", mu=0, sigma=1, shape=(6, 2), initval=np.random.randn(6, 2))
    kappa = pm.Gamma("kappa", alpha=3.0, beta=2.0, shape=6)
    B = pm.Deterministic("B", pt.dot(W, W.T) + pt.diag(kappa))
    coreg = pm.gp.cov.Coregion(input_dim=5, B=B, active_dims=[4])

    cov_x1 = pm.gp.cov.Matern52(5, ls=ls1, active_dims=[1]) 
    cov_x2 = pm.gp.cov.Matern52(5, ls=ls2, active_dims=[2])
    cov_x3 = pm.gp.cov.Matern52(5, ls=ls3, active_dims=[3])
    cov_product = kernel*cov_x1*cov_x2*cov_x3*coreg
    
    sigma = pm.HalfNormal("sigma", sigma=0.06)

    # Define a Multi-output GP
    mogp = pm.gp.Marginal(cov_func=cov_product)
    y_ = mogp.marginal_likelihood("f", X, Y, sigma=sigma)

I will report the results as soon as I get them. The sampling seems to take a while. Any suggestions on making this faster?

5

Sorry it took a while to get back to you. This is a tricky question. I think your model is going to be a little slow inherently, but I would definitely recommend moving away from those Gamma(0.2, 0.2) priors. I’d recommend using InverseGamma priors as per this post. You can set the upper and lower bounds using pm.find_constrained_prior or the maxent funcion in Preliz so you don’t have to choose alpha and beta directly. The lengthscale priors can have a very big effect on sampling speed.

1 Like