From cb082898d335aea8be110e615a9666065669da4e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 18:56:37 +0000 Subject: [PATCH 1/4] added trace_sum for efficiency --- GPy/models/sparse_GP.py | 6 +++--- GPy/util/linalg.py | 7 +++++++ 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index f70938c9..4846bf8a 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -3,7 +3,7 @@ import numpy as np import pylab as pb -from ..util.linalg import mdot, jitchol, chol_inv, pdinv +from ..util.linalg import mdot, jitchol, chol_inv, pdinv, trace_dot from ..util.plot import gpplot from .. import kern from GP import GP @@ -107,7 +107,7 @@ class sparse_GP(GP): self.C = np.dot(tmp,tmp.T) #self.C = mdot(self.Lmi.T, self.Bi, self.Lmi) #self.E = mdot(self.C, self.psi1VVpsi1/sf2, self.C.T) - tmp = np.dot(self.C,self.psi1V/sf) + tmp = np.dot(self.C/sf,self.psi1V) self.E = np.dot(tmp,tmp.T) # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertin inputs case @@ -156,7 +156,7 @@ class sparse_GP(GP): 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 * np.sum(self.Bi*self.A)/beta + dbeta += - 0.5 * self.D * trace_dot(self.Bi,self.A)/beta dbeta += np.sum((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) * self.psi1VVpsi1 )/beta self.partial_for_likelihood = -dbeta*self.likelihood.precision**2 diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 26105789..cf023284 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -14,6 +14,13 @@ import types #import scipy.lib.lapack.flapack import scipy as sp +def trace_dot(a,b): + """ + efficiently compute the trace of the matrix product of a and b + """ + assert a.shape==b.T.shape + return np.dot(a.flatten(),b.T.flatten()) + def mdot(*args): """Multiply all the arguments using matrix product rules. The output is equivalent to multiplying the arguments one by one From d3c87feffadbbc281f16cb470b0c5cae6731e923 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 11 Mar 2013 19:43:46 +0000 Subject: [PATCH 2/4] some messing with the linear algebra in sparse_GP. This should be more efficient... let's hope nothing breaks --- GPy/models/sparse_GP.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 4846bf8a..f5279f86 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -102,17 +102,16 @@ class sparse_GP(GP): self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B) self.psi1V = np.dot(self.psi1, self.V) - self.psi1VVpsi1 = np.dot(self.psi1V, self.psi1V.T) tmp = np.dot(self.Lmi.T, self.LBi.T) self.C = np.dot(tmp,tmp.T) - #self.C = mdot(self.Lmi.T, self.Bi, self.Lmi) - #self.E = mdot(self.C, self.psi1VVpsi1/sf2, self.C.T) - tmp = np.dot(self.C/sf,self.psi1V) - self.E = np.dot(tmp,tmp.T) + self.Cpsi1V = np.dot(self.C,self.psi1V) + self.Cpsi1VVpsi1 = np.dot(self.Cpsi1V,self.psi1V.T) + self.E = np.dot(self.Cpsi1VVpsi1,self.C)/sf2 # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertin inputs case self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten() - self.dL_dpsi1 = mdot(self.V, self.psi1V.T,self.C).T + #self.dL_dpsi1 = mdot(self.V, self.psi1V.T,self.C).T + self.dL_dpsi1 = np.dot(self.Cpsi1V,self.V.T) if self.likelihood.is_heteroscedastic: if self.has_uncertain_inputs: self.dL_dpsi2 = 0.5 * self.likelihood.precision[:,None,None] * self.D * self.Kmmi[None,:,:] # dB @@ -139,7 +138,7 @@ class sparse_GP(GP): # Compute dL_dKmm self.dL_dKmm = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB self.dL_dKmm += -0.5 * self.D * (- self.C/sf2 - 2.*mdot(self.C, self.psi2_beta_scaled, self.Kmmi) + self.Kmmi) # dC - self.dL_dKmm += np.dot(np.dot(self.E*sf2, self.psi2_beta_scaled) - np.dot(self.C, self.psi1VVpsi1), self.Kmmi) + 0.5*self.E # dD + self.dL_dKmm += np.dot(np.dot(self.E*sf2, self.psi2_beta_scaled) - self.Cpsi1VVpsi1, self.Kmmi) + 0.5*self.E # dD #the partial derivative vector for the likelihood if self.likelihood.Nparams ==0: @@ -157,7 +156,7 @@ class sparse_GP(GP): 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.sum((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) * self.psi1VVpsi1 )/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 @@ -198,7 +197,7 @@ class sparse_GP(GP): 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.sum(self.psi1VVpsi1 * self.C) + D = 0.5*np.trace(self.Cpsi1VVpsi1) return A+B+C+D def _log_likelihood_gradients(self): From 7e4b460cdbd4e37597efc9aa6403ed084b09e020 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 12 Mar 2013 09:18:15 +0000 Subject: [PATCH 3/4] more messing with the linear algebra in sparse_GP --- GPy/models/sparse_GP.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index f5279f86..be451c12 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -137,8 +137,9 @@ class sparse_GP(GP): # Compute dL_dKmm self.dL_dKmm = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB - self.dL_dKmm += -0.5 * self.D * (- self.C/sf2 - 2.*mdot(self.C, self.psi2_beta_scaled, self.Kmmi) + self.Kmmi) # dC - self.dL_dKmm += np.dot(np.dot(self.E*sf2, self.psi2_beta_scaled) - self.Cpsi1VVpsi1, self.Kmmi) + 0.5*self.E # dD + #self.dL_dKmm += -0.5 * self.D * (- self.C/sf2 - 2.*mdot(self.C, self.psi2_beta_scaled, self.Kmmi) + self.Kmmi) # dC + #self.dL_dKmm += np.dot(np.dot(self.E*sf2, self.psi2_beta_scaled) - self.Cpsi1VVpsi1, self.Kmmi) + 0.5*self.E # dD + self.dL_dKmm += 0.5*(self.D*(self.C/sf2 -self.Kmmi) + self.E) + np.dot(np.dot(self.D*self.C + self.E*sf2,self.psi2_beta_scaled) - self.Cpsi1VVpsi1,self.Kmmi) # d(C+D) #the partial derivative vector for the likelihood if self.likelihood.Nparams ==0: From bddeb998bf691e51514ba845cbecfa10484ff0fe Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 12 Mar 2013 10:04:02 +0000 Subject: [PATCH 4/4] typo in comments --- GPy/models/sparse_GP.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index be451c12..a348c9f4 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -80,7 +80,7 @@ class sparse_GP(GP): #The rather complex computations of psi2_beta_scaled if self.likelihood.is_heteroscedastic: - assert self.likelihood.D == 1 #TODO: what is the likelihood is heterscedatic and there are multiple independent outputs? + assert self.likelihood.D == 1 #TODO: what if the likelihood is heterscedatic and there are multiple independent outputs? if self.has_uncertain_inputs: self.psi2_beta_scaled = (self.psi2*(self.likelihood.precision.flatten().reshape(self.N,1,1)/sf2)).sum(0) else: