From dd27b285d6376b56f77d5969dfb13961e78d5444 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 15 Apr 2013 17:53:26 +0100 Subject: [PATCH 1/5] rbf computation of psi2 now works if there's only one datum --- GPy/kern/rbf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index a654cd0f..9ff7a93e 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -239,13 +239,13 @@ class rbf(kernpart): psi2 = np.empty((N,M,M)) psi2_Zdist_sq = self._psi2_Zdist_sq - half_log_psi2_denom = 0.5*np.log(self._psi2_denom).squeeze() + _psi2_denom = self._psi2_denom.squeeze().reshape(N,self.D) + half_log_psi2_denom = 0.5*np.log(self._psi2_denom).squeeze().reshape(N,self.D) variance_sq = float(np.square(self.variance)) if self.ARD: lengthscale2 = self.lengthscale2 else: lengthscale2 = np.ones(Q)*self.lengthscale2 - _psi2_denom = self._psi2_denom.squeeze() code = """ double tmp; From 403dc4246299efcba9a865e7f95587df1e2b6dd3 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 16 Apr 2013 11:06:21 +0100 Subject: [PATCH 2/5] improved stability of sparse GP computations Specifically in computing dL_dKmm --- GPy/models/sparse_GP.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 7da36835..32e868a5 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -99,6 +99,7 @@ class sparse_GP(GP): self.V = (self.likelihood.precision/self.scale_factor)*self.likelihood.Y #Compute A = L^-1 psi2 beta L^-T + #self. A = mdot(self.Lmi,self.psi2_beta_scaled,self.Lmi.T) tmp = linalg.lapack.flapack.dtrtrs(self.Lm,self.psi2_beta_scaled.T,lower=1)[0] self.A = linalg.lapack.flapack.dtrtrs(self.Lm,np.asarray(tmp.T,order='F'),lower=1)[0] @@ -142,7 +143,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_old = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB + tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.A),lower=1,trans=1)[0] + self.dL_dKmm = -0.5*self.D*sf2*linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0] #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) From 445e6c1d30b1eea7df9de4959fe5314573046008 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 16 Apr 2013 11:09:33 +0100 Subject: [PATCH 3/5] comments only --- GPy/models/sparse_GP.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 32e868a5..4d9edacc 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -144,10 +144,10 @@ class sparse_GP(GP): # Compute dL_dKmm #self.dL_dKmm_old = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB - tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.A),lower=1,trans=1)[0] - self.dL_dKmm = -0.5*self.D*sf2*linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0] #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 + tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.A),lower=1,trans=1)[0] + self.dL_dKmm = -0.5*self.D*sf2*linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0] #dA 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 From 9acc6e972313f576aef2f43d19ccfe467f02b928 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 16 Apr 2013 11:13:53 +0100 Subject: [PATCH 4/5] changed printing behaviour in cholesky to kill last line --- GPy/util/linalg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index d82bb50f..f88099a4 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -61,7 +61,7 @@ def jitchol(A,maxtries=5): raise linalg.LinAlgError, "not pd: negative diagonal elements" jitter= diagA.mean()*1e-6 for i in range(1,maxtries+1): - print 'Warning: adding jitter of '+str(jitter) + print '\rWarning: adding jitter of {:.10e} '.format(jitter), try: return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True) except: @@ -89,7 +89,7 @@ def jitchol_old(A,maxtries=5): raise linalg.LinAlgError, "not pd: negative diagonal elements" jitter= diagA.mean()*1e-6 for i in range(1,maxtries+1): - print 'Warning: adding jitter of '+str(jitter) + print '\rWarning: adding jitter of {:.10e} '.format(jitter), try: return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True) except: From 350497c72606f188f83b68588140f0058190559b Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 16 Apr 2013 11:25:51 +0100 Subject: [PATCH 5/5] adjusted plotting behaviour in X1d --- GPy/examples/dimensionality_reduction.py | 62 +++++++++---- GPy/models/mrd.py | 112 +++++++++++++++++------ 2 files changed, 128 insertions(+), 46 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 283bdc76..2c7d6bea 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -118,13 +118,13 @@ def mrd_simulation(plot_sim=False): # Y2 -= Y2.mean(0) # make_params = lambda ard: np.hstack([[1], ard, [1, .3]]) - D1, D2, D3, N, M, Q = 6, 7, 8, 150, 18, 5 - x = np.linspace(0, 2 * np.pi, N)[:, None] + D1, D2, D3, N, M, Q = 50, 100, 8, 200, 2, 5 + x = np.linspace(0, 8 * np.pi, N)[:, None] s1 = np.vectorize(lambda x: np.sin(x)) s2 = np.vectorize(lambda x: np.cos(x)) s3 = np.vectorize(lambda x:-np.exp(-np.cos(2 * x))) - sS = np.vectorize(lambda x: np.sin(2 * x)) + sS = np.vectorize(lambda x: x * np.sin(2 * x)) s1 = s1(x) s2 = s2(x) @@ -144,20 +144,30 @@ def mrd_simulation(plot_sim=False): S2 = np.hstack([s2, sS]) S3 = np.hstack([s3, sS]) + from GPy.models import mrd + from GPy import kern + reload(mrd); reload(kern) + +# k = kern.rbf(2, ARD=True) + kern.bias(2) + kern.white(2) +# Y1 = np.random.multivariate_normal(np.zeros(N), k.K(S1), D1).T +# Y2 = np.random.multivariate_normal(np.zeros(N), k.K(S2), D2).T +# Y3 = np.random.multivariate_normal(np.zeros(N), k.K(S3), D3).T + Y1 = S1.dot(np.random.randn(S1.shape[1], D1)) Y2 = S2.dot(np.random.randn(S2.shape[1], D2)) Y3 = S3.dot(np.random.randn(S3.shape[1], D3)) - Y1 += .1 * np.random.randn(*Y1.shape) - Y2 += .1 * np.random.randn(*Y2.shape) - Y3 += .1 * np.random.randn(*Y3.shape) + Y1 += .5 * np.random.randn(*Y1.shape) + Y2 += .5 * np.random.randn(*Y2.shape) + Y3 += .5 * np.random.randn(*Y3.shape) - Y1 -= Y1.mean(0) - Y2 -= Y2.mean(0) - Y3 -= Y3.mean(0) - Y1 /= Y1.std(0) - Y2 /= Y2.std(0) - Y3 /= Y3.std(0) +# Y1 -= Y1.mean(0) +# Y2 -= Y2.mean(0) +# Y3 -= Y3.mean(0) + + # Y1 /= Y1.std(0) + # Y2 /= Y2.std(0) + # Y3 /= Y3.std(0) Slist = [s1, s2, sS] Ylist = [Y1, Y2] @@ -173,21 +183,33 @@ def mrd_simulation(plot_sim=False): ax.plot(x, S, label=lab) ax.legend() for i, Y in enumerate(Ylist): - ax = fig.add_subplot(2, len(Ylist), len(Slist) + i) + ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) ax.imshow(Y) ax.set_title("Y{}".format(i + 1)) pylab.draw() pylab.tight_layout() - from GPy.models import mrd - from GPy import kern - reload(mrd); reload(kern) - k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q) - m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, init="single", _debug=False) + # k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q) + k = kern.linear(Q, ARD=True) + kern.bias(Q) + kern.white(Q) + m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", _debug=False) m.ensure_default_constraints() - # cstr = "noise|white|variance" - # m.unconstrain(cstr); m.constrain_bounded(cstr, 1e-10, 1.) + for i, Y in enumerate(Ylist): + m.set('{}_noise'.format(i + 1), Y.var() / 100.) + +# import ipdb;ipdb.set_trace() + cstr = "variance" + m.unconstrain(cstr); m.constrain_bounded(cstr, 1e-15, 1.) + +# print "initializing beta" +# cstr = "noise" +# m.unconstrain(cstr); m.constrain_fixed(cstr) +# m.optimize('scg', messages=1, max_f_eval=200) +# +# print "releasing beta" +# cstr = "noise" +# m.unconstrain(cstr); m.constrain_positive(cstr) + m.auto_scale_factor = True diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 81e46774..943db420 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -5,11 +5,12 @@ Created on 10 Apr 2013 ''' from GPy.core import model from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM -import numpy from GPy.models.sparse_GP import sparse_GP +from GPy.util.linalg import PCA +from scipy import linalg +import numpy import itertools import pylab -from GPy.util.linalg import PCA class MRD(model): """ @@ -23,8 +24,10 @@ class MRD(model): :type names: [str] :param Q: latent dimensionality (will raise :type Q: int - :param init: initialisation method for the latent space - :type init: 'PCA'|'random' + :param initx: initialisation method for the latent space + :type initx: 'PCA'|'random' + :param initz: initialisation method for inducing inputs + :type initz: 'permute'|'random' :param X: Initial latent space :param X_variance: @@ -38,6 +41,7 @@ class MRD(model): :param kernel: kernel to use """ + def __init__(self, *Ylist, **kwargs): if kwargs.has_key("_debug"): self._debug = kwargs['_debug'] @@ -55,24 +59,30 @@ class MRD(model): del kwargs['kernel'] else: k = lambda: None - if kwargs.has_key('init'): - init = kwargs['init'] - del kwargs['init'] + if kwargs.has_key('initx'): + initx = kwargs['initx'] + del kwargs['initx'] else: - init = "PCA" + initx = "PCA" + if kwargs.has_key('initz'): + initz = kwargs['initz'] + del kwargs['initz'] + else: + initz = "permute" try: self.Q = kwargs["Q"] except KeyError: raise ValueError("Need Q for MRD") try: self.M = kwargs["M"] + del kwargs["M"] except KeyError: self.M = 10 self._init = True - X = self._init_X(init, Ylist) - Z = numpy.random.permutation(X.copy())[:self.M] - self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, **kwargs) for Y in Ylist] + X = self._init_X(initx, Ylist) + Z = self._init_Z(initz, X) + self.bgplvms = [Bayesian_GPLVM(Y, kernel=k(), X=X, Z=Z, M=self.M, **kwargs) for Y in Ylist] del self._init self.gref = self.bgplvms[0] @@ -96,6 +106,26 @@ class MRD(model): if not self._init: raise AttributeError("bgplvm list not initialized") @property + def Z(self): + return self.gref.Z + @Z.setter + def Z(self, Z): + try: + self.propagate_param(Z=Z) + except AttributeError: + if not self._init: + raise AttributeError("bgplvm list not initialized") + @property + def X_variance(self): + return self.gref.X_variance + @X_variance.setter + def X_variance(self, X_var): + try: + self.propagate_param(X_variance=X_var) + except AttributeError: + if not self._init: + raise AttributeError("bgplvm list not initialized") + @property def Ylist(self): return [g.likelihood.Y for g in self.bgplvms] @Ylist.setter @@ -120,6 +150,11 @@ class MRD(model): for g in self.bgplvms: g.__setattr__(key, val) + def randomize(self, initx='concat', initz='permute', *args, **kw): + super(MRD, self).randomize(*args, **kw) + self._init_X(initx, self.Ylist) + self._init_Z(initz, self.X) + def _get_param_names(self): # X_names = sum([['X_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], []) # S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.Q)] for n in range(self.N)], []) @@ -154,14 +189,14 @@ class MRD(model): params = numpy.hstack([X, X_var, Z, numpy.hstack(thetas)]) return params - def _set_var_params(self, g, X, X_var, Z): - g.X = X - g.X_variance = X_var - g.Z = Z - - def _set_kern_params(self, g, p): - g.kern._set_params(p[:g.kern.Nparam]) - g.likelihood._set_params(p[g.kern.Nparam:]) +# def _set_var_params(self, g, X, X_var, Z): +# g.X = X.reshape(self.N, self.Q) +# g.X_variance = X_var.reshape(self.N, self.Q) +# g.Z = Z.reshape(self.M, self.Q) +# +# def _set_kern_params(self, g, p): +# g.kern._set_params(p[:g.kern.Nparam]) +# g.likelihood._set_params(p[g.kern.Nparam:]) def _set_params(self, x): start = 0; end = self.NQ @@ -225,25 +260,50 @@ class MRD(model): self.X = X return X + + def _init_Z(self, init="permute", X=None): + if X is None: + X = self.X + if init in "permute": + Z = numpy.random.permutation(X.copy())[:self.M] + elif init in "random": + Z = numpy.random.randn(self.M, self.Q) * X.var() + self.Z = Z + return Z + def plot_X_1d(self, colors=None): - if colors is None: - colors = pylab.gca()._get_lines.color_cycle - fig = pylab.figure(num="MRD X 1d", figsize=(4 * len(self.bgplvms), (2 * self.X.shape[1]))) + fig = pylab.figure(num="MRD X 1d", figsize=(min(8, (3 * len(self.bgplvms))), min(12, (2 * self.X.shape[1])))) fig.clf() ax1 = fig.add_subplot(self.X.shape[1], 1, 1) + if colors is None: + colors = ax1._get_lines.color_cycle ax1.plot(self.X, c='k', alpha=.3) plots = ax1.plot(self.X.T[0], c=colors.next()) + ax1.fill_between(numpy.arange(self.X.shape[0]), + self.X.T[0] - 2 * numpy.sqrt(self.gref.X_variance.T[0]), + self.X.T[0] + 2 * numpy.sqrt(self.gref.X_variance.T[0]), + facecolor=plots[-1].get_color(), + alpha=.3) + ax1.text(1, 1, r"$\mathbf{{X_{}}}".format(1), + horizontalalignment='right', + verticalalignment='top', + transform=ax1.transAxes) for i in range(self.X.shape[1] - 1): ax = fig.add_subplot(self.X.shape[1], 1, i + 2) ax.plot(self.X, c='k', alpha=.3) plots.extend(ax.plot(self.X.T[i + 1], c=colors.next())) + ax.fill_between(numpy.arange(self.X.shape[0]), + self.X.T[i + 1] - 2 * numpy.sqrt(self.gref.X_variance.T[i + 1]), + self.X.T[i + 1] + 2 * numpy.sqrt(self.gref.X_variance.T[i + 1]), + facecolor=plots[-1].get_color(), + alpha=.3) if i < self.X.shape[1] - 2: ax.set_xticklabels('') ax1.set_xticklabels('') - ax1.legend(plots, [r"$\mathbf{{X_{}}}$".format(i + 1) for i in range(self.X.shape[1])], - bbox_to_anchor=(0., 1 + .01 * self.X.shape[1], - 1., 1. + .01 * self.X.shape[1]), loc=3, - ncol=self.X.shape[1], mode="expand", borderaxespad=0.) +# ax1.legend(plots, [r"$\mathbf{{X_{}}}$".format(i + 1) for i in range(self.X.shape[1])], +# bbox_to_anchor=(0., 1 + .01 * self.X.shape[1], +# 1., 1. + .01 * self.X.shape[1]), loc=3, +# ncol=self.X.shape[1], mode="expand", borderaxespad=0.) pylab.draw() fig.tight_layout(h_pad=.01, rect=(0, 0, 1, .95)) return fig