From 3a0a192362cde1ae238665c229fa6b84e98faece Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 25 Feb 2014 15:48:35 +0000 Subject: [PATCH] more efficient computations in stationary --- GPy/kern/__init__.py | 31 -------------- GPy/kern/_src/stationary.py | 84 +++++++++++++++++++++++++++++++------ 2 files changed, 71 insertions(+), 44 deletions(-) diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index 826b15aa..030cba66 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -3,40 +3,9 @@ from _src.rbf import RBF from _src.linear import Linear from _src.static import Bias, White from _src.brownian import Brownian -from _src.stationary import Exponential, Matern32, Matern52, ExpQuad from _src.sympykern import Sympykern -#from _src.kern import kern_test -#import coregionalize -#import exponential -#import eq_ode1 -#import finite_dimensional -#import fixed -#import gibbs -#import hetero -#import hierarchical -#import independent_outputs -#import linear -#import Matern32 -#import Matern52 -#import mlp -#import ODE_1 -#import periodic_exponential -#import periodic_Matern32 -#import periodic_Matern52 -#import poly -#import prod_orthogonal -#import prod -#import rational_quadratic -#import rbfcos -#import rbf -#import rbf_inv -#import spline -#import symmetric -#import white -======= from _src.stationary import Exponential, Matern32, Matern52, ExpQuad, RatQuad, Cosine from _src.mlp import MLP from _src.periodic import PeriodicExponential, PeriodicMatern32, PeriodicMatern52 from _src.independent_outputs import IndependentOutputs, Hierarchical from _src.coregionalize import Coregionalize ->>>>>>> da4686dd3c8db8639b0c3c6e30609d0b3fa59130 diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index a2a83929..b998969c 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -5,6 +5,7 @@ from kern import Kern from ...core.parameterization import Param from ...core.parameterization.transformations import Logexp +from ...util.linalg import tdot from ... import util import numpy as np from scipy import integrate @@ -33,10 +34,10 @@ class Stationary(Kern): self.add_parameters(self.variance, self.lengthscale) def K_of_r(self, r): - raise NotImplementedError, "implement the covaraiance functino and a fn of r to use this class" + raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" def dK_dr(self, r): - raise NotImplementedError, "implement the covaraiance functino and a fn of r to use this class" + raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" def K(self, X, X2=None): r = self._scaled_dist(X, X2) @@ -47,8 +48,44 @@ class Stationary(Kern): X2 = X return X[:, None, :] - X2[None, :, :] + def _unscaled_dist(self, X, X2=None): + """ + Compute the square distance between each row of X and X2, or between + each pair of rows of X if X2 is None. + """ + if X2 is None: + Xsq = np.sum(np.square(X),1) + return np.sqrt(-2.*tdot(X) + (Xsq[:,None] + Xsq[None,:])) + else: + X1sq = np.sum(np.square(X),1) + X2sq = np.sum(np.square(X2),1) + return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) + def _scaled_dist(self, X, X2=None): - return np.sqrt(np.sum(np.square(self._dist(X, X2) / self.lengthscale), -1)) + """ + Efficiently compute the scaled distance, r. + + r = \sum_{q=1}^Q (x_q - x'q)^2/l_q^2 + + Note that if thre is only one lengthscale, l comes outside the sum. In + this case we compute the unscaled distance first (in a separate + function for caching) and divide by lengthscale afterwards + + """ + if self.ARD: + if X2 is None: + Xl = X/self.lengthscale + Xsq = np.sum(np.square(Xl),1) + return np.sqrt(np.sqrt(-2.*tdot(Xl) +(Xsq[:,None] + Xsq[None,:]))) + else: + X1l = X/self.lengthscale + X2l = X2/self.lengthscale + X1sq = np.sum(np.square(X1l),1) + X2sq = np.sum(np.square(X2l),1) + return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) + else: + return self._unscaled_dist(X, X2)/self.lengthscale + def Kdiag(self, X): ret = np.empty(X.shape[0]) @@ -65,16 +102,22 @@ class Stationary(Kern): rinv = self._inv_dist(X, X2) dL_dr = self.dK_dr(r) * dL_dK - x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 if self.ARD: + x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0) else: + x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum() self.variance.gradient = np.sum(K * dL_dK)/self.variance def _inv_dist(self, X, X2=None): + """ + Compute the elementwise inverse of the distance matrix, expecpt on the + diagonal, where we return zero (the distance on the diagonal is zero). + This term appears in derviatives. + """ dist = self._scaled_dist(X, X2) if X2 is None: nondiag = util.diag.offdiag_view(dist) @@ -84,12 +127,25 @@ class Stationary(Kern): return 1./np.where(dist != 0., dist, np.inf) def gradients_X(self, dL_dK, X, X2=None): + """ + Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X + """ r = self._scaled_dist(X, X2) - dL_dr = self.dK_dr(r) * dL_dK invdist = self._inv_dist(X, X2) - ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2 + dL_dr = self.dK_dr(r) * dL_dK + #The high-memory numpy way: ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2 + #if X2 is None: + #ret *= 2. + + #the lower memory way with a loop + tmp = invdist*dL_dr if X2 is None: - ret *= 2. + tmp *= 2. + X2 = X + ret = np.empty(X.shape, dtype=np.float64) + [np.copyto(ret[:,q], np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), 1)) for q in xrange(self.input_dim)] + ret /= self.lengthscale**2 + return ret def gradients_X_diag(self, dL_dKdiag, X): @@ -241,17 +297,19 @@ class RatQuad(Stationary): self.add_parameters(self.power) def K_of_r(self, r): - return self.variance*(1. + r**2/2.)**(-self.power) + r2 = np.power(r, 2.) + return self.variance*np.power(1. + r2/2., -self.power) def dK_dr(self, r): - return -self.variance*self.power*r*(1. + r**2/2)**(-self.power - 1.) + r2 = np.power(r, 2.) + return -self.variance*self.power*r*np.power(1. + r2/2., - self.power - 1.) def update_gradients_full(self, dL_dK, X, X2=None): super(RatQuad, self).update_gradients_full(dL_dK, X, X2) r = self._scaled_dist(X, X2) - r2 = r**2 - dpow = -2.**self.power*(r2 + 2.)**(-self.power)*np.log(0.5*(r2+2.)) - self.power.gradient = np.sum(dL_dK*dpow) - + r2 = np.power(r, 2.) + dK_dpow = -self.variance * np.power(2., self.power) * np.power(r2 + 2., -self.power) * np.log(0.5*(r2+2.)) + grad = np.sum(dL_dK*dK_dpow) + self.power.gradient = grad