Check how GPy handles multi-task learning and prediction.

[jsjol]

cf: http://nbviewer.jupyter.org/github/SheffieldML/notebook/blob/master/GPy/coregionalized_regression_tutorial.ipynb

Also look at: GPy.models.gp_kronecker_gaussian_regression.py

[jsjol]

I think coregionalized regression is not what we want to do. Instead its seems much simpler to create a kernel using the active_dims argument. Example: k1 = GPy.kern.Linear(input_dim=1, active_dims=[0]) # works on the first column of X, index=0 k2 = GPy.kern.ExpQuad(input_dim=1, lengthscale=3, active_dims=[1]) # works on the second column of X, index=1 k = k1 * k2

[jsjol]

Computationally, it would be super convenient if the spatial covariance can be factorized like e.g. the RBF kernel. Then the covariance matrix can be written as Kronecker products of the respective coordinates covariance functions. See Wilson (2014) p. 63 for the theory. Check if the code in GPy/kern/src/grid_kerns.py is applicable.

[jsjol]

The method of Stegle et al. (2011), implemented in GPy.models.gp_kronecker_gaussian_regression.py is the way to go. The main computation is symmetric eigenvalue decomposition of spatial covariance and feature covariance, treated separately.

It will work out-of-the-box for arbitrary spatial covariance matrix. Very minor adjustments are required if you want to leverage special structure (which is of special interested for the spatial covariance). For example,

• Factorized kernel: eigendecomposition respects Kronecker structure (cf. Wilson (2014)). So, the eigendecomposition could be computed for each dimension separately.
• Sparse inverse covariance: model the inverse covariance directly. Compute the eigendecomposition of the inverse covariance using methods for sparse eigendecomposition (e.g. Arnoldi method).

[jsjol]

Now I think that to utilize the Kronecker structure, it would be more appropriate to use one of GPy's methods for gridded regression.

[jsjol]

I believe this is resolved by my extensions to GPy.