From 56ecd4782a576c4b471379b594aaa9639f90e799 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 22 Apr 2013 11:59:32 +0100 Subject: [PATCH] made the basic GP class use dtrtrs where possible --- GPy/models/GP.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/GPy/models/GP.py b/GPy/models/GP.py index cfda0cfe..a46a35d0 100644 --- a/GPy/models/GP.py +++ b/GPy/models/GP.py @@ -9,6 +9,7 @@ from ..core import model from ..util.linalg import pdinv,mdot from ..util.plot import gpplot,x_frame1D,x_frame2D, Tango from ..likelihoods import EP +from scipy import linalg class GP(model): """ @@ -78,10 +79,13 @@ class GP(model): #the gradient of the likelihood wrt the covariance matrix if self.likelihood.YYT is None: - alpha = np.dot(self.Ki,self.likelihood.Y) + #alpha = np.dot(self.Ki,self.likelihood.Y) + alpha,info = linalg.lapack.flapack.dpotrs(self.L,np.asfortranarray(self.likelihood.Y),lower=1) self.dL_dK = 0.5*(np.dot(alpha,alpha.T)-self.D*self.Ki) else: - tmp = mdot(self.Ki, self.likelihood.YYT, self.Ki) + #tmp = mdot(self.Ki, self.likelihood.YYT, self.Ki) + tmp,info = linalg.lapack.flapack.dpotrs(self.L,np.asfortranarray(self.likelihood.YYT),lower=1) + tmp,info = linalg.lapack.flapack.dpotrs(self.L,np.asfortranarray(tmp.T),lower=1) self.dL_dK = 0.5*(tmp - self.D*self.Ki) def _get_params(self): @@ -105,10 +109,13 @@ class GP(model): Computes the model fit using YYT if it's available """ if self.likelihood.YYT is None: - return -0.5*np.sum(np.square(np.dot(self.Li,self.likelihood.Y))) + #return -0.5*np.sum(np.square(np.dot(self.Li,self.likelihood.Y))) + tmp,info = linalg.lapack.flapack.dtrtrs(self.L,np.asfortranarray(self.likelihood.Y),lower=1) + return -0.5*np.sum(np.square(tmp)) else: return -0.5*np.sum(np.multiply(self.Ki, self.likelihood.YYT)) + def log_likelihood(self): """ The log marginal likelihood of the GP. @@ -136,8 +143,11 @@ class GP(model): for normalization or likelihood """ Kx = self.kern.K(self.X,_Xnew, slices1=self.Xslices,slices2=slices) - mu = np.dot(np.dot(Kx.T,self.Ki),self.likelihood.Y) - KiKx = np.dot(self.Ki,Kx) + #mu = np.dot(np.dot(Kx.T,self.Ki),self.likelihood.Y) + tmp,info = linalg.lapack.flapack.dpotrs(self.L,np.asfortranarray(self.likelihood.Y),lower=1) + mu = np.dot(Kx.T,tmp) + #KiKx = np.dot(self.Ki,Kx) + KiKx,info = linalg.lapack.flapack.dpotrs(self.L,np.asfortranarray(Kx),lower=1) if full_cov: Kxx = self.kern.K(_Xnew, slices1=slices,slices2=slices) var = Kxx - np.dot(KiKx.T,Kx)