Setting up multi-output GP with partial target-correlations

[lucky-luke-98]

Hello,

I have an issue regarding the creation of a Gaussian Process for a multi-output experiment.
In my case, I have a block-design setup for the descriptives. Futher, I have 10 noisy targets.

Due to the nature of my targets, I know that some of them are correlated to another, and some of them are not.
In fact, I am assuming that the targets can be grouped into correlation groups: [[Y1, Y2, Y3], [Y4, Y5, Y6, Y7], [Y8], [Y9], [Y10]], where the grouping means that they are correlated to each other (in accordance to a given correlation threshold). This means that, e.g., target 10 does not correlated to other targets.

In botorch's issues and discussions, I am always reading that taking a SingleTaskGP is the go-to way to model Multi-Output GPs that do not have correlating targets. Further, targets that do correlate with each other should be modeled using a MultiTaskGP.

My idea was to separate the creation of each correlation group and treat them based on having correlation to other targets (-> MultiTaskGP) or not having correlation to other targets (-> SingleTaskGP).
Afterwards, to create one unified model, I would append them in a ModelListGP.

However, I am not sure about this procedure, since the documentation of ModelListGP mentions, that it is a useful tool for handing the grouping of non-block-design models and that block-design models should be handled in a SingleTaskGP or MultiTaskGP.
Further, following my code, I run into problems with this procedure when performing a Multi-Step-Lookahead (qMultiStepLookahead). Here, the fantasize method of the ModelListGP assumes that the number of outputs (self.num_outputs) is equal to the number of models in the ModelList. In my case this assumption is not true.
So, I am assuming that there is another, perhaps easier way to create the GP for my experiment.

I would be thankful for any recommendations!
Thanks in advance!

[esantorella]

When you say correlated, do you mean correlation in the underlying function value f, in the noise term, or in both? That is, can you run one evaluation that measures both Y1 and Y2 and produces correlated measurement errors, which would be different if you repeated the measurement?

[lucky-luke-98]

I only have correlation in the underlying function. The noise is of statistic nature which I tried to model using a fixed noise tensor as train_yvar argument. Similar to the deprecated FixedNoiseGP.

[esantorella]

Ah. You may have discovered a bug in the fantasize method. Your approach of using MultiTaskGPs and SingleTaskGPs within a ModelListGP sounds right to me. There's nothing wrong with using a ModelListGP with block design; it's more computationally efficient to use a SingleTaskGP when that's feasible, but in your case it probably isn't. And ModelListGP is supposed to support multi-task models as of #2154. However, it doesn't look like that PR updated or tested the fantasize method. Could you post a reproducible example?

[lucky-luke-98]

Thanks for your help. Great to have my design choices confirmed by you! I will send a reproducible example by the upcoming Sunday.

[esantorella]

Thanks! I opened #2398 for this so we can properly track it as a bug.

[Balandat]

As an aside: Another option could be to use a full multi-task model, but to encode the correlation structure in the cross-task covariance matrix parameterization where we force the elements that we know are zero to zero. This would require some thinking about how to best do it, and it also likely would require using some iterative structure-exploiting solves as the overall covariance matrix could become huge. The main benefit of it would be simplicity - no handling of multiple models or anything like that.

[lucky-luke-98]

Since I have a strict deadline with my project, I am thankful for every alternative suggestion.
This idea sounds interesting. However, I do not know what you mean by using an iterative structure. Could you elaborate on this?

My first intuition would be to utilize the gpytorch.priors.MultivariateNormalPrior, create the correlation matrix of my targets, set coefficients under a given threshold to 0, and initialize the MultivariateNormalPrior given the adapted correlation matrix. Then I can hand this Prior to the MultiTaskGP and proceed with this model to the MultiStepLookahead.

Does that sound reasonable to you?

[lucky-luke-98]

FYI:

I followed this approach leading to a performance boost in my prediction accuracy.
However, setting up a CustomPrior for the MultiTaskGP also leads to a vast increase in computational time for the fitting (at least for the MultivariateNormalPrior).

[...], and it also likely would require using some iterative structure-exploiting solves as the overall covariance matrix could become huge.

Would be interested on how to improve the computational time in an iterative manner.
My current approach looks something like this:

# computation of correlation between targets
real_corr_matrix = torch.corrcoef(train_y.T)
# restricting coefficient under a specific threshold to 0
altered_correlation_matrix = replace_tensor_entries(input=real_corr_matrix, threshold=0.5, replace_value=0.0)

# creation of custom prior
custom_prior = MultivariateNormalPrior(
    loc=torch.zeros(num_tasks), 
    covariance_matrix=altered_correlation_matrix
)

# creation of model
model = MultiTaskGP(
    train_X=train_X_augmented,
    train_Y=train_Y,
    task_feature=-1,
    train_Yvar=train_yvar,
    task_covar_prior=custom_prior
)
mll = ExactMarginalLogLikelihood(likelihood=model.likelihood, model=model)

model.train()
mll.train()

# model fitting
fit_gpytorch_mll(mll, max_attempts=30)

replace_tensor_entries() is a custom function that replaces values that are under a given absolute threshold with the replace_value.