diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 75820407..6615b40c 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -176,12 +176,13 @@ def bgplvm_simulation_matlab_compare(): Y = sim_data['Y'] S = sim_data['S'] mu = sim_data['mu'] - M, [_, Q] = 20, mu.shape + M, (_, Q) = 30, mu.shape from GPy.models import mrd from GPy import kern reload(mrd); reload(kern) - k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) + #k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) + k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, # X=mu, # X_variance=S, diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 413d8a07..14c789b8 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -76,8 +76,13 @@ class sparse_GP(GP): 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) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,self.psi2_beta_scaled.T,lower=1) - self.A, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1) + evals, evecs = linalg.eigh(self.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) else: tmp = self.psi1*(np.sqrt(self.likelihood.precision.flatten().reshape(1,self.N))/sf) self.psi2_beta_scaled = tdot(tmp) @@ -86,8 +91,14 @@ class sparse_GP(GP): else: if self.has_uncertain_inputs: self.psi2_beta_scaled = (self.psi2*(self.likelihood.precision/sf2)).sum(0) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,self.psi2_beta_scaled.T,lower=1) - self.A, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1) + evals, evecs = linalg.eigh(self.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) + self.psi2_beta_scaled = tdot(tmp) + 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) self.psi2_beta_scaled = tdot(tmp) diff --git a/GPy/util/datasets/BGPLVMSimulation.mat b/GPy/util/datasets/BGPLVMSimulation.mat new file mode 100644 index 00000000..c1cff0a0 Binary files /dev/null and b/GPy/util/datasets/BGPLVMSimulation.mat differ diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index b19aa2b6..42a98fea 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -72,7 +72,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 '\rWarning: adding jitter of {:.10e} '.format(jitter), + print 'Warning: adding jitter of {:.10e}'.format(jitter) try: return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True) except: