From 83275c03e1232c3433257eb1e451331496b8088c Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 10 May 2013 11:19:27 +0100 Subject: [PATCH 1/3] termination rule update --- GPy/inference/SCG.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GPy/inference/SCG.py b/GPy/inference/SCG.py index 70b11940..d76fcf9b 100644 --- a/GPy/inference/SCG.py +++ b/GPy/inference/SCG.py @@ -49,6 +49,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto function_eval = 1 fnow = fold gradnew = gradf(x, *optargs) # Initial gradient. + current_grad = np.dot(gradnew, gradnew) gradold = gradnew.copy() d = -gradnew # Initial search direction. success = True # Force calculation of directional derivs. @@ -110,7 +111,7 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto iteration += 1 if display: print '\r', - print 'i: {0:>5g} f:{1:> 12e} b:{2:> 12e} |g|:{3:> 12e}'.format(iteration, fnow, beta, np.dot(gradnew, gradnew)), + print 'i: {0:>5g} f:{1:> 12e} b:{2:> 12e} |g|:{3:> 12e}'.format(iteration, fnow, beta, current_grad), # print 'Iteration:', iteration, ' Objective:', fnow, ' Scale:', beta, '\r', sys.stdout.flush() @@ -125,8 +126,9 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto fold = fnew gradold = gradnew gradnew = gradf(x, *optargs) + current_grad = np.dot(gradnew, gradnew) # If the gradient is zero then we are done. - if (np.dot(gradnew, gradnew)) <= gtol: + if current_grad <= gtol: status = 'converged' return x, flog, function_eval, status From 214eab5f2e495233eb5421b5c2798f80a836c4db Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 10 May 2013 11:21:19 +0100 Subject: [PATCH 2/3] MRD updates and minor changes --- GPy/models/Bayesian_GPLVM.py | 10 ++++++---- GPy/models/mrd.py | 2 +- GPy/models/sparse_GP.py | 6 +++--- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index b1b6cbcd..f2393df8 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -136,14 +136,16 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): self._savedparams.append([self.f_call, self._get_params()]) self._savedgradients.append([self.f_call, self._log_likelihood_gradients()]) self._savedpsiKmm.append([self.f_call, [self.Kmm, self.dL_dKmm]]) - sf2 = self.scale_factor ** 2 +# 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) +# B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2) + B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A)) 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 = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5 * self.M * np.log(sf2)) +# B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A) * sf2) + B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A)) + C = -self.D * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.M * np.log(sf2)) D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V)) self._savedABCD.append([self.f_call, A, B, C, D]) diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index b72f65fc..23aa81b3 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -273,7 +273,7 @@ class MRD(model): def _handle_plotting(self, fig_num, axes, plotf): if axes is None: - fig = pylab.figure(num=fig_num, figsize=(4 * len(self.bgplvms), 3 * len(self.bgplvms))) + fig = pylab.figure(num=fig_num, figsize=(4 * len(self.bgplvms), 3)) for i, g in enumerate(self.bgplvms): if axes is None: ax = fig.add_subplot(1, len(self.bgplvms), i + 1) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index f099fe53..3a45488b 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -104,7 +104,7 @@ class sparse_GP(GP): tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) self.A = tdot(tmp) else: -# tmp = self.psi1 * (np.sqrt(self.likelihood.precision) / sf) + # tmp = self.psi1 * (np.sqrt(self.likelihood.precision) / sf) tmp = self.psi1 * (np.sqrt(self.likelihood.precision)) tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) self.A = tdot(tmp) @@ -163,7 +163,7 @@ class sparse_GP(GP): else: # likelihood is not heterscedatic self.partial_for_likelihood = -0.5 * self.N * self.D * self.likelihood.precision + 0.5 * self.likelihood.trYYT * 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 * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision * sf2) self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision) self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V))) @@ -177,7 +177,7 @@ class sparse_GP(GP): # B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2) B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A)) 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 + 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) B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A)) # C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5 * self.M * np.log(sf2)) From 8a843378a0ba52d4f3ebd0d34c8ec6d7a622833a Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 10 May 2013 11:21:53 +0100 Subject: [PATCH 3/3] Example update to run oil dataset --- GPy/examples/dimensionality_reduction.py | 32 ++++++++++++------------ 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 8e836299..e7c775dd 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -63,7 +63,7 @@ def GPLVM_oil_100(optimize=True): m.plot_latent(labels=m.data_labels) return m -def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300, plot=False): +def BGPLVM_oil(optimize=True, N=100, Q=10, M=20, max_f_eval=300, plot=False): data = GPy.util.datasets.oil() # create simple GP model @@ -72,19 +72,19 @@ def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300, plot=False): m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=kernel, M=M) m.data_labels = data['Y'][:N].argmax(axis=1) + m.constrain('variance', logexp_clipped()) + m.constrain('length', logexp_clipped()) + m['lengt'] = 100. + m.ensure_default_constraints() + # optimize if optimize: - m.constrain_fixed('noise', 1. / Y.var()) - m.constrain('variance', logexp_clipped()) - m['lengt'] = 1000 + m.unconstrain('noise'); m.constrain_fixed('noise', Y.var() / 100.) + m.optimize('scg', messages=1, max_f_eval=150) - m.ensure_default_constraints() - m.optimize('scg', messages=1, max_f_eval=max(80, max_f_eval)) m.unconstrain('noise') - m.constrain_positive('noise') - m.optimize('scg', messages=1, max_f_eval=max(0, max_f_eval - 80)) - else: - m.ensure_default_constraints() + m.constrain('noise', logexp_clipped()) + m.optimize('scg', messages=1, max_f_eval=max_f_eval) if plot: y = m.likelihood.Y[0, :] @@ -92,7 +92,7 @@ def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300, plot=False): plt.sca(latent_axes) m.plot_latent() data_show = GPy.util.visualize.vector_show(y) - lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :], m, data_show, latent_axes=latent_axes) # , sense_axes=sense_axes) + lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :], m, data_show, latent_axes=latent_axes) # , sense_axes=sense_axes) raw_input('Press enter to finish') plt.close('all') # # plot @@ -182,7 +182,7 @@ def bgplvm_simulation_matlab_compare(): Y = sim_data['Y'] S = sim_data['S'] mu = sim_data['mu'] - M, [_, Q] = 20, mu.shape + M, [_, Q] = 3, mu.shape from GPy.models import mrd from GPy import kern @@ -192,7 +192,7 @@ def bgplvm_simulation_matlab_compare(): m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, # X=mu, # X_variance=S, - _debug=True) + _debug=False) m.ensure_default_constraints() m.auto_scale_factor = True m['noise'] = Y.var() / 100. @@ -223,7 +223,7 @@ def bgplvm_simulation(burnin='scg', plot_sim=False, Y = Ylist[0] - k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q) + k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q) # k = kern.white(Q, .00001) + kern.bias(Q) m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True) # m.set('noise',) @@ -358,7 +358,7 @@ def mrd_simulation(plot_sim=False): # import ipdb; ipdb.set_trace() # np.seterrcall(ipdbonerr) - return m # , mtest + return m # , mtest def mrd_silhouette(): @@ -371,7 +371,7 @@ def brendan_faces(): # optimize m.ensure_default_constraints() - # m.optimize(messages=1, max_f_eval=10000) + m.optimize(messages=1, max_f_eval=10000) ax = m.plot_latent() y = m.likelihood.Y[0, :]