some code to debug the sprase GP gradients with.

GPy/core/model.py

@@ -363,7 +363,8 @@ class model(parameterised):
                     grad_string = "{0:^{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name,r,d,g, ng, c0 = cols[0]+9, c1 = cols[1], c2 = cols[2], c3 = cols[3], c4 = cols[4])
                     print grad_string
 
-            print ''
-            
+            if verbose:
+                print ''
+
             return False
         return True

GPy/models/__init__.py

@@ -3,7 +3,7 @@
 
 
 from GP_regression import GP_regression
-from sparse_GP_regression import sparse_GP_regression
+from sparse_GP_regression import sparse_GP_regression, sgp_debugB, sgp_debugC, sgp_debugE
 from GPLVM import GPLVM
 from warped_GP import warpedGP
 from GP_EP import GP_EP

GPy/models/sparse_GP_regression.py

@@ -197,3 +197,127 @@ class sparse_GP_regression(GP_regression):
                 pb.errorbar(self.X[:,0], pb.ylim()[0]+np.zeros(self.N), xerr=2*np.sqrt(self.X_uncertainty.flatten()))
         if self.Q==2:
             pb.plot(self.Z[:,0],self.Z[:,1],'wo')
+
+class sgp_debugB(sparse_GP_regression):
+    def _computations(self):
+        self.V = self.beta*self.Y
+        self.psi1V = np.dot(self.psi1, self.V)
+        self.psi1VVpsi1 = np.dot(self.psi1V, self.psi1V.T)
+        self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
+        self.A = mdot(self.Lmi, self.beta*self.psi2, self.Lmi.T)
+        self.B = np.eye(self.M) + self.A
+        self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B)
+        self.LLambdai = np.dot(self.LBi, self.Lmi)
+        self.trace_K = self.psi0 - np.trace(self.A)/self.beta
+        self.LBL_inv = mdot(self.Lmi.T, self.Bi, self.Lmi)
+        self.C = mdot(self.LLambdai, self.psi1V)
+        self.G =  mdot(self.LBL_inv, self.psi1VVpsi1, self.LBL_inv.T)
+
+        # Compute dL_dpsi
+        self.dL_dpsi0 = - 0.5 * self.D * self.beta * np.ones(self.N)
+        self.dL_dpsi1 = np.zeros_like(self.psi1)
+        self.dL_dpsi2 = - 0.5 * self.beta * (self.D*( - self.Kmmi))
+
+        # Compute dL_dKmm
+        self.dL_dKmm = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi) # dB
+
+    def log_likelihood(self):
+        A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.beta))
+        B = -0.5*self.beta*self.D*self.trace_K
+        C = -0.5*self.D * self.B_logdet
+        D = -0.5*self.beta*self.trYYT
+        E = +0.5*np.sum(self.psi1VVpsi1 * self.LBL_inv)
+        return B
+
+    def dL_dbeta(self):
+        dA_dbeta =   0.5 * self.N*self.D/self.beta
+        dB_dbeta = - 0.5 * self.D * self.trace_K
+        dC_dbeta = - 0.5 * self.D * np.sum(self.Bi*self.A)/self.beta
+        dD_dbeta = - 0.5 * self.trYYT
+        tmp = mdot(self.LBi.T, self.LLambdai, self.psi1V)
+        dE_dbeta = (np.sum(np.square(self.C)) - 0.5 * np.sum(self.A * np.dot(tmp, tmp.T)))/self.beta
+        return np.squeeze(dB_dbeta)
+
+
+class sgp_debugC(sparse_GP_regression):
+    def _computations(self):
+        self.V = self.beta*self.Y
+        self.psi1V = np.dot(self.psi1, self.V)
+        self.psi1VVpsi1 = np.dot(self.psi1V, self.psi1V.T)
+        self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
+        self.A = mdot(self.Lmi, self.beta*self.psi2, self.Lmi.T)
+        self.B = np.eye(self.M) + self.A
+        self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B)
+        self.LLambdai = np.dot(self.LBi, self.Lmi)
+        self.trace_K = self.psi0 - np.trace(self.A)/self.beta
+        self.LBL_inv = mdot(self.Lmi.T, self.Bi, self.Lmi)
+        self.C = mdot(self.LLambdai, self.psi1V)
+        self.G =  mdot(self.LBL_inv, self.psi1VVpsi1, self.LBL_inv.T)
+
+        # Compute dL_dpsi
+        self.dL_dpsi0 = np.zeros(self.N)
+        self.dL_dpsi1 = np.zeros_like(self.psi1)
+        self.dL_dpsi2 = - 0.5 * self.beta * (self.D*(self.LBL_inv))
+
+        # Compute dL_dKmm
+        self.dL_dKmm = -0.5 * self.D * (- self.LBL_inv - 2.*self.beta*mdot(self.LBL_inv, self.psi2, self.Kmmi) + self.Kmmi) # dC
+
+    def log_likelihood(self):
+        A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.beta))
+        B = -0.5*self.beta*self.D*self.trace_K
+        C = -0.5*self.D * self.B_logdet
+        D = -0.5*self.beta*self.trYYT
+        E = +0.5*np.sum(self.psi1VVpsi1 * self.LBL_inv)
+        return C
+
+    def dL_dbeta(self):
+        dA_dbeta =   0.5 * self.N*self.D/self.beta
+        dB_dbeta = - 0.5 * self.D * self.trace_K
+        dC_dbeta = - 0.5 * self.D * np.sum(self.Bi*self.A)/self.beta
+        dD_dbeta = - 0.5 * self.trYYT
+        tmp = mdot(self.LBi.T, self.LLambdai, self.psi1V)
+        dE_dbeta = (np.sum(np.square(self.C)) - 0.5 * np.sum(self.A * np.dot(tmp, tmp.T)))/self.beta
+        return np.squeeze(dC_dbeta)
+
+
+class sgp_debugE(sparse_GP_regression):
+    def _computations(self):
+        self.V = self.beta*self.Y
+        self.psi1V = np.dot(self.psi1, self.V)
+        self.psi1VVpsi1 = np.dot(self.psi1V, self.psi1V.T)
+        self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
+        self.A = mdot(self.Lmi, self.beta*self.psi2, self.Lmi.T)
+        self.B = np.eye(self.M) + self.A
+        self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B)
+        self.LLambdai = np.dot(self.LBi, self.Lmi)
+        self.trace_K = self.psi0 - np.trace(self.A)/self.beta
+        self.LBL_inv = mdot(self.Lmi.T, self.Bi, self.Lmi)
+        self.C = mdot(self.LLambdai, self.psi1V)
+        self.G =  mdot(self.LBL_inv, self.psi1VVpsi1, self.LBL_inv.T)
+
+        # Compute dL_dpsi
+        self.dL_dpsi0 = np.zeros(self.N)
+        self.dL_dpsi1 = np.zeros_like(self.psi1)
+        self.dL_dpsi2 = - 0.5 * self.beta * (self.G)
+
+        # Compute dL_dKmm
+        self.dL_dKmm =  np.dot(np.dot(self.G,self.beta*self.psi2) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE
+
+    def log_likelihood(self):
+        A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.beta))
+        B = -0.5*self.beta*self.D*self.trace_K
+        C = -0.5*self.D * self.B_logdet
+        D = -0.5*self.beta*self.trYYT
+        E = +0.5*np.sum(self.psi1VVpsi1 * self.LBL_inv)
+        return E
+
+    def dL_dbeta(self):
+        dA_dbeta =   0.5 * self.N*self.D/self.beta
+        dB_dbeta = - 0.5 * self.D * self.trace_K
+        dC_dbeta = - 0.5 * self.D * np.sum(self.Bi*self.A)/self.beta
+        dD_dbeta = - 0.5 * self.trYYT
+        tmp = mdot(self.LBi.T, self.LLambdai, self.psi1V)
+        dE_dbeta = (np.sum(np.square(self.C)) - 0.5 * np.sum(self.A * np.dot(tmp, tmp.T)))/self.beta
+        return np.squeeze(dE_dbeta)
+
+