I've set up a multi task Gaussian process problem with the help of coregionalization (more precise: an intrinsic coregionalization model (ICM)). My features space constists of ~20 input dimensions for six outputs. Following this paper the kernel matrix $\mathbf{K}$ corresponding to a dataset $\mathbf{X}$ takes the form
$$\mathbf{K}(\mathbf{X, X}) = \mathbf{B} \otimes k(\mathbf{X, X})$$
$k(\mathbf{X, X})$ is here a GP kernel of choice possibly able to do Automatic Relevance Determination (ARD). In my case, $k(\mathbf{X, X})$ is a rational quadratic kernel, exploiting its lengthscale for ARD. $\mathbf{B}$ is the coregionalization matrix, $\otimes$ the kronecker product.
According to the paper, I understand that $k(\mathbf{X, X})$ models the input dependence, independently of the particular set of outputs, while $\mathbf{B}$ models the dependence between the outputs, indepently of the inputs.
Training the model and modelling the outputs works. I am now interested in the feature relevance / input sensitivity of the input features for each output individually. I assume that the lengthscales of $k(\mathbf{X, X})$, after training, somehow reflect the input feature importance for all outputs together.
I am wondering if there is a way to combine the information carried by $\mathbf{B}$ with the info carried by $k(\mathbf{X, X})$ to obtain a set of most important features for each output individually?
(The technical implementation is done with the python GPy.models.GPCoregionalizedRegression class of package GPy).
Best Answer
Unless I misunderstand, with this type of model the implicit assumption is that all of your outputs are determined by the same "underlying" or "latent" Gaussian Process.
The "relevance" of the input features should consequently be identical for each output, since the underlying covariance structure with respect to $X$ will be the same (the latent GP).