merged

GPy/core/sparse_gp.py

@@ -31,7 +31,7 @@ class SparseGP(GP):
 
     """
 
-    def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, name='sparse gp'):
+    def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, name='sparse gp', Y_metadata=None):
 
         #pick a sensible inference method
         if inference_method is None:
@@ -45,7 +45,7 @@ class SparseGP(GP):
         self.Z = Param('inducing inputs', Z)
         self.num_inducing = Z.shape[0]
 
-        GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name)
+        GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata)
 
         self.add_parameter(self.Z, index=0)
 
@@ -53,7 +53,7 @@ class SparseGP(GP):
         return isinstance(self.X, VariationalPosterior)
 
     def parameters_changed(self):
-        self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y)
+        self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y, self.Y_metadata)
         self.likelihood.update_gradients(self.grad_dict['dL_dthetaL'])
         if isinstance(self.X, VariationalPosterior):
             #gradients wrt kernel
@@ -75,7 +75,6 @@ class SparseGP(GP):
             target += self.kern.gradient
             self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None)
             self.kern.gradient += target
-
             #gradients wrt Z
             self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z)
             self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X)

GPy/examples/coreg_example.py

@@ -2,6 +2,55 @@ import numpy as np
 import pylab as pb
 import GPy
 pb.ion()
+pb.close('all')
+
+X1 = np.arange(3)[:,None]
+X2 = np.arange(4)[:,None]
+I1 = np.zeros_like(X1)
+I2 = np.ones_like(X2)
+
+_X = np.vstack([ X1, X2 ])
+_I = np.vstack([ I1, I2 ])
+
+X = np.hstack([ _X, _I ])
+
+Y1 = np.sin(X1/8.)
+Y2 = np.cos(X2/8.)
+
+Bias = GPy.kern.Bias(1,active_dims=[0])
+Coreg = GPy.kern.Coregionalize(1,2,active_dims=[1])
+K = Bias.prod(Coreg,name='X')
+
+#K.coregion.W = 0
+#print K.coregion.W
+
+#print Bias.K(_X,_X)
+#print K.K(X,X)
+
+#pb.matshow(K.K(X,X))
+
+
+Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")]
+kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=2,kernels_list=Mlist,name='H')
+kern.B.W = 0
+kern.B.kappa = 1.
+#kern.B.W.fix()
+#kern.B.kappa.fix()
+#m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern)
+m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], kernel=kern)
+#m.optimize()
+m.checkgrad(verbose=1)
+
+fig = pb.figure()
+ax0 = fig.add_subplot(211)
+ax1 = fig.add_subplot(212)
+slices = GPy.util.multioutput.get_slices([Y1,Y2])
+m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
+#m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
+
+
+
+"""
 
 X1 = 100 * np.random.rand(100)[:,None]
 X2 = 100 * np.random.rand(100)[:,None]
@@ -28,3 +77,4 @@ slices = GPy.util.multioutput.get_slices([Y1,Y2])
 m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
 m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
 
