From 6c2975079517364f00b2345f0ef9b3d2f5a14103 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 31 May 2013 16:59:54 +0100 Subject: [PATCH] Took out all the asserts and using pure broadcasting method of diagonal now --- GPy/examples/laplace_approximations.py | 4 +- GPy/likelihoods/Laplace.py | 70 ++++++-------------------- GPy/models/GP.py | 3 +- 3 files changed, 20 insertions(+), 57 deletions(-) diff --git a/GPy/examples/laplace_approximations.py b/GPy/examples/laplace_approximations.py index 5103eefb..14ff44a0 100644 --- a/GPy/examples/laplace_approximations.py +++ b/GPy/examples/laplace_approximations.py @@ -39,8 +39,8 @@ 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.array([0.5])[:, None] + #X = np.linspace(0.0, 10.0, 100)[:, 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] diff --git a/GPy/likelihoods/Laplace.py b/GPy/likelihoods/Laplace.py index af74755f..74d37d48 100644 --- a/GPy/likelihoods/Laplace.py +++ b/GPy/likelihoods/Laplace.py @@ -69,9 +69,7 @@ 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 = mdot(np.diagflat(self.W_12), self.Bi, np.diagflat(self.W_12)) #same as rasms R - Wi_K_inew = self.W_12*self.Bi*self.W_12.T #same as rasms R - assert np.all(Wi_K_i == Wi_K_inew) + 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 @@ -150,9 +148,7 @@ class Laplace(likelihood): #((L.T*w)_i + I)f_hat = y_tilde L = jitchol(self.K) Li = chol_inv(L) - Lt_W = np.dot(L.T, np.diagflat(self.W)) #FIXME: Can make Faster - Lt_Wnew = L.T*self.W.T - assert np.all(Lt_Wnew == Lt_W) + Lt_W = L.T*self.W.T ##Check it isn't singular! if cond(Lt_W) > epsilon: @@ -164,25 +160,15 @@ class Laplace(likelihood): #f.T(Ki + W)f f_Ki_W_f = (np.dot(self.f_hat.T, cho_solve((L, True), self.f_hat)) - + mdot(self.f_hat.T, np.diagflat(self.W), self.f_hat) - ) - f_Ki_W_fnew = (np.dot(self.f_hat.T, cho_solve((L, True), self.f_hat)) + mdot(self.f_hat.T, self.W*self.f_hat) ) - assert np.all(f_Ki_W_f == f_Ki_W_fnew) - y_W_f = mdot((Y_tilde.T, np.diagflat(self.W)), self.f_hat) - y_W_fnew = mdot(Y_tilde.T*self.W.T, self.f_hat) - assert np.all(y_W_f == y_W_fnew) + y_W_f = mdot(Y_tilde.T*self.W.T, self.f_hat) - y_W_y = mdot((Y_tilde.T, np.diagflat(self.W)), Y_tilde) - y_W_ynew = mdot(Y_tilde.T, self.W*Y_tilde) - assert np.all(y_W_y == y_W_ynew) + y_W_y = mdot(Y_tilde.T, self.W*Y_tilde) - ln_W_det = det_ln_diag(np.diagflat(self.W)) - ln_W_detnew = np.log(self.W).sum() - assert np.all(ln_W_det == ln_W_detnew) + ln_W_det = np.log(self.W).sum() #FIXME: Revisit this Z_tilde = (- self.NORMAL_CONST @@ -203,15 +189,13 @@ class Laplace(likelihood): #+ y_W_f #+ self.ln_z_hat #) - self.Z_tilde = 0 + #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 = inv(np.diagflat(self.W)) # Damn - Sigma_tildenew = np.diagflat(1.0/self.W) - assert np.all(self.Sigma_tilde == Sigma_tildenew) + self.Sigma_tilde = np.diagflat(1.0/self.W) #Convert to float as its (1, 1) and Z must be a scalar self.Z = np.float64(Z_tilde) @@ -251,23 +235,15 @@ 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, (np.diagflat(self.W_12), self.Bi, np.diagflat(self.W_12)), self.K) # Funky, order matters on stability! - Ki_W_inew = self.K - mdot(self.K, self.W_12*self.Bi*self.W_12.T, self.K) - assert np.all(self.Ki_W_i == Ki_W_inew) + 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) - b = np.dot(np.diagflat(self.W), self.f_hat) + self.likelihood_function.dlik_df(self.data, self.f_hat, extra_data=self.extra_data) - bnew = self.W*self.f_hat + self.likelihood_function.dlik_df(self.data, self.f_hat, extra_data=self.extra_data) - assert np.all(b == bnew) + 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), mdot((np.diagflat(self.W_12), self.K), b)) - solve_cholnew = cho_solve((self.B_chol, True), np.dot(self.W_12*self.K, b)) - assert np.all(solve_chol == solve_cholnew) + solve_chol = cho_solve((self.B_chol, True), np.dot(self.W_12*self.K, b)) - a = b - mdot(np.diagflat(self.W_12), solve_chol) - anew = b - self.W_12*solve_chol - assert np.all(a == anew) + a = b - self.W_12*solve_chol self.Ki_f = a self.f_Ki_f = np.dot(self.f_hat.T, self.Ki_f) @@ -291,10 +267,6 @@ class Laplace(likelihood): """ #W is diagonal so its sqrt is just the sqrt of the diagonal elements W_12 = np.sqrt(W) - # FIXME Take this out when you've done multiinput, Weirdly this is - # better when its W_12.T*K*W_12 which shouldnt make a difference - # because K is symmetrical - assert np.allclose(W_12*K*W_12.T, np.dot(np.diagflat(W_12), np.dot(K, np.diagflat(W_12)))) B = np.eye(self.N) + W_12*K*W_12.T L = jitchol(B) return (B, L, W_12) @@ -360,9 +332,7 @@ class Laplace(likelihood): # This is a property only held by non-log-concave likelihoods B, L, W_12 = self._compute_B_statistics(K, W) - W_f = np.dot(np.diagflat(W), f) - W_fnew = W*f - assert np.all(W_f == W_fnew) + W_f = W*f grad = self.likelihood_function.dlik_df(self.data, f, extra_data=self.extra_data) #Find K_i_f b = W_f + grad @@ -370,21 +340,13 @@ class Laplace(likelihood): #a should be equal to Ki*f now so should be able to use it c = np.dot(K, W_f) + f*(1-step_size) + step_size*np.dot(K, grad) - solve_L = cho_solve((L, True), np.dot(np.diagflat(W_12), c)) - solve_Lnew = cho_solve((L, True), W_12*c) - assert np.all(solve_L == solve_Lnew) + solve_L = cho_solve((L, True), W_12*c) - f = c - np.dot(K, np.dot(np.diagflat(W_12), solve_L)) - fnew = c - np.dot(K, W_12*solve_L) - assert np.all(f == fnew) + f = c - np.dot(K, W_12*solve_L) - solve_L = cho_solve((L, True), np.dot(np.diagflat(W_12), np.dot(K, b))) - solve_Lnew = cho_solve((L, True), W_12*np.dot(K, b)) - assert np.all(solve_L == solve_Lnew) + solve_L = cho_solve((L, True), W_12*np.dot(K, b)) - a = b - np.dot(np.diagflat(W_12), solve_L) - anew = b - W_12*solve_L - assert np.all(a == anew) + a = b - W_12*solve_L tmp_old_obj = old_obj old_obj = new_obj diff --git a/GPy/models/GP.py b/GPy/models/GP.py index 787429de..0ba20d7b 100644 --- a/GPy/models/GP.py +++ b/GPy/models/GP.py @@ -152,8 +152,9 @@ class GP(model): #Need to pass in a matrix of ones to get access to raw dK_dthetaK values without being chained fake_dL_dKs = np.ones(self.dL_dK.shape) #FIXME: Check this is right... #fake_dL_dKs = np.eye(self.dL_dK.shape[0]) #FIXME: Check this is right... + + #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) - #THIS SHOULD NOT BE (1,num_k_params) matrix it should be (N,N,num_k_params) dL_dthetaK = self.likelihood._Kgradients(dK_dthetaK=dK_dthetaK) dL_dthetaL = self.likelihood._gradients(partial=np.diag(self.dL_dK))