From 59d866907b12c408f0194defa0887ad6077ec323 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Thu, 16 May 2013 13:50:16 +0100 Subject: [PATCH 1/9] changes in SGD --- GPy/inference/SGD.py | 47 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/GPy/inference/SGD.py b/GPy/inference/SGD.py index bfc6ee15..3b967466 100644 --- a/GPy/inference/SGD.py +++ b/GPy/inference/SGD.py @@ -18,7 +18,7 @@ class opt_SGD(Optimizer): """ - def __init__(self, start, iterations = 10, learning_rate = 1e-4, momentum = 0.9, model = None, messages = False, batch_size = 1, self_paced = False, center = True, iteration_file = None, **kwargs): + def __init__(self, start, iterations = 10, learning_rate = 1e-4, momentum = 0.9, model = None, messages = False, batch_size = 1, self_paced = False, center = True, iteration_file = None, learning_rate_adaptation=None, **kwargs): self.opt_name = "Stochastic Gradient Descent" self.model = model @@ -33,6 +33,13 @@ class opt_SGD(Optimizer): self.center = center self.param_traces = [('noise',[])] self.iteration_file = iteration_file + self.learning_rate_adaptation = learning_rate_adaptation + if self.learning_rate_adaptation != None: + if self.learning_rate_adaptation == 'annealing': + self.learning_rate_0 = self.learning_rate + else: + self.learning_rate_0 = self.learning_rate.mean() + # if len([p for p in self.model.kern.parts if p.name == 'bias']) == 1: # self.param_traces.append(('bias',[])) # if len([p for p in self.model.kern.parts if p.name == 'linear']) == 1: @@ -204,6 +211,7 @@ class opt_SGD(Optimizer): ci = self.shift_constraints(j) f, fp = f_fp(self.x_opt[j]) + step[j] = self.momentum * step[j] + self.learning_rate[j] * fp self.x_opt[j] -= step[j] self.restore_constraints(ci) @@ -216,9 +224,32 @@ class opt_SGD(Optimizer): return f, step, self.model.N + def adapt_learning_rate(self, t): + if self.learning_rate_adaptation == 'adagrad': + if t > 5: + g = np.array(self.grads) + l2_g = np.sqrt(np.square(g).sum(0)) + self.learning_rate = 0.001/l2_g + else: + self.learning_rate = np.zeros_like(self.learning_rate) + elif self.learning_rate_adaptation == 'annealing': + self.learning_rate = self.learning_rate_0/(1+float(t+1)/2) + elif self.learning_rate_adaptation == 'semi_pesky': + if t == 0: + self.hbar_t = 0.0 + self.tau_t = 1000.0 + self.gbar_t = 0.0 + g_t = self.model.grads + self.gbar_t = (1-1/self.tau_t)*self.gbar_t + 1/self.tau_t * g_t + self.hbar_t = (1-1/self.tau_t)*self.hbar_t + 1/self.tau_t * np.dot(g_t.T, g_t) + self.learning_rate = np.dot(self.gbar_t.T, self.gbar_t) / self.hbar_t + self.tau_t = self.tau_t*(1-self.learning_rate) + 1 + print self.learning_rate + self.learning_rate *= np.ones_like(self.x_opt) + def opt(self, f_fp=None, f=None, fp=None): self.x_opt = self.model._get_params_transformed() - self.model.grads = np.zeros_like(self.x_opt) + self.grads = [] X, Y = self.model.X.copy(), self.model.likelihood.Y.copy() @@ -235,6 +266,7 @@ class opt_SGD(Optimizer): step = np.zeros_like(num_params) for it in range(self.iterations): + self.model.grads = np.zeros_like(self.x_opt) # TODO this is ugly if it == 0 or self.self_paced is False: features = np.random.permutation(Y.shape[1]) @@ -272,16 +304,17 @@ class opt_SGD(Optimizer): sys.stdout.write(status) sys.stdout.flush() self.param_traces['noise'].append(noise) - NLL.append(f) - self.fopt_trace.append(f) + NLL.append(f) + self.fopt_trace.append(NLL[-1]) # fig = plt.figure('traces') # plt.clf() # plt.plot(self.param_traces['noise']) # for k in self.param_traces.keys(): # self.param_traces[k].append(self.model.get(k)[0]) - + self.grads.append(self.model.grads.tolist()) + self.adapt_learning_rate(it) # should really be a sum(), but earlier samples in the iteration will have a very crappy ll self.f_opt = np.mean(NLL) self.model.N = N @@ -293,7 +326,7 @@ class opt_SGD(Optimizer): sigma = self.model.likelihood._variance self.model.likelihood._variance = None # invalidate cache self.model.likelihood._set_params(sigma) - + self.trace.append(self.f_opt) if self.iteration_file is not None: f = open(self.iteration_file + "iteration%d.pickle" % it, 'w') @@ -303,6 +336,6 @@ class opt_SGD(Optimizer): if self.messages != 0: sys.stdout.write('\r' + ' '*len(status)*2 + ' \r') - status = "SGD Iteration: {0: 3d}/{1: 3d} f: {2: 2.3f}\n".format(it+1, self.iterations, self.f_opt) + status = "SGD Iteration: {0: 3d}/{1: 3d} f: {2: 2.3f} max eta: {3: 1.5f}\n".format(it+1, self.iterations, self.f_opt, self.learning_rate.max()) sys.stdout.write(status) sys.stdout.flush() From 5183a18a1f94e3d51506b446575172e004109348 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Fri, 17 May 2013 12:26:08 +0100 Subject: [PATCH 2/9] minor SGD changes --- GPy/inference/SGD.py | 43 ++++++++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/GPy/inference/SGD.py b/GPy/inference/SGD.py index 3b967466..701a8c65 100644 --- a/GPy/inference/SGD.py +++ b/GPy/inference/SGD.py @@ -233,19 +233,40 @@ class opt_SGD(Optimizer): else: self.learning_rate = np.zeros_like(self.learning_rate) elif self.learning_rate_adaptation == 'annealing': - self.learning_rate = self.learning_rate_0/(1+float(t+1)/2) + self.learning_rate = self.learning_rate_0/(1+float(t+1)/10) elif self.learning_rate_adaptation == 'semi_pesky': - if t == 0: - self.hbar_t = 0.0 - self.tau_t = 1000.0 - self.gbar_t = 0.0 + if self.model.__class__.__name__ == 'Bayesian_GPLVM': + if t == 0: + N = self.model.N + Q = self.model.Q + M = self.model.M + + iip_pos = np.arange(2*N*Q,2*N*Q+M*Q) + mu_pos = np.arange(0,N*Q) + S_pos = np.arange(N*Q,2*N*Q) + self.vbparam_dict = {'iip': [iip_pos], + 'mu': [mu_pos], + 'S': [S_pos]} + + for k in self.vbparam_dict.keys(): + hbar_t = 0.0 + tau_t = 1000.0 + gbar_t = 0.0 + self.vbparam_dict[k].append(hbar_t) + self.vbparam_dict[k].append(tau_t) + self.vbparam_dict[k].append(gbar_t) + g_t = self.model.grads - self.gbar_t = (1-1/self.tau_t)*self.gbar_t + 1/self.tau_t * g_t - self.hbar_t = (1-1/self.tau_t)*self.hbar_t + 1/self.tau_t * np.dot(g_t.T, g_t) - self.learning_rate = np.dot(self.gbar_t.T, self.gbar_t) / self.hbar_t - self.tau_t = self.tau_t*(1-self.learning_rate) + 1 - print self.learning_rate - self.learning_rate *= np.ones_like(self.x_opt) + + for k in self.vbparam_dict.keys(): + pos, hbar_t, tau_t, gbar_t = self.vbparam_dict[k] + + gbar_t = (1-1/tau_t)*gbar_t + 1/tau_t * g_t[pos] + hbar_t = (1-1/tau_t)*hbar_t + 1/tau_t * np.dot(g_t[pos].T, g_t[pos]) + self.learning_rate[pos] = np.dot(gbar_t.T, gbar_t) / hbar_t + tau_t = tau_t*(1-self.learning_rate[pos]) + 1 + self.vbparam_dict[k] = [pos, hbar_t, tau_t, gbar_t] + def opt(self, f_fp=None, f=None, fp=None): self.x_opt = self.model._get_params_transformed() From ddd3ece3ceb225b33f0ac59215d38832fc67ca14 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Fri, 17 May 2013 12:29:13 +0100 Subject: [PATCH 3/9] test --- GPy/kern/kern.py | 50 +++++++++++++++++++++++++----------------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index c682fdcc..0e425e38 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -315,31 +315,33 @@ class kern(parameterised): # compute the "cross" terms # TODO: input_slices needed + crossterms = 0 for p1, p2 in itertools.combinations(self.parts, 2): - # white doesn;t combine with anything - if p1.name == 'white' or p2.name == 'white': - pass - # rbf X bias - elif p1.name == 'bias' and p2.name == 'rbf': - target += p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :]) - elif p2.name == 'bias' and p1.name == 'rbf': - target += p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :]) - # linear X bias - elif p1.name == 'bias' and p2.name == 'linear': - tmp = np.zeros((mu.shape[0], Z.shape[0])) - p2.psi1(Z, mu, S, tmp) - target += p1.variance * (tmp[:, :, None] + tmp[:, None, :]) - elif p2.name == 'bias' and p1.name == 'linear': - tmp = np.zeros((mu.shape[0], Z.shape[0])) - p1.psi1(Z, mu, S, tmp) - target += p2.variance * (tmp[:, :, None] + tmp[:, None, :]) - # rbf X linear - elif p1.name == 'linear' and p2.name == 'rbf': - raise NotImplementedError # TODO - elif p2.name == 'linear' and p1.name == 'rbf': - raise NotImplementedError # TODO - else: - raise NotImplementedError, "psi2 cannot be computed for this kernel" + prod = np.multiply + # # white doesn;t combine with anything + # if p1.name == 'white' or p2.name == 'white': + # pass + # # rbf X bias + # elif p1.name == 'bias' and p2.name == 'rbf': + # target += p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :]) + # elif p2.name == 'bias' and p1.name == 'rbf': + # target += p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :]) + # # linear X bias + # elif p1.name == 'bias' and p2.name == 'linear': + # tmp = np.zeros((mu.shape[0], Z.shape[0])) + # p2.psi1(Z, mu, S, tmp) + # target += p1.variance * (tmp[:, :, None] + tmp[:, None, :]) + # elif p2.name == 'bias' and p1.name == 'linear': + # tmp = np.zeros((mu.shape[0], Z.shape[0])) + # p1.psi1(Z, mu, S, tmp) + # target += p2.variance * (tmp[:, :, None] + tmp[:, None, :]) + # # rbf X linear + # elif p1.name == 'linear' and p2.name == 'rbf': + # raise NotImplementedError # TODO + # elif p2.name == 'linear' and p1.name == 'rbf': + # raise NotImplementedError # TODO + # else: + # raise NotImplementedError, "psi2 cannot be computed for this kernel" return target def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S): From 19407293dcaf7b642ab4335fea985aa7ee398ec7 Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Fri, 17 May 2013 16:29:32 +0100 Subject: [PATCH 4/9] cross-terms --- GPy/examples/dimensionality_reduction.py | 4 +- GPy/kern/bias.py | 5 +- GPy/kern/kern.py | 130 ++++++----------------- GPy/kern/kernpart.py | 2 - GPy/kern/white.py | 4 +- 5 files changed, 40 insertions(+), 105 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index dcda4f42..9b51947b 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -17,11 +17,11 @@ def BGPLVM(seed=default_seed): D = 4 # generate GPLVM-like data X = np.random.rand(N, Q) - k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) + k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001) K = k.K(X) Y = np.random.multivariate_normal(np.zeros(N), K, D).T - k = GPy.kern.linear(Q, ARD=True) + GPy.kern.white(Q) + k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True) + GPy.kern.white(Q) # k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q) # k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) # k = GPy.kern.rbf(Q, ARD = False) + GPy.kern.white(Q, 0.00001) diff --git a/GPy/kern/bias.py b/GPy/kern/bias.py index b5883f87..09f0afa9 100644 --- a/GPy/kern/bias.py +++ b/GPy/kern/bias.py @@ -55,8 +55,9 @@ class bias(kernpart): target += self.variance def psi1(self, Z, mu, S, target): - target += self.variance - + self._psi1 = self.variance + target += self._psi1 + def psi2(self, Z, mu, S, target): target += self.variance**2 diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index 0e425e38..c9582ac8 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -316,32 +316,19 @@ class kern(parameterised): # compute the "cross" terms # TODO: input_slices needed crossterms = 0 + for p1, p2 in itertools.combinations(self.parts, 2): - prod = np.multiply - # # white doesn;t combine with anything - # if p1.name == 'white' or p2.name == 'white': - # pass - # # rbf X bias - # elif p1.name == 'bias' and p2.name == 'rbf': - # target += p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :]) - # elif p2.name == 'bias' and p1.name == 'rbf': - # target += p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :]) - # # linear X bias - # elif p1.name == 'bias' and p2.name == 'linear': - # tmp = np.zeros((mu.shape[0], Z.shape[0])) - # p2.psi1(Z, mu, S, tmp) - # target += p1.variance * (tmp[:, :, None] + tmp[:, None, :]) - # elif p2.name == 'bias' and p1.name == 'linear': - # tmp = np.zeros((mu.shape[0], Z.shape[0])) - # p1.psi1(Z, mu, S, tmp) - # target += p2.variance * (tmp[:, :, None] + tmp[:, None, :]) - # # rbf X linear - # elif p1.name == 'linear' and p2.name == 'rbf': - # raise NotImplementedError # TODO - # elif p2.name == 'linear' and p1.name == 'rbf': - # raise NotImplementedError # TODO - # else: - # raise NotImplementedError, "psi2 cannot be computed for this kernel" + + # TODO psi1 this must be faster/better/precached/more nice + tmp1 = np.zeros((mu.shape[0], Z.shape[0])) + p1.psi1(Z, mu, S, tmp1) + tmp2 = np.zeros((mu.shape[0], Z.shape[0])) + p2.psi1(Z, mu, S, tmp2) + + prod = np.multiply(tmp1, tmp2) + crossterms += prod[:,:,None] + prod[:, None, :] + + target += crossterms return target def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S): @@ -350,71 +337,34 @@ class kern(parameterised): # compute the "cross" terms # TODO: better looping, input_slices - for i1, i2 in itertools.combinations(range(len(self.parts)), 2): + for i1, i2 in itertools.permutations(range(len(self.parts)), 2): p1, p2 = self.parts[i1], self.parts[i2] # ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2] ps1, ps2 = self.param_slices[i1], self.param_slices[i2] - # white doesn;t combine with anything - if p1.name == 'white' or p2.name == 'white': - pass - # rbf X bias - elif p1.name == 'bias' and p2.name == 'rbf': - p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2]) - p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2._psi1 * 2., Z, mu, S, target[ps1]) - elif p2.name == 'bias' and p1.name == 'rbf': - p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target[ps1]) - p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1._psi1 * 2., Z, mu, S, target[ps2]) - # linear X bias - elif p1.name == 'bias' and p2.name == 'linear': - p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2]) # [ps1]) - psi1 = np.zeros((mu.shape[0], Z.shape[0])) - p2.psi1(Z, mu, S, psi1) - p1.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps1]) - elif p2.name == 'bias' and p1.name == 'linear': - p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target[ps1]) - psi1 = np.zeros((mu.shape[0], Z.shape[0])) - p1.psi1(Z, mu, S, psi1) - p2.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps2]) - # rbf X linear - elif p1.name == 'linear' and p2.name == 'rbf': - raise NotImplementedError # TODO - elif p2.name == 'linear' and p1.name == 'rbf': - raise NotImplementedError # TODO - else: - raise NotImplementedError, "psi2 cannot be computed for this kernel" + tmp = np.zeros((mu.shape[0], Z.shape[0])) + p1.psi1(Z, mu, S, tmp) + p2.dpsi1_dtheta((tmp[:,None,:]*dL_dpsi2).sum(1)*2., Z, mu, S, target[ps2]) return self._transform_gradients(target) def dpsi2_dZ(self, dL_dpsi2, Z, mu, S): target = np.zeros_like(Z) [p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] + #target *= 2 # compute the "cross" terms # TODO: we need input_slices here. - for p1, p2 in itertools.combinations(self.parts, 2): - # white doesn;t combine with anything - if p1.name == 'white' or p2.name == 'white': - pass - # rbf X bias - elif p1.name == 'bias' and p2.name == 'rbf': - p2.dpsi1_dX(dL_dpsi2.sum(1).T * p1.variance, Z, mu, S, target) - elif p2.name == 'bias' and p1.name == 'rbf': - p1.dpsi1_dZ(dL_dpsi2.sum(1).T * p2.variance, Z, mu, S, target) - # linear X bias - elif p1.name == 'bias' and p2.name == 'linear': - p2.dpsi1_dZ(dL_dpsi2.sum(1).T * p1.variance, Z, mu, S, target) - elif p2.name == 'bias' and p1.name == 'linear': - p1.dpsi1_dZ(dL_dpsi2.sum(1).T * p2.variance, Z, mu, S, target) - # rbf X linear - elif p1.name == 'linear' and p2.name == 'rbf': - raise NotImplementedError # TODO - elif p2.name == 'linear' and p1.name == 'rbf': - raise NotImplementedError # TODO - else: - raise NotImplementedError, "psi2 cannot be computed for this kernel" + for p1, p2 in itertools.permutations(self.parts, 2): + if p1.name == 'linear' and p2.name == 'linear': + raise NotImplementedError("We don't handle linear/linear cross-terms") + tmp = np.zeros((mu.shape[0], Z.shape[0])) + p1.psi1(Z, mu, S, tmp) + tmp2 = np.zeros_like(target) + p2.dpsi1_dZ((tmp[:,None,:]*dL_dpsi2).sum(1).T, Z, mu, S, tmp2) + target += tmp2 - return target * 2. + return target * 2 def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S): target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1])) @@ -422,27 +372,13 @@ class kern(parameterised): # compute the "cross" terms # TODO: we need input_slices here. - for p1, p2 in itertools.combinations(self.parts, 2): - # white doesn;t combine with anything - if p1.name == 'white' or p2.name == 'white': - pass - # rbf X bias - elif p1.name == 'bias' and p2.name == 'rbf': - p2.dpsi1_dmuS(dL_dpsi2.sum(1).T * p1.variance * 2., Z, mu, S, target_mu, target_S) - elif p2.name == 'bias' and p1.name == 'rbf': - p1.dpsi1_dmuS(dL_dpsi2.sum(1).T * p2.variance * 2., Z, mu, S, target_mu, target_S) - # linear X bias - elif p1.name == 'bias' and p2.name == 'linear': - p2.dpsi1_dmuS(dL_dpsi2.sum(1).T * p1.variance * 2., Z, mu, S, target_mu, target_S) - elif p2.name == 'bias' and p1.name == 'linear': - p1.dpsi1_dmuS(dL_dpsi2.sum(1).T * p2.variance * 2., Z, mu, S, target_mu, target_S) - # rbf X linear - elif p1.name == 'linear' and p2.name == 'rbf': - raise NotImplementedError # TODO - elif p2.name == 'linear' and p1.name == 'rbf': - raise NotImplementedError # TODO - else: - raise NotImplementedError, "psi2 cannot be computed for this kernel" + for p1, p2 in itertools.permutations(self.parts, 2): + if p1.name == 'linear' and p2.name == 'linear': + raise NotImplementedError("We don't handle linear/linear cross-terms") + + tmp = np.zeros((mu.shape[0], Z.shape[0])) + p1.psi1(Z, mu, S, tmp) + p2.dpsi1_dmuS((tmp[:,None,:]*dL_dpsi2).sum(1).T*2., Z, mu, S, target_mu, target_S) return target_mu, target_S diff --git a/GPy/kern/kernpart.py b/GPy/kern/kernpart.py index 30a1cc3d..7de150e9 100644 --- a/GPy/kern/kernpart.py +++ b/GPy/kern/kernpart.py @@ -54,5 +54,3 @@ class kernpart(object): raise NotImplementedError def dK_dX(self,X,X2,target): raise NotImplementedError - - diff --git a/GPy/kern/white.py b/GPy/kern/white.py index be6aad45..d5701cd9 100644 --- a/GPy/kern/white.py +++ b/GPy/kern/white.py @@ -18,7 +18,8 @@ class white(kernpart): self.Nparam = 1 self.name = 'white' self._set_params(np.array([variance]).flatten()) - + self._psi1 = 0 # TODO: more elegance here + def _get_params(self): return self.variance @@ -81,4 +82,3 @@ class white(kernpart): def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S): pass - From 6b63237adf5568a2cdaf779d6d4ce922132bd224 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 17 May 2013 17:17:30 +0100 Subject: [PATCH 5/9] SCG printing prettyfied --- GPy/examples/classification.py | 38 +++++++++++++++++----------------- GPy/inference/SCG.py | 25 +++++++++++++--------- 2 files changed, 34 insertions(+), 29 deletions(-) diff --git a/GPy/examples/classification.py b/GPy/examples/classification.py index d6697d7c..9168db7c 100644 --- a/GPy/examples/classification.py +++ b/GPy/examples/classification.py @@ -9,8 +9,8 @@ import pylab as pb import numpy as np import GPy -default_seed=10000 -def crescent_data(seed=default_seed): #FIXME +default_seed = 10000 +def crescent_data(seed=default_seed): # FIXME """Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood. :param model_type: type of model to fit ['Full', 'FITC', 'DTC']. @@ -27,10 +27,10 @@ def crescent_data(seed=default_seed): #FIXME # Likelihood object distribution = GPy.likelihoods.likelihood_functions.probit() - likelihood = GPy.likelihoods.EP(data['Y'],distribution) + likelihood = GPy.likelihoods.EP(data['Y'], distribution) - m = GPy.models.GP(data['X'],likelihood,kernel) + m = GPy.models.GP(data['X'], likelihood, kernel) m.ensure_default_constraints() m.update_likelihood_approximation() @@ -54,10 +54,10 @@ def oil(): # Likelihood object distribution = GPy.likelihoods.likelihood_functions.probit() - likelihood = GPy.likelihoods.EP(data['Y'][:, 0:1],distribution) + likelihood = GPy.likelihoods.EP(data['Y'][:, 0:1], distribution) # Create GP model - m = GPy.models.GP(data['X'],likelihood=likelihood,kernel=kernel) + m = GPy.models.GP(data['X'], likelihood=likelihood, kernel=kernel) # Contrain all parameters to be positive m.constrain_positive('') @@ -85,17 +85,17 @@ def toy_linear_1d_classification(seed=default_seed): # Likelihood object distribution = GPy.likelihoods.likelihood_functions.probit() - likelihood = GPy.likelihoods.EP(Y,distribution) + likelihood = GPy.likelihoods.EP(Y, distribution) # Model definition - m = GPy.models.GP(data['X'],likelihood=likelihood,kernel=kernel) + m = GPy.models.GP(data['X'], likelihood=likelihood, kernel=kernel) m.ensure_default_constraints() # Optimize m.update_likelihood_approximation() # Parameters optimization: m.optimize() - #m.pseudo_EM() #FIXME + # m.pseudo_EM() #FIXME # Plot pb.subplot(211) @@ -121,20 +121,20 @@ def sparse_toy_linear_1d_classification(seed=default_seed): # Likelihood object distribution = GPy.likelihoods.likelihood_functions.probit() - likelihood = GPy.likelihoods.EP(Y,distribution) + likelihood = GPy.likelihoods.EP(Y, distribution) - Z = np.random.uniform(data['X'].min(),data['X'].max(),(10,1)) + Z = np.random.uniform(data['X'].min(), data['X'].max(), (10, 1)) # Model definition - m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z,normalize_X=False) - m.set('len',2.) + m = GPy.models.sparse_GP(data['X'], likelihood=likelihood, kernel=kernel, Z=Z, normalize_X=False) + m.set('len', 2.) m.ensure_default_constraints() # Optimize m.update_likelihood_approximation() # Parameters optimization: m.optimize() - #m.EPEM() #FIXME + # m.EPEM() #FIXME # Plot pb.subplot(211) @@ -162,15 +162,15 @@ def sparse_crescent_data(inducing=10, seed=default_seed): # Likelihood object distribution = GPy.likelihoods.likelihood_functions.probit() - likelihood = GPy.likelihoods.EP(data['Y'],distribution) + likelihood = GPy.likelihoods.EP(data['Y'], distribution) - sample = np.random.randint(0,data['X'].shape[0],inducing) - Z = data['X'][sample,:] + sample = np.random.randint(0, data['X'].shape[0], inducing) + Z = data['X'][sample, :] # create sparse GP EP model - m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z) + m = GPy.models.sparse_GP(data['X'], likelihood=likelihood, kernel=kernel, Z=Z) m.ensure_default_constraints() - m.set('len',10.) + m.set('len', 10.) m.update_likelihood_approximation() diff --git a/GPy/inference/SCG.py b/GPy/inference/SCG.py index f9520813..e6ef25c0 100644 --- a/GPy/inference/SCG.py +++ b/GPy/inference/SCG.py @@ -39,6 +39,9 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto function_eval number of fn evaluations status: string describing convergence status """ + print " SCG" + print ' {0:{mi}s} {1:11s} {2:11s} {3:11s}'.format("I", "F", "Scale", "|g|", mi=len(str(maxiters))) + if xtol is None: xtol = 1e-6 if ftol is None: @@ -46,18 +49,18 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto if gtol is None: gtol = 1e-5 sigma0 = 1.0e-4 - fold = f(x, *optargs) # Initial function value. + fold = f(x, *optargs) # Initial function value. function_eval = 1 fnow = fold - gradnew = gradf(x, *optargs) # Initial gradient. + 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. - nsuccess = 0 # nsuccess counts number of successes. - beta = 1.0 # Initial scale parameter. - betamin = 1.0e-15 # Lower bound on scale. - betamax = 1.0e100 # Upper bound on scale. + d = -gradnew # Initial search direction. + success = True # Force calculation of directional derivs. + nsuccess = 0 # nsuccess counts number of successes. + beta = 1.0 # Initial scale parameter. + betamin = 1.0e-15 # Lower bound on scale. + betamax = 1.0e100 # Upper bound on scale. status = "Not converged" flog = [fold] @@ -82,6 +85,8 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto # Increase effective curvature and evaluate step size alpha. delta = theta + beta * kappa if delta <= 0: + if display: + print "" delta = beta * kappa beta = beta - theta / kappa @@ -107,12 +112,12 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto fnow = fold # Store relevant variables - flog.append(fnow) # Current function value + flog.append(fnow) # Current function value iteration += 1 if display: print '\r', - print 'Iter: {0:>0{mi}g} Obj:{1:> 12e} Scale:{2:> 12e} |g|:{3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len(str(maxiters))), + print '{0:>0{mi}g} {1:> 12e} {2:> 12e} {3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len(str(maxiters))), # print 'Iteration:', iteration, ' Objective:', fnow, ' Scale:', beta, '\r', sys.stdout.flush() From 58b206a24541fe55183676d11b225116a61c9bb7 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Sun, 19 May 2013 19:16:25 +0100 Subject: [PATCH 6/9] enabled DSYR on feeble machines i.e. where numpy is compiled without proper blas linkage --- GPy/util/linalg.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index fa7de8c1..abfb1900 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -236,7 +236,7 @@ def tdot(*args, **kwargs): else: return tdot_numpy(*args,**kwargs) -def DSYR(A,x,alpha=1.): +def DSYR_blas(A,x,alpha=1.): """ Performs a symmetric rank-1 update operation: A <- A + alpha * np.dot(x,x.T) @@ -258,6 +258,26 @@ def DSYR(A,x,alpha=1.): x_, byref(INCX), A_, byref(LDA)) symmetrify(A,upper=True) +def DSYR_numpy(A,x,alpha=1.): + """ + Performs a symmetric rank-1 update operation: + A <- A + alpha * np.dot(x,x.T) + + Arguments + --------- + :param A: Symmetric NxN np.array + :param x: Nx1 np.array + :param alpha: scalar + """ + A += alpha*np.dot(x[:,None],x[None,:]) + + +def DSYR(*args, **kwargs): + if _blas_available: + return DSYR_blas(*args,**kwargs) + else: + return DSYR_numpy(*args,**kwargs) + def symmetrify(A,upper=False): """ Take the square matrix A and make it symmetrical by copting elements from the lower half to the upper From 3f22e61d2d09f17684cfc35a8b5d62ca9f273ad3 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Sun, 19 May 2013 19:17:55 +0100 Subject: [PATCH 7/9] tidied the computation of A in sparse_GP --- GPy/models/sparse_GP.py | 41 ++++++++++++----------------------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 2aafab16..ef7b445d 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -70,39 +70,22 @@ class sparse_GP(GP): self.Lm = jitchol(self.Kmm) # The rather complex computations of self.A - if self.likelihood.is_heteroscedastic: - assert self.likelihood.D == 1 # TODO: what if the likelihood is heterscedatic and there are multiple independent outputs? - if self.has_uncertain_inputs: -# psi2_beta_scaled = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1) / sf2)).sum(0) - psi2_beta_scaled = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0) - evals, evecs = linalg.eigh(psi2_beta_scaled) - clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable - if not np.allclose(evals, clipped_evals): - print "Warning: clipping posterior eigenvalues" - tmp = evecs * np.sqrt(clipped_evals) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) - self.A = tdot(tmp) + if self.has_uncertain_inputs: + if self.likelihood.is_heteroscedastic: + psi2_beta = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0) else: -# tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)) / sf) - tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N))) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) - self.A = tdot(tmp) + psi2_beta = self.psi2.sum(0) * self.likelihood.precision + evals, evecs = linalg.eigh(psi2_beta) + clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable + tmp = evecs * np.sqrt(clipped_evals) else: - if self.has_uncertain_inputs: -# psi2_beta_scaled = (self.psi2 * (self.likelihood.precision / sf2)).sum(0) - psi2_beta_scaled = (self.psi2 * (self.likelihood.precision)).sum(0) - evals, evecs = linalg.eigh(psi2_beta_scaled) - clipped_evals = np.clip(evals, 0., 1e15) # TODO: make clipping configurable - if not np.allclose(evals, clipped_evals): - print "Warning: clipping posterior eigenvalues" - tmp = evecs * np.sqrt(clipped_evals) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) - self.A = tdot(tmp) + if self.likelihood.is_heteroscedastic: + tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N))) else: - # 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) + tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) + self.A = tdot(tmp) + # factor B # self.B = np.eye(self.M) / sf2 + self.A From 025f31a43fb51423732122f3e1c0121e85016bf9 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 20 May 2013 09:38:08 +0100 Subject: [PATCH 8/9] more tidying in sparse_GP --- GPy/models/sparse_GP.py | 29 +++++------------------------ 1 file changed, 5 insertions(+), 24 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index ef7b445d..aeff83c2 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -16,9 +16,9 @@ class sparse_GP(GP): :param X: inputs :type X: np.ndarray (N x Q) :param likelihood: a likelihood instance, containing the observed data - :type likelihood: GPy.likelihood.(Gaussian | EP) - :param kernel : the kernel/covariance function. See link kernels - :type kernel: a GPy kernel + :type likelihood: GPy.likelihood.(Gaussian | EP | Laplace) + :param kernel : the kernel (covariance function). See link kernels + :type kernel: a GPy.kern.kern instance :param X_variance: The uncertainty in the measurements of X (Gaussian variance) :type X_variance: np.ndarray (N x Q) | None :param Z: inducing inputs (optional, see note) @@ -30,8 +30,6 @@ class sparse_GP(GP): """ def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False): -# self.scale_factor = 100.0 # a scaling factor to help keep the algorithm stable -# self.auto_scale_factor = False self.Z = Z self.M = Z.shape[0] self.likelihood = likelihood @@ -63,8 +61,6 @@ class sparse_GP(GP): self.psi2 = None def _computations(self): -# sf = self.scale_factor -# sf2 = sf ** 2 # factor Kmm self.Lm = jitchol(self.Kmm) @@ -88,7 +84,6 @@ class sparse_GP(GP): # factor B -# self.B = np.eye(self.M) / sf2 + self.A self.B = np.eye(self.M) + self.A self.LB = jitchol(self.B) @@ -104,8 +99,6 @@ class sparse_GP(GP): # Compute dL_dKmm tmp = tdot(self._LBi_Lmi_psi1V) self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.D * np.eye(self.M) + tmp) -# tmp = -0.5 * self.DBi_plus_BiPBi / sf2 -# tmp += -0.5 * self.B * sf2 * self.D tmp = -0.5 * self.DBi_plus_BiPBi tmp += -0.5 * self.B * self.D tmp += self.D * np.eye(self.M) @@ -115,6 +108,7 @@ class sparse_GP(GP): self.dL_dpsi0 = -0.5 * self.D * (self.likelihood.precision * np.ones([self.N, 1])).flatten() self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.V.T) dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.D * np.eye(self.M) - self.DBi_plus_BiPBi) + if self.likelihood.is_heteroscedastic: if self.has_uncertain_inputs: self.dL_dpsi2 = self.likelihood.precision[:, None, None] * dL_dpsi2_beta[None, :, :] @@ -141,7 +135,6 @@ 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) 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))) @@ -149,16 +142,12 @@ class sparse_GP(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.likelihood.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)) 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) 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)) 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)) return A + B + C + D @@ -168,14 +157,6 @@ class sparse_GP(GP): self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.Nparam]) self.likelihood._set_params(p[self.Z.size + self.kern.Nparam:]) self._compute_kernel_matrices() - # if self.auto_scale_factor: - # self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision) - # if self.auto_scale_factor: - # if self.likelihood.is_heteroscedastic: - # self.scale_factor = max(100,np.sqrt(self.psi2_beta_scaled.sum(0).mean())) - # else: - # self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision) -# self.scale_factor = 100. self._computations() def _get_params(self): @@ -188,7 +169,7 @@ class sparse_GP(GP): """ Approximates a non-gaussian likelihood using Expectation Propagation - For a Gaussian (or direct: TODO) likelihood, no iteration is required: + For a Gaussian likelihood, no iteration is required: this function does nothing """ if self.has_uncertain_inputs: From 99017e5e50ff6c01bf8ae3651cfbb7e213ed1bd6 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 20 May 2013 10:11:27 +0100 Subject: [PATCH 9/9] last stability changes --- GPy/core/transformations.py | 8 +++++--- GPy/examples/dimensionality_reduction.py | 11 +++++++++-- GPy/inference/SCG.py | 2 -- GPy/models/Bayesian_GPLVM.py | 3 --- 4 files changed, 14 insertions(+), 10 deletions(-) diff --git a/GPy/core/transformations.py b/GPy/core/transformations.py index f7e59ab6..fcbfb548 100644 --- a/GPy/core/transformations.py +++ b/GPy/core/transformations.py @@ -39,8 +39,8 @@ class logexp(transformation): return '(+ve)' class logexp_clipped(transformation): - max_bound = 1e300 - min_bound = 1e-10 + max_bound = 1e250 + min_bound = 1e-9 log_max_bound = np.log(max_bound) log_min_bound = np.log(min_bound) def __init__(self, lower=1e-6): @@ -49,11 +49,13 @@ class logexp_clipped(transformation): def f(self, x): exp = np.exp(np.clip(x, self.log_min_bound, self.log_max_bound)) f = np.log(1. + exp) + if np.isnan(f).any(): + import ipdb;ipdb.set_trace() return f def finv(self, f): return np.log(np.exp(np.clip(f, self.min_bound, self.max_bound)) - 1.) def gradfactor(self, f): - ef = np.exp(f) + ef = np.exp(f) # np.clip(f, self.min_bound, self.max_bound)) gf = (ef - 1.) / ef return np.where(f < self.lower, 0, gf) def initialize(self, f): diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 18995c50..4f713c1f 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -273,8 +273,8 @@ def bgplvm_simulation(optimize='scg', pylab.figure(); pylab.axis(); m.kern.plot_ARD() return m -def mrd_simulation(plot_sim=False): - D1, D2, D3, N, M, Q = 150, 250, 300, 700, 3, 7 +def mrd_simulation(optimize=True, plot_sim=False): + D1, D2, D3, N, M, Q = 150, 250, 30, 300, 3, 7 slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim) from GPy.models import mrd @@ -292,6 +292,13 @@ def mrd_simulation(plot_sim=False): m.constrain('variance|noise', logexp_clipped()) m.ensure_default_constraints() + # DEBUG + np.seterr("raise") + + if optimize: + print "Optimizing Model:" + m.optimize('scg', messages=1, max_iters=3e3) + return m def brendan_faces(): diff --git a/GPy/inference/SCG.py b/GPy/inference/SCG.py index e6ef25c0..f190d002 100644 --- a/GPy/inference/SCG.py +++ b/GPy/inference/SCG.py @@ -85,8 +85,6 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto # Increase effective curvature and evaluate step size alpha. delta = theta + beta * kappa if delta <= 0: - if display: - print "" delta = beta * kappa beta = beta - theta / kappa diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index 5511a1b9..464d7425 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -171,9 +171,6 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): self.dbound_dZtheta = sparse_GP._log_likelihood_gradients(self) return np.hstack((self.dbound_dmuS.flatten(), self.dbound_dZtheta)) - def _log_likelihood_normal_gradients(self): - Si, _, _, _ = pdinv(self.X_variance) - def plot_latent(self, which_indices=None, *args, **kwargs): if which_indices is None: