From a342dc2f2329e63d6beb4f9dce2f68f2f0fd9b2d Mon Sep 17 00:00:00 2001 From: Nicolas Date: Tue, 12 Mar 2013 16:43:51 +0000 Subject: [PATCH 1/3] errors handled in Mat32 --- GPy/kern/periodic_Matern32.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py index 898dff7b..662c1506 100644 --- a/GPy/kern/periodic_Matern32.py +++ b/GPy/kern/periodic_Matern32.py @@ -1,6 +1,7 @@ from kernpart import kernpart import numpy as np from GPy.util.linalg import mdot, pdinv +from GPy.util.decorators import silence_errors class periodic_Matern32(kernpart): """ @@ -39,12 +40,16 @@ class periodic_Matern32(kernpart): def f(x): return alpha*np.cos(omega*x+phase) return f + + @silence_errors def _cos_factorization(self,alpha,omega,phase): r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None] r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None] r = np.sqrt(r1**2 + r2**2) psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2)) return r,omega[:,0:1], psi + + @silence_errors def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2): Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) ) Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + np.cos(phi1-phi2.T)*(self.upper-self.lower) @@ -55,6 +60,7 @@ class periodic_Matern32(kernpart): def _get_params(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscale,self.period)) + def _set_params(self,x): """set the value of the parameters.""" assert x.size==3 @@ -101,6 +107,7 @@ class periodic_Matern32(kernpart): FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target) + @silence_errors def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X @@ -172,6 +179,7 @@ class periodic_Matern32(kernpart): #np.add(target[:,:,2],dK_dper, target[:,:,2]) target[2] += np.sum(dK_dper*dL_dK) + @silence_errors def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal covariance matrix with respect to the parameters""" FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) From 7a9b6ad1131fd2ae5bf62104f6f2b209b82ee121 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Tue, 12 Mar 2013 16:50:12 +0000 Subject: [PATCH 2/3] The warnings are now handeled properly in the periodic kernels --- GPy/kern/periodic_Matern32.py | 4 ++++ GPy/kern/periodic_Matern52.py | 9 +++++++++ GPy/kern/periodic_exponential.py | 10 +++++++++- 3 files changed, 22 insertions(+), 1 deletion(-) diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py index 662c1506..95684a02 100644 --- a/GPy/kern/periodic_Matern32.py +++ b/GPy/kern/periodic_Matern32.py @@ -1,3 +1,7 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + from kernpart import kernpart import numpy as np from GPy.util.linalg import mdot, pdinv diff --git a/GPy/kern/periodic_Matern52.py b/GPy/kern/periodic_Matern52.py index c533961f..07cb11ea 100644 --- a/GPy/kern/periodic_Matern52.py +++ b/GPy/kern/periodic_Matern52.py @@ -1,6 +1,11 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + from kernpart import kernpart import numpy as np from GPy.util.linalg import mdot, pdinv +from GPy.util.decorators import silence_errors class periodic_Matern52(kernpart): """ @@ -40,6 +45,7 @@ class periodic_Matern52(kernpart): return alpha*np.cos(omega*x+phase) return f + @silence_errors def _cos_factorization(self,alpha,omega,phase): r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None] r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None] @@ -57,6 +63,7 @@ class periodic_Matern52(kernpart): def _get_params(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscale,self.period)) + def _set_params(self,x): """set the value of the parameters.""" assert x.size==3 @@ -105,6 +112,7 @@ class periodic_Matern52(kernpart): FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target) + @silence_errors def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X @@ -184,6 +192,7 @@ class periodic_Matern52(kernpart): #np.add(target[:,:,2],dK_dper, target[:,:,2]) target[2] += np.sum(dK_dper*dL_dK) + @silence_errors def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters""" FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) diff --git a/GPy/kern/periodic_exponential.py b/GPy/kern/periodic_exponential.py index b966bbef..0018a8f9 100644 --- a/GPy/kern/periodic_exponential.py +++ b/GPy/kern/periodic_exponential.py @@ -1,6 +1,11 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + from kernpart import kernpart import numpy as np from GPy.util.linalg import mdot, pdinv +from GPy.util.decorators import silence_errors class periodic_exponential(kernpart): """ @@ -40,6 +45,7 @@ class periodic_exponential(kernpart): return alpha*np.cos(omega*x+phase) return f + @silence_errors def _cos_factorization(self,alpha,omega,phase): r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None] r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None] @@ -57,6 +63,7 @@ class periodic_exponential(kernpart): def _get_params(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscale,self.period)) + def _set_params(self,x): """set the value of the parameters.""" assert x.size==3 @@ -101,6 +108,7 @@ class periodic_exponential(kernpart): FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target) + @silence_errors def dK_dtheta(self,dL_dK,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X @@ -166,6 +174,7 @@ class periodic_exponential(kernpart): target[1] += np.sum(dK_dlen*dL_dK) target[2] += np.sum(dK_dper*dL_dK) + @silence_errors def dKdiag_dtheta(self,dL_dKdiag,X,target): """derivative of the diagonal of the covariance matrix with respect to the parameters""" FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) @@ -225,4 +234,3 @@ class periodic_exponential(kernpart): target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag) target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag) target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag) - From 525ef56dcac62a99109cf7217472130af1f58521 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 12 Mar 2013 17:04:48 +0000 Subject: [PATCH 3/3] increased stability of _compuations in sparse_GP --- GPy/models/sparse_GP.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index a348c9f4..3d44ad6b 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -153,14 +153,26 @@ class sparse_GP(GP): #self.partial_for_likelihood += -np.diag(np.dot((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) , self.psi1VVpsi1 ))*self.likelihood.precision #dD else: #likelihood is not heterscedatic - beta = self.likelihood.precision - dbeta = 0.5 * self.N*self.D/beta - 0.5 * np.sum(np.square(self.likelihood.Y)) - dbeta += - 0.5 * self.D * (self.psi0.sum() - np.trace(self.A)/beta*sf2) - dbeta += - 0.5 * self.D * trace_dot(self.Bi,self.A)/beta - dbeta += np.trace(self.Cpsi1VVpsi1)/beta - 0.5 * trace_dot(np.dot(self.C,self.psi2_beta_scaled) , self.Cpsi1VVpsi1 )/beta - self.partial_for_likelihood = -dbeta*self.likelihood.precision**2 + self.partial_for_likelihood = - 0.5 * self.N*self.D*self.likelihood.precision + 0.5 * np.sum(np.square(self.likelihood.Y))*self.likelihood.precision**2 + self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum()*self.likelihood.precision**2 - np.trace(self.A)*self.likelihood.precision*sf2) + self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self.likelihood.precision + self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1)) + + def log_likelihood(self): + """ Compute the (lower bound on the) log marginal likelihood """ + sf2 = self.scale_factor**2 + if self.likelihood.is_heteroscedastic: + A = -0.5*self.N*self.D*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y) + B = -0.5*self.D*(np.sum(self.likelihood.precision.flatten()*self.psi0) - np.trace(self.A)*sf2) + else: + A = -0.5*self.N*self.D*(np.log(2.*np.pi) + np.log(self.likelihood._variance)) -0.5*self.likelihood.precision*self.likelihood.trYYT + B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2) + C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2)) + D = 0.5*np.trace(self.Cpsi1VVpsi1) + return A+B+C+D + def _set_params(self, p): self.Z = p[:self.M*self.Q].reshape(self.M, self.Q) self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam]) @@ -188,18 +200,6 @@ class sparse_GP(GP): #self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0) self._set_params(self._get_params()) # update the GP - def log_likelihood(self): - """ Compute the (lower bound on the) log marginal likelihood """ - sf2 = self.scale_factor**2 - if self.likelihood.is_heteroscedastic: - A = -0.5*self.N*self.D*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y) - B = -0.5*self.D*(np.sum(self.likelihood.precision.flatten()*self.psi0) - np.trace(self.A)*sf2) - else: - A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.likelihood.precision)) -0.5*self.likelihood.precision*self.likelihood.trYYT - B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2) - C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2)) - D = 0.5*np.trace(self.Cpsi1VVpsi1) - return A+B+C+D def _log_likelihood_gradients(self): return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood)))