Checkgrads with explicit and implicit components half the time

GPy/examples/laplace_approximations.py

@@ -39,11 +39,11 @@ def debug_student_t_noise_approx():
     plot = False
     real_var = 0.1
     #Start a function, any function
-    #X = np.linspace(0.0, 10.0, 100)[:, None]
+    X = np.linspace(0.0, 10.0, 15)[:, None]
-    X = np.array([0.5])[:, None]
+    #X = np.array([0.5])[:, None]
     Y = np.sin(X) + np.random.randn(*X.shape)*real_var
 
-    X_full = np.linspace(0.0, 10.0, 500)[:, None]
+    X_full = np.linspace(0.0, 10.0, 15)[:, None]
     Y_full = np.sin(X_full)
 
     Y = Y/Y.max()
@@ -83,7 +83,8 @@ def debug_student_t_noise_approx():
         #plt.plot(X_full, Y_full)
     #print m
 
-    edited_real_sd = initial_var_guess #real_sd
+    #edited_real_sd = initial_var_guess #real_sd
+    edited_real_sd = real_sd
 
     print "Clean student t, rasm"
     t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
@@ -94,7 +95,7 @@ def debug_student_t_noise_approx():
     #m.constrain_fixed('rbf_l', 1.8651)
     #m.constrain_fixed('t_noise_variance', real_sd)
     m.constrain_positive('rbf')
-    m.constrain_fixed('t_noi', real_sd)
+    #m.constrain_fixed('t_noi', real_sd)
     m.ensure_default_constraints()
     m.update_likelihood_approximation()
     m.optimize(messages=True)
@@ -148,7 +149,7 @@ def student_t_approx():
     #Yc = Yc/Yc.max()
 
     #Add student t random noise to datapoints
-    deg_free = 10
+    deg_free = 8
     real_sd = np.sqrt(real_var)
     print "Real noise: ", real_sd
 
@@ -202,8 +203,6 @@ def student_t_approx():
     plt.title('Gaussian corrupt')
     print m
 
-    import ipdb; ipdb.set_trace() ### XXX BREAKPOINT
-
     plt.figure(2)
     plt.suptitle('Student-t likelihood')
     edited_real_sd = real_sd #initial_var_guess
@@ -236,33 +235,35 @@ def student_t_approx():
     plt.ylim(-2.5, 2.5)
     plt.title('Student-t rasm corrupt')
 
-    print "Clean student t, ncg"
+    return m
-    t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
-    stu_t_likelihood = GPy.likelihoods.Laplace(Y, t_distribution, rasm=False)
-    m = GPy.models.GP(X, stu_t_likelihood, kernel3)
-    m.ensure_default_constraints()
-    m.update_likelihood_approximation()
-    m.optimize()
-    print(m)
-    plt.subplot(221)
-    m.plot()
-    plt.plot(X_full, Y_full)
-    plt.ylim(-2.5, 2.5)
-    plt.title('Student-t ncg clean')
 
-    print "Corrupt student t, ncg"
+    #print "Clean student t, ncg"
-    t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
+    #t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
-    corrupt_stu_t_likelihood = GPy.likelihoods.Laplace(Yc.copy(), t_distribution, rasm=False)
+    #stu_t_likelihood = GPy.likelihoods.Laplace(Y, t_distribution, rasm=False)
-    m = GPy.models.GP(X, corrupt_stu_t_likelihood, kernel5)
+    #m = GPy.models.GP(X, stu_t_likelihood, kernel3)
-    m.ensure_default_constraints()
+    #m.ensure_default_constraints()
-    m.update_likelihood_approximation()
+    #m.update_likelihood_approximation()
-    m.optimize()
+    #m.optimize()
-    print(m)
+    #print(m)
-    plt.subplot(223)
+    #plt.subplot(221)
-    m.plot()
+    #m.plot()
-    plt.plot(X_full, Y_full)
+    #plt.plot(X_full, Y_full)
-    plt.ylim(-2.5, 2.5)
+    #plt.ylim(-2.5, 2.5)
-    plt.title('Student-t ncg corrupt')
+    #plt.title('Student-t ncg clean')
+
+    #print "Corrupt student t, ncg"
+    #t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
+    #corrupt_stu_t_likelihood = GPy.likelihoods.Laplace(Yc.copy(), t_distribution, rasm=False)
+    #m = GPy.models.GP(X, corrupt_stu_t_likelihood, kernel5)
+    #m.ensure_default_constraints()
+    #m.update_likelihood_approximation()
+    #m.optimize()
+    #print(m)
+    #plt.subplot(223)
+    #m.plot()
+    #plt.plot(X_full, Y_full)
+    #plt.ylim(-2.5, 2.5)
+    #plt.title('Student-t ncg corrupt')
 
 
     ###with a student t distribution, since it has heavy tails it should work well

GPy/likelihoods/Laplace.py

@@ -8,9 +8,6 @@ from ..util.linalg import pdinv, mdot, jitchol, chol_inv, det_ln_diag, pddet
 from scipy.linalg.lapack import dtrtrs
 import random
 #import pylab as plt
-np.random.seed(50)
-random.seed(50)
-
 
 class Laplace(likelihood):
     """Laplace approximation to a posterior"""
@@ -45,7 +42,7 @@ class Laplace(likelihood):
         self.is_heteroscedastic = True
         self.Nparams = 0
 
-        self.NORMAL_CONST = -((0.5 * self.N) * np.log(2 * np.pi))
+        self.NORMAL_CONST = ((0.5 * self.N) * np.log(2 * np.pi))
 
         #Initial values for the GP variables
         self.Y = np.zeros((self.N, 1))
@@ -72,26 +69,36 @@ class Laplace(likelihood):
         #FIXME: Careful of side effects! And make sure W and K are up to date!
         d3lik_d3fhat = self.likelihood_function.d3lik_d3f(self.data, self.f_hat)
         dL_dfhat = -0.5*(np.diag(self.Ki_W_i)[:, None]*d3lik_d3fhat)
+
         Wi_K_i = self.W_12*self.Bi*self.W_12.T #same as rasms R
 
         I_KW_i = np.eye(self.N) - np.dot(self.K, Wi_K_i)
         return dL_dfhat, I_KW_i, Wi_K_i
 
-    def _Kgradients(self, dK_dthetaK):
+    def _Kgradients(self, dK_dthetaK, X):
         """
         Gradients with respect to prior kernel parameters
         """
         dL_dfhat, I_KW_i, Wi_K_i = self._shared_gradients_components()
         dlp = self.likelihood_function.dlik_df(self.data, self.f_hat)
 
-        dL_dthetaK = np.zeros(dK_dthetaK.shape)
-        for thetaK_i, dK_dthetaK_i in enumerate(dK_dthetaK):
-            #Explicit
-            f_Ki_dK_dtheta_Ki_f = mdot(self.Ki_f.T, dK_dthetaK_i, self.Ki_f)
-            dL_dthetaK[thetaK_i] = 0.5*f_Ki_dK_dtheta_Ki_f - 0.5*np.trace(Wi_K_i*dK_dthetaK_i)
         #Implicit
-            df_hat_dthetaK = mdot(I_KW_i, dK_dthetaK_i, dlp)
+        impl = mdot(dlp, dL_dfhat.T, I_KW_i)
-            dL_dthetaK[thetaK_i] += np.dot(dL_dfhat.T, df_hat_dthetaK)
+        expl_a = - mdot(self.Ki_f, self.Ki_f.T)
+        expl_b = Wi_K_i
+        expl = 0.5*expl_a - 0.5*expl_b
+        dL_dthetaK_exp = dK_dthetaK(expl, X)
+        dL_dthetaK_imp = dK_dthetaK(impl, X)
+        dL_dthetaK = -(dL_dthetaK_imp + dL_dthetaK_exp)
+
+        #dL_dthetaK = np.zeros(dK_dthetaK.shape)
+        #for thetaK_i, dK_dthetaK_i in enumerate(dK_dthetaK):
+            ##Explicit
+            #f_Ki_dK_dtheta_Ki_f = mdot(self.Ki_f.T, dK_dthetaK_i, self.Ki_f)
+            #dL_dthetaK[thetaK_i] = 0.5*f_Ki_dK_dtheta_Ki_f - 0.5*np.trace(Wi_K_i*dK_dthetaK_i)
+            ##Implicit
+            #df_hat_dthetaK = mdot(I_KW_i, dK_dthetaK_i, dlp)
+            #dL_dthetaK[thetaK_i] += np.dot(dL_dfhat.T, df_hat_dthetaK)
 
         return dL_dthetaK
 
@@ -99,13 +106,12 @@ class Laplace(likelihood):
         """
         Gradients with respect to likelihood parameters
         """
