First attempt at making coregionalise work with the sparse model

Gradients are failing! have implemented prod_othogonal.dKdiag_dtheta

GPy/examples/regression.py

@@ -143,6 +143,43 @@ def coregionalisation_toy():
     return m
 
 
+def coregionalisation_sparse():
+    """
+    A simple demonstration of coregionalisation on two sinusoidal functions
+    """
+    X1 = np.random.rand(500,1)*8
+    X2 = np.random.rand(300,1)*5
+    index = np.vstack((np.zeros_like(X1),np.ones_like(X2)))
+    X = np.hstack((np.vstack((X1,X2)),index))
+    Y1 = np.sin(X1) + np.random.randn(*X1.shape)*0.05
+    Y2 = -np.sin(X2) + np.random.randn(*X2.shape)*0.05
+    Y = np.vstack((Y1,Y2))
+
+    Z = np.hstack((np.random.rand(25,1)*8,np.random.randint(0,2,25)[:,None]))
+
+    k1 = GPy.kern.rbf(1)
+    k2 = GPy.kern.coregionalise(2,2)
+    k = k1.prod_orthogonal(k2) + GPy.kern.white(2,0.001)
+
+    m = GPy.models.sparse_GP_regression(X,Y,kernel=k,Z=Z)
+    m.constrain_fixed('rbf_var',1.)
+    m.constrain_positive('kappa')
+    m.constrain_fixed('iip')
+    m.ensure_default_constraints()
+    #m.optimize()
+
+    pb.figure()
+    Xtest1 = np.hstack((np.linspace(0,9,100)[:,None],np.zeros((100,1))))
+    Xtest2 = np.hstack((np.linspace(0,9,100)[:,None],np.ones((100,1))))
+    mean, var,low,up = m.predict(Xtest1)
+    GPy.util.plot.gpplot(Xtest1[:,0],mean,low,up)
+    mean, var,low,up = m.predict(Xtest2)
+    GPy.util.plot.gpplot(Xtest2[:,0],mean,low,up)
+    pb.plot(X1[:,0],Y1[:,0],'rx',mew=2)
+    pb.plot(X2[:,0],Y2[:,0],'gx',mew=2)
+    return m
+
+
 
 def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000):
     """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""

GPy/kern/coregionalise.py

@@ -46,8 +46,8 @@ class coregionalise(kernpart):
             index2 = index
         else:
             index2 = np.asarray(index2,dtype=np.int)
-        ii,jj = np.meshgrid(index,index2)
+        ii,jj = np.meshgrid(index2,index)
-        target += self.B[ii,jj].T
+        target += self.B[ii,jj]
 
     def Kdiag(self,index,target):
         target += np.diag(self.B)[np.asarray(index,dtype=np.int).flatten()]
@@ -58,26 +58,47 @@ class coregionalise(kernpart):
             index2 = index
         else:
             index2 = np.asarray(index2,dtype=np.int)
-        ii,jj = np.meshgrid(index,index2)
+        ii,jj = np.meshgrid(index2,index)
         PK = np.zeros((self.R,self.R))
-        dkappa = np.zeros(self.Nout)
         partial_small = np.zeros_like(self.B)
         for i in range(self.Nout):
             for j in range(self.Nout):
-                partial_small[j,i] = np.sum(partial[(ii==i)*(jj==j)])
+                partial_small[i,j] = np.sum(partial[(ii==i)*(jj==j)])
-        #print partial_small
         dkappa = np.diag(partial_small)
 
-        ##target += (((X2[:, None, :] * self.variances)) * partial[:,:, None]).sum(0)
         dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0)
 
         target += np.hstack([dW.flatten(),dkappa])
 
     def dKdiag_dtheta(self,partial,index,target):
-        raise NotImplementedError
+        index = np.asarray(index,dtype=np.int).flatten()
+        partial_small = np.zeros(self.Nout)
+        for i in range(self.Nout):
+            partial_small[i] += np.sum(partial[index==i])
+        dW = 2.*self.W*partial_small[:,None]
+        dkappa = partial_small
+        target += np.hstack([dW.flatten(),dkappa])
 
     def dK_dX(self,partial,X,X2,target):
         pass
 
+    def dKdiag_dthetai_(self,partial,index,target):
+        index = np.asarray(index,dtype=np.int)
+        index2 = index
+        ii,jj = np.meshgrid(index2,index)
+        PK = np.zeros((self.R,self.R))
+        partial_small = np.zeros_like(self.B)
+        for i in range(self.Nout):
+            for j in range(self.Nout):
+                partial_small[j,i] = np.sum(partial[np.diag((ii==i)*(jj==j))])
+        #print partial_small
+        dkappa = np.diag(partial_small)
+
+        ##target += (((X2[:, None, :] * self.variances)) * partial[:,:, None]).sum(0)
+        partial_small = np.diag(np.diag(partial_small))
+        #dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0)
+        dW = 2.
+
+        target += np.hstack([dW.flatten(),dkappa])
 
 

GPy/kern/product_orthogonal.py

@@ -46,14 +46,6 @@ class product_orthogonal(kernpart):
         self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2)
         target += target1 * target2
 
-    def Kdiag(self,X,target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        target1 = np.zeros((X.shape[0],))
-        target2 = np.zeros((X.shape[0],))
-        self.k1.Kdiag(X[:,0:self.k1.D],target1)
-        self.k2.Kdiag(X[:,self.k1.D:],target2)
-        target += target1 * target2
-
     def dK_dtheta(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters."""
         if X2 is None: X2 = X
@@ -70,6 +62,22 @@ class product_orthogonal(kernpart):
         target[:self.k1.Nparam] += k1_target
         target[self.k1.Nparam:] += k2_target
 
+    def Kdiag(self,X,target):
+        """Compute the diagonal of the covariance matrix associated to X."""
+        target1 = np.zeros((X.shape[0],))
+        target2 = np.zeros((X.shape[0],))
+        self.k1.Kdiag(X[:,:self.k1.D],target1)
+        self.k2.Kdiag(X[:,self.k1.D:],target2)
+        target += target1 * target2
+
+    def dKdiag_dtheta(self,partial,X,target):
+        K1 = np.zeros(X.shape[0])
+        K2 = np.zeros(X.shape[0])
+        self.k1.Kdiag(X[:,:self.k1.D],K1)
+        self.k2.Kdiag(X[:,self.k1.D:],K2)
+        self.k1.dKdiag_dtheta(partial*K2,X[:,:self.k1.D],target[:self.k1.Nparam])
+        self.k2.dKdiag_dtheta(partial*K1,X[:,self.k1.D:],target[self.k1.Nparam:])
+
     def dK_dX(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to X."""
         if X2 is None: X2 = X