diff --git a/GPy/inference/SGD.py b/GPy/inference/SGD.py index 701a8c65..7cb7566f 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, learning_rate_adaptation=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, actual_iter=None, schedule=None, **kwargs): self.opt_name = "Stochastic Gradient Descent" self.model = model @@ -34,12 +34,14 @@ class opt_SGD(Optimizer): self.param_traces = [('noise',[])] self.iteration_file = iteration_file self.learning_rate_adaptation = learning_rate_adaptation + self.actual_iter = actual_iter 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() + self.schedule = schedule # 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: @@ -224,48 +226,36 @@ class opt_SGD(Optimizer): return f, step, self.model.N - def adapt_learning_rate(self, t): + def adapt_learning_rate(self, t, D): 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 + if t > 0: + g_k = self.model.grads + self.s_k += np.square(g_k) + t0 = 100.0 + self.learning_rate = 0.1/(t0 + np.sqrt(self.s_k)) + + import pdb; pdb.set_trace() else: self.learning_rate = np.zeros_like(self.learning_rate) + self.s_k = np.zeros_like(self.x_opt) + elif self.learning_rate_adaptation == 'annealing': - self.learning_rate = self.learning_rate_0/(1+float(t+1)/10) + #self.learning_rate = self.learning_rate_0/(1+float(t+1)/10) + self.learning_rate = np.ones_like(self.learning_rate) * self.schedule[t] + + elif self.learning_rate_adaptation == 'semi_pesky': if self.model.__class__.__name__ == 'Bayesian_GPLVM': + g_t = self.model.grads if t == 0: - N = self.model.N - Q = self.model.Q - M = self.model.M + self.hbar_t = 0.0 + self.tau_t = 100.0 + self.gbar_t = 0.0 - 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 - - 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] + 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.ones_like(self.learning_rate)*(np.dot(self.gbar_t.T, self.gbar_t) / self.hbar_t) + tau_t = self.tau_t*(1-self.learning_rate) + 1 def opt(self, f_fp=None, f=None, fp=None): @@ -274,8 +264,8 @@ class opt_SGD(Optimizer): X, Y = self.model.X.copy(), self.model.likelihood.Y.copy() - self.model.likelihood.YYT = None - self.model.likelihood.trYYT = None + self.model.likelihood.YYT = 0 + self.model.likelihood.trYYT = 0 self.model.likelihood._bias = 0.0 self.model.likelihood._scale = 1.0 @@ -287,6 +277,9 @@ class opt_SGD(Optimizer): step = np.zeros_like(num_params) for it in range(self.iterations): + if self.actual_iter != None: + it = self.actual_iter + self.model.grads = np.zeros_like(self.x_opt) # TODO this is ugly if it == 0 or self.self_paced is False: @@ -316,6 +309,7 @@ class opt_SGD(Optimizer): self.model.likelihood.trYYT = np.trace(self.model.likelihood.YYT) Nj = N f, fp = f_fp(self.x_opt) + self.model.grads = fp.copy() step = self.momentum * step + self.learning_rate * fp self.x_opt -= step @@ -326,6 +320,7 @@ class opt_SGD(Optimizer): sys.stdout.flush() self.param_traces['noise'].append(noise) + self.adapt_learning_rate(it+count, D) NLL.append(f) self.fopt_trace.append(NLL[-1]) # fig = plt.figure('traces') @@ -335,7 +330,6 @@ class opt_SGD(Optimizer): # 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 diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index f5a3c2e2..74f526d0 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -41,8 +41,40 @@ class MRD(model): :param kernel: kernel to use """ - #TODO allow different kernels for different outputs - #def __init__(self, *Ylist, **kwargs): + def __init__(self,likelihood_list,Q,M=10,names=None,kernels=None,initX='PCA',initz='permute',_debug=False, **kwargs): + if names is None: + self.names = ["{}".format(i + 1) for i in range(len(likelihood_list))] + + #sort out the kernels + if kernels is None: + kernels = [None]*len(likelihood_list) + elif isinstance(kernels,kern.kern): + kernels = [kernels.copy() for i in range(len(likelihood_list))] + else: + assert len(kernels)==len(likelihood_list), "need one kernel per output" + assert all([isinstance(k, kern.kern) for k in kernels]), "invalid kernel object detected!" + + self.Q = Q + self.M = M + self.N = self.gref.N + self.NQ = self.N * self.Q + self.MQ = self.M * self.Q + + self._init = True + X = self._init_X(initx, likelihood_list) + Z = self._init_Z(initz, X) + self.bgplvms = [Bayesian_GPLVM(l, k, X=X, Z=Z, M=self.M, **kwargs) for l,k in zip(likelihood_list,kernels)] + + del self._init + + self.gref = self.bgplvms[0] + nparams = numpy.array([0] + [sparse_GP._get_params(g).size - g.Z.size for g in self.bgplvms]) + self.nparams = nparams.cumsum() + + model.__init__(self) # @UndefinedVariable + + + def __init__(self, *likelihood_list, **kwargs): if kwargs.has_key("_debug"): self._debug = kwargs['_debug'] diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 8605e2a8..93877cbc 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -234,7 +234,7 @@ class sparse_GP(GP): Kxx = self.kern.Kdiag(Xnew, which_parts=which_parts) var = Kxx - np.sum(Kx * np.dot(Kmmi_LmiBLmi, Kx), 0) else: - assert which_parts=='all', "swithching out parts of variational kernels is not implemented" + # assert which_parts=='all', "swithching out parts of variational kernels is not implemented" Kx = self.kern.psi1(self.Z, Xnew, X_variance_new)#, which_parts=which_parts) TODO: which_parts mu = np.dot(Kx, self.Cpsi1V) if full_cov: diff --git a/GPy/util/visualize.py b/GPy/util/visualize.py index 4ecf6c4f..2b515e91 100644 --- a/GPy/util/visualize.py +++ b/GPy/util/visualize.py @@ -86,7 +86,6 @@ class lvm(data_show): def modify(self, vals): """When latent values are modified update the latent representation and ulso update the output visualization.""" y = self.model.predict(vals)[0] - print y self.data_visualize.modify(y) self.latent_handle.set_data(vals[self.latent_index[0]], vals[self.latent_index[1]]) self.axes.figure.canvas.draw()