-        return np.zeros(1)
+        #return np.zeros(1)
-        #return np.zeros(0)
         dL_dfhat, I_KW_i, Wi_K_i = self._shared_gradients_components()
         dlik_dthetaL, dlik_grad_dthetaL, dlik_hess_dthetaL = self.likelihood_function._gradients(self.data, self.f_hat)
 
         num_params = len(dlik_dthetaL)
-        dL_dthetaL = np.zeros((1, num_params)) # make space for one derivative for each likelihood parameter
+        dL_dthetaL = np.zeros(num_params) # make space for one derivative for each likelihood parameter
         for thetaL_i in range(num_params):
             #Explicit
             #dL_dthetaL[thetaL_i] = np.sum(dlik_dthetaL[thetaL_i]) - 0.5*np.trace(np.dot(Ki_W_i.T, np.diagflat(dlik_hess_dthetaL[thetaL_i])))
@@ -143,8 +149,6 @@ class Laplace(likelihood):
         $$\tilde{\Sigma} = W^{-1}$$
 
         """
-        epsilon = 1e14
-
         #Wi(Ki + W) = WiKi + I = KW_i + I = L_Lt_W_i + I = Wi_Lit_Li + I = Lt_W_i_Li + I
         #dtritri -> L -> L_i
         #dtrtrs -> L.T*W, L_i -> (L.T*W)_i*L_i
@@ -153,54 +157,38 @@ class Laplace(likelihood):
         Li = chol_inv(L)
         Lt_W = L.T*self.W.T
 
-        ##Check it isn't singular!
-        if cond(Lt_W) > epsilon:
-            print "WARNING: L_inv.T * W matrix is singular,\nnumerical stability may be a problem"
-
         Lt_W_i_Li = dtrtrs(Lt_W, Li, lower=False)[0]
         self.Wi__Ki_W = Lt_W_i_Li + np.eye(self.N)
 
         Y_tilde = np.dot(self.Wi__Ki_W, self.f_hat)
 
-        #f.T(Ki + W)f
+        ln_W_det = det_ln_diag(self.W)
-        f_Ki_W_f = (np.dot(self.f_hat.T, cho_solve((L, True), self.f_hat))
+        yf_W_yf = mdot((Y_tilde - self.f_hat).T, np.diagflat(self.W), (Y_tilde - self.f_hat))
-                    + mdot(self.f_hat.T, self.W*self.f_hat)
-                    )
 
-        y_W_f = mdot(Y_tilde.T*self.W.T, self.f_hat)
+        #Z_tilde = (+ self.NORMAL_CONST
-
-
-        y_W_y = mdot(Y_tilde.T, self.W*Y_tilde)
-
-        ln_W_det = np.log(self.W).sum()
-
-        #FIXME: Revisit this
-        Z_tilde = (- self.NORMAL_CONST
-                   + 0.5*self.ln_K_det
-                   + 0.5*ln_W_det
-                   + 0.5*self.ln_Ki_W_i_det
-                   + 0.5*f_Ki_W_f
-                   + 0.5*y_W_y
-                   - y_W_f
-                   + self.ln_z_hat
-                   )
-        #Z_tilde = (self.NORMAL_CONST
-                   #- 0.5*self.ln_K_det
-                   #- 0.5*ln_W_det
-                   #- 0.5*self.ln_Ki_W_i_det
-                   #- 0.5*f_Ki_W_f
-                   #- 0.5*y_W_y
-                   #+ y_W_f
                    #+ self.ln_z_hat
+                   #+ 0.5*self.ln_I_KW_det
+                   #- 0.5*ln_W_det
+                   #+ 0.5*self.f_Ki_f
+                   #+ 0.5*yf_W_yf
                    #)
-        #self.Z_tilde = 0
-
-        ##Check it isn't singular!
-        if cond(self.W) > epsilon:
-            print "WARNING: Transformed covariance matrix is singular,\nnumerical stability may be a problem"
 
         self.Sigma_tilde = np.diagflat(1.0/self.W)
 
+        Ki, _, _, K_det = pdinv(self.K)
+        ln_det_K_Wi__Bi = self.ln_I_KW_det + pddet(self.Sigma_tilde + self.K)
+        W = np.diagflat(self.W)
+        Wi = self.Sigma_tilde
+        W12i = np.sqrt(Wi)
+        D = Ki - mdot((Ki + W), W12i, self.Bi, W12i, (Ki + W))
+        fDf = mdot(self.f_hat.T, D, self.f_hat)
+        l = self.likelihood_function.link_function(self.data, self.f_hat, extra_data=self.extra_data)
+        Z_tilde = (+ self.NORMAL_CONST
+                   + l
+                   + 0.5*ln_det_K_Wi__Bi
+                   - 0.5*fDf
+                  )
+
         #Convert to float as its (1, 1) and Z must be a scalar
         self.Z = np.float64(Z_tilde)
         self.Y = Y_tilde
@@ -239,10 +227,6 @@ class Laplace(likelihood):
         self.B, self.B_chol, self.W_12 = self._compute_B_statistics(self.K, self.W)
         self.Bi, _, _, B_det = pdinv(self.B)
 
-        self.Ki_W_i = self.K - mdot(self.K, self.W_12*self.Bi*self.W_12.T, self.K)
-
-        self.ln_Ki_W_i_det = np.linalg.det(self.Ki_W_i)
-
         #Do the computation again at f to get Ki_f which is useful
         b = self.W*self.f_hat + self.likelihood_function.dlik_df(self.data, self.f_hat, extra_data=self.extra_data)
         solve_chol = cho_solve((self.B_chol, True), np.dot(self.W_12*self.K, b))
@@ -250,12 +234,14 @@ class Laplace(likelihood):
         self.Ki_f = a
 
         self.f_Ki_f = np.dot(self.f_hat.T, self.Ki_f)
-        self.ln_K_det = pddet(self.K)
+        self.Ki_W_i = self.K - mdot(self.K, self.W_12*self.Bi*self.W_12.T, self.K)
-        #_, _, _, self.ln_K_det = pdinv(self.K)
 
+        #For det, |I + KW| == |I + W_12*K*W_12|
+        self.ln_I_KW_det = pddet(np.eye(self.N) + self.W_12*self.K*self.W_12.T)
+
+        #self.ln_I_KW_det = pddet(np.eye(self.N) + np.dot(self.K, self.W))
         self.ln_z_hat = (- 0.5*self.f_Ki_f
-                         - 0.5*self.ln_K_det
+                         - self.ln_I_KW_det
-                         + 0.5*self.ln_Ki_W_i_det
                          + self.likelihood_function.link_function(self.data, self.f_hat, extra_data=self.extra_data)
                          )
 
@@ -289,7 +275,7 @@ class Laplace(likelihood):
         #ONLY WORKS FOR 1D DATA
         def obj(f):
             res = -1 * (self.likelihood_function.link_function(self.data[:, 0], f, extra_data=self.extra_data) - 0.5 * np.dot(f.T, np.dot(self.Ki, f))
-                        + self.NORMAL_CONST)
+                        - self.NORMAL_CONST)
             return float(res)
 
         def obj_grad(f):

GPy/models/GP.py

@@ -141,6 +141,8 @@ class GP(model):
 
         Note, we use the chain rule: dL_dtheta = dL_dK * d_K_dtheta
         """
+        self.likelihood.fit_full(self.kern.K(self.X))
+        self.likelihood._set_params(self.likelihood._get_params())
         dL_dthetaK = self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X)
         if isinstance(self.likelihood, Laplace):
             #Reapproximate incase it hasnt been done...
@@ -155,8 +157,9 @@ class GP(model):
             #BUG: THIS SHOULD NOT BE (1,num_k_params) matrix it should be (N,N,num_k_params)
             #dK_dthetaK = self.kern.dK_dtheta(dL_dK=fake_dL_dKs, X=self.X)
 
-            #dL_dthetaK = self.likelihood._Kgradients(dK_dthetaK=dK_dthetaK)
+            dK_dthetaK = self.kern.dK_dtheta
-            dL_dthetaL = 0 #self.likelihood._gradients(partial=np.diag(self.dL_dK))
+            dL_dthetaK = self.likelihood._Kgradients(dK_dthetaK, self.X)
+            dL_dthetaL = self.likelihood._gradients(partial=np.diag(self.dL_dK))
             #print "dL_dthetaK after: ",dL_dthetaK
             #print "Stacked dL_dthetaK, dL_dthetaL: ", np.hstack((dL_dthetaK, dL_dthetaL))
         else:

GPy/util/linalg.py

@@ -34,7 +34,7 @@ def det_ln_diag(A):
 
 def pddet(A):
     """
-    Determinant of a positive definite matrix
+    Determinant of a positive definite matrix, only symmetric matricies though
     """
     L = jitchol(A)
     logdetA = 2*sum(np.log(np.diag(L)))