+"""

GPy/inference/latent_function_inference/dtc.py

@@ -91,8 +91,12 @@ class vDTC(object):
     def __init__(self):
         self.const_jitter = 1e-6
 
-    def inference(self, kern, X, Z, likelihood, Y):
+    def inference(self, kern, X, X_variance, Z, likelihood, Y):
-        #assert X_variance is None, "cannot use X_variance with DTC. Try varDTC."
+        assert X_variance is None, "cannot use X_variance with DTC. Try varDTC."
+
+        #TODO: MAX! fix this!
+        from ...util.misc import param_to_array
+        Y = param_to_array(Y)
 
         num_inducing, _ = Z.shape
         num_data, output_dim = Y.shape
@@ -105,14 +109,15 @@ class vDTC(object):
         Kmm = kern.K(Z)
         Knn = kern.Kdiag(X)
         Knm = kern.K(X, Z)
-        KnmY = np.dot(Knm.T,Y)
+        U = Knm
+        Uy = np.dot(U.T,Y)
 
         #factor Kmm 
         Kmmi, L, Li, _ = pdinv(Kmm)
 
         # Compute A
-        LiKmnbeta = np.dot(Li, Knm.T)*np.sqrt(beta)
+        LiUTbeta = np.dot(Li, U.T)*np.sqrt(beta)
-        A_ = tdot(LiKmnbeta)
+        A_ = tdot(LiUTbeta)
         trace_term = -0.5*(np.sum(Knn)*beta - np.trace(A_))
         A = A_ + np.eye(num_inducing)
 
@@ -120,7 +125,7 @@ class vDTC(object):
         LA = jitchol(A)
 
         # back substutue to get b, P, v
-        tmp, _ = dtrtrs(L, KnmY, lower=1)
+        tmp, _ = dtrtrs(L, Uy, lower=1)
         b, _ = dtrtrs(LA, tmp*beta, lower=1)
         tmp, _ = dtrtrs(LA, b, lower=1, trans=1)
         v, _ = dtrtrs(L, tmp, lower=1, trans=1)
@@ -140,18 +145,19 @@ class vDTC(object):
         LAL = Li.T.dot(A).dot(Li)
         dL_dK = Kmmi - 0.5*(vvT_P + LAL)
 
-        # Compute dL_dKnm
+        # Compute dL_dU
         vY = np.dot(v.reshape(-1,1),Y.T)
-        dL_dKmn = vY - np.dot(vvT_P - Kmmi, Knm.T)
+        #dL_dU = vY - np.dot(vvT_P, U.T)
-        dL_dKmn *= beta
+        dL_dU = vY - np.dot(vvT_P - Kmmi, U.T)
+        dL_dU *= beta
 
         #compute dL_dR
-        Knmv = np.dot(Knm, v)
+        Uv = np.dot(U, v)
-        dL_dR = 0.5*(np.sum(Knm*np.dot(Knm,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Knmv*Y, 1) + np.sum(np.square(Knmv), 1) )*beta**2
+        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) )*beta**2
         dL_dR -=beta*trace_term/num_data
 
         dL_dthetaL = likelihood.exact_inference_gradients(dL_dR)
-        grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn) + -0.5*beta, 'dL_dKnm':dL_dKmn.T, 'dL_dthetaL':dL_dthetaL}
+        grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn) + -0.5*beta, 'dL_dKnm':dL_dU.T, 'dL_dthetaL':dL_dthetaL}
 
         #construct a posterior object
         post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L)

GPy/inference/latent_function_inference/exact_gaussian_inference.py

@@ -52,6 +52,6 @@ class ExactGaussianInference(object):
 
         dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi)
 
-        dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK))
+        dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK),Y_metadata)
 
         return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL}

GPy/inference/latent_function_inference/var_dtc.py

@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 from posterior import Posterior
-from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify
+from ...util.linalg import mdot, jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify
 from ...util import diag
 from ...core.parameterization.variational import VariationalPosterior
 import numpy as np
@@ -74,7 +74,7 @@ class VarDTC(object):
         trYYT = self.get_trYYT(Y)
 
         # do the inference:
-        het_noise = beta.size < 1
+        het_noise = beta.size > 1
         num_inducing = Z.shape[0]
         num_data = Y.shape[0]
         # kernel computations, using BGPLVM notation
@@ -134,16 +134,16 @@ class VarDTC(object):
 
         # log marginal likelihood
         log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise,
-            psi0, A, LB, trYYT, data_fit)
+            psi0, A, LB, trYYT, data_fit, Y)
 
         #put the gradients in the right places
         dL_dR = _compute_dL_dR(likelihood,
             het_noise, uncertain_inputs, LB,
             _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A,
             psi0, psi1, beta,
-            data_fit, num_data, output_dim, trYYT)
+            data_fit, num_data, output_dim, trYYT, Y)
 
-        dL_dthetaL = likelihood.exact_inference_gradients(dL_dR)
+        dL_dthetaL = likelihood.exact_inference_gradients(dL_dR,Y_metadata)
 
         if uncertain_inputs:
             grad_dict = {'dL_dKmm': dL_dKmm,
@@ -387,7 +387,7 @@ def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, C
     return dL_dpsi0, dL_dpsi1, dL_dpsi2
 
 
-def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, psi0, psi1, beta, data_fit, num_data, output_dim, trYYT):
+def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, psi0, psi1, beta, data_fit, num_data, output_dim, trYYT, Y):
     # the partial derivative vector for the likelihood
     if likelihood.size == 0:
         # save computation here.
@@ -396,19 +396,20 @@ def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf,
         if uncertain_inputs:
             raise NotImplementedError, "heteroscedatic derivates with uncertain inputs not implemented"
         else:
-            from ...util.linalg import chol_inv
+            #from ...util.linalg import chol_inv
-            LBi = chol_inv(LB)
+            #LBi = chol_inv(LB)
+            LBi, _ = dtrtrs(LB,np.eye(LB.shape[0]))
+
             Lmi_psi1, nil = dtrtrs(Lm, psi1.T, lower=1, trans=0)
             _LBi_Lmi_psi1, _ = dtrtrs(LB, Lmi_psi1, lower=1, trans=0)
 
-            dL_dR = -0.5 * beta + 0.5 * likelihood.V**2
+            dL_dR = -0.5 * beta + 0.5 * (beta*Y)**2
             dL_dR += 0.5 * output_dim * (psi0 - np.sum(Lmi_psi1**2,0))[:,None] * beta**2
 
             dL_dR += 0.5*np.sum(mdot(LBi.T,LBi,Lmi_psi1)*Lmi_psi1,0)[:,None]*beta**2
 
-            dL_dR += -np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T * likelihood.Y * beta**2
+            dL_dR += -np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T * Y * beta**2
             dL_dR += 0.5*np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T**2 * beta**2
-
     else:
         # likelihood is not heteroscedatic
         dL_dR = -0.5 * num_data * output_dim * beta + 0.5 * trYYT * beta ** 2
@@ -416,11 +417,11 @@ def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf,
         dL_dR += beta * (0.5 * np.sum(A * DBi_plus_BiPBi) - data_fit)
     return dL_dR
 
-def _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit):
+def _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit,Y):
-#compute log marginal likelihood
+    #compute log marginal likelihood
     if het_noise:
-        lik_1 = -0.5 * num_data * output_dim * np.log(2. * np.pi) + 0.5 * np.sum(np.log(beta)) - 0.5 * np.sum(likelihood.V * likelihood.Y)
+        lik_1 = -0.5 * num_data * output_dim * np.log(2. * np.pi) + 0.5 * np.sum(np.log(beta)) - 0.5 * np.sum(beta * Y**2)
-        lik_2 = -0.5 * output_dim * (np.sum(beta * psi0) - np.trace(A))
+        lik_2 = -0.5 * output_dim * (np.sum(beta.flatten() * psi0) - np.trace(A))
     else:
         lik_1 = -0.5 * num_data * output_dim * (np.log(2. * np.pi) - np.log(beta)) - 0.5 * beta * trYYT
         lik_2 = -0.5 * output_dim * (np.sum(beta * psi0) - np.trace(A))

GPy/likelihoods/gaussian.py

@@ -50,7 +50,11 @@ class Gaussian(Likelihood):
         if isinstance(gp_link, link_functions.Identity):
             self.log_concave = True
 
-    def gaussian_variance(self, Y, Y_metadata=None):
+    def betaY(self,Y,Y_metadata=None):
+        #TODO: ~Ricardo this does not live here
+        return Y/self.gaussian_variance(Y_metadata)
+
+    def gaussian_variance(self, Y_metadata=None):
         return self.variance
 
     def update_gradients(self, grad):

GPy/likelihoods/mixed_noise.py

@@ -18,6 +18,17 @@ class MixedNoise(Likelihood):
         self.likelihoods_list = likelihoods_list
         self.log_concave = False
 
+    def gaussian_variance(self, Y_metadata):
+        assert all([isinstance(l, Gaussian) for l in self.likelihoods_list])
+        ind = Y_metadata['output_index'].flatten()
+        variance = np.zeros(ind.size)
+        for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))):
+            variance[ind==j] = lik.variance
+        return variance[:,None]
+
+    def betaY(self,Y,Y_metadata):
+        return Y/self.gaussian_variance(Y_metadata=Y_metadata)
+
     def update_gradients(self, gradients):
         self.gradient = gradients
 
@@ -32,13 +43,9 @@ class MixedNoise(Likelihood):
             _variance = np.array([self.likelihoods_list[j].variance for j in ind ])
             if full_cov:
                 var += np.eye(var.shape[0])*_variance
-                #d = 2*np.sqrt(np.diag(var))
-                #low, up = mu - d, mu + d
             else:
                 var += _variance
-                #d = 2*np.sqrt(var)
+            return mu, var
-                #low, up = mu - d, mu + d
-            return mu, var#, low, up
         else:
             raise NotImplementedError
 
@@ -51,12 +58,13 @@ class MixedNoise(Likelihood):
 
 
     def covariance_matrix(self, Y, Y_metadata):
-        assert all([isinstance(l, Gaussian) for l in self.likelihoods_list])
+        #assert all([isinstance(l, Gaussian) for l in self.likelihoods_list])
-        ind = Y_metadata['output_index'].flatten()
+        #ind = Y_metadata['output_index'].flatten()
-        variance = np.zeros(Y.shape[0])
+        #variance = np.zeros(Y.shape[0])
-        for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))):
+        #for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))):
-            variance[ind==j] = lik.variance
+        #    variance[ind==j] = lik.variance
-        return np.diag(variance)
+        #return np.diag(variance)
+        return np.diag(self.gaussian_variance(Y_metadata).flatten())
 
 
     def samples(self, gp, Y_metadata):

GPy/models/__init__.py

@@ -15,4 +15,5 @@ from mrd import MRD
 from gradient_checker import GradientChecker
 from ss_gplvm import SSGPLVM
 from gp_coregionalized_regression import GPCoregionalizedRegression
+from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression
 #.py file not included!!! #from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression

GPy/models/sparse_gp_coregionalized_regression.py

@@ -0,0 +1,66 @@
+# Copyright (c) 2012 - 2014 the GPy Austhors (see AUTHORS.txt)
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import numpy as np
+from ..core import SparseGP
+from ..inference.latent_function_inference import VarDTC
+from .. import likelihoods
+from .. import kern
+from .. import util
+
+class SparseGPCoregionalizedRegression(SparseGP):
+    """
+    Sparse Gaussian Process model for heteroscedastic multioutput regression
+
+    This is a thin wrapper around the SparseGP class, with a set of sensible defaults
+
+    :param X_list: list of input observations corresponding to each output
+    :type X_list: list of numpy arrays
+    :param Y_list: list of observed values related to the different noise models
+    :type Y_list: list of numpy arrays
+    :param Z_list: list of inducing inputs (optional)
+    :type Z_list: empty list | list of numpy arrays
+    :param kernel: a GPy kernel, defaults to RBF ** Coregionalized
+    :type kernel: None | GPy.kernel defaults
+    :likelihoods_list: a list of likelihoods, defaults to list of Gaussian likelihoods
+    :type likelihoods_list: None | a list GPy.likelihoods
+    :param num_inducing: number of inducing inputs, defaults to 10 per output (ignored if Z_list is not empty)
+    :type num_inducing: integer | list of integers
+
+    :param name: model name
+    :type name: string
+    :param W_rank: number tuples of the corregionalization parameters 'W' (see coregionalize kernel documentation)
+    :type W_rank: integer
+    :param kernel_name: name of the kernel
+    :type kernel_name: string
+    """
+
+    def __init__(self, X_list, Y_list, Z_list=[], kernel=None, likelihoods_list=None, num_inducing=10, X_variance=None, name='SGPCR',W_rank=1,kernel_name='X'):
+
+        #Input and Output
+        X,Y,self.output_index = util.multioutput.build_XY(X_list,Y_list)
+        Ny = len(Y_list)
+
+        #Kernel
+        if kernel is None:
+            kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name)
+
+        #Likelihood
+        likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list)
+
+        #Inducing inputs list
+        if len(Z_list):
+            assert len(Z_list) == self.output_dim, 'Number of outputs do not match length of inducing inputs list.'
+        else:
+            if isinstance(num_inducing,np.int):
+                num_inducing = [num_inducing] * Ny
+            num_inducing = np.asarray(num_inducing)
+            assert num_inducing.size == Ny, 'Number of outputs do not match length of inducing inputs list.'
+            for ni,Xi in zip(num_inducing,X_list):
+                i = np.random.permutation(Xi.shape[0])[:ni]
+                Z_list.append(Xi[i].copy())
+
+        Z, _, Iz = util.multioutput.build_XY(Z_list)
+
+        super(SparseGPCoregionalizedRegression, self).__init__(X, Y, Z, kernel, likelihood, inference_method=VarDTC(), Y_metadata={'output_index':self.output_index})
+        self['.*inducing'][:,-1].fix()

GPy/plotting/matplot_dep/models_plots.py

@@ -7,6 +7,7 @@ import Tango
 from base_plots import gpplot, x_frame1D, x_frame2D
 from ...util.misc import param_to_array
 from ...models.gp_coregionalized_regression import GPCoregionalizedRegression
+from ...models.sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression
 
 
 def plot_fit(model, plot_limits=None, which_data_rows='all',
@@ -86,7 +87,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
             lower = m - 2*np.sqrt(v)
             upper = m + 2*np.sqrt(v)
         else:
-            meta = {'output_index': Xgrid[:,-1:].astype(np.int)} if isinstance(model,GPCoregionalizedRegression) else None
+            if isinstance(model,GPCoregionalizedRegression) or isinstance(model,SparseGPCoregionalizedRegression):
+                meta = {'output_index': Xgrid[:,-1:].astype(np.int)}
+            else:
+                meta = None
             m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta)
             lower, upper = model.predict_quantiles(Xgrid, Y_metadata=meta)