Merge branch 'devel' of github.com:SheffieldML/GPy into devel

GPy/core/transformations.py

@@ -39,8 +39,8 @@ class logexp(transformation):
         return '(+ve)'
 
 class logexp_clipped(transformation):
-    max_bound = 1e300
+    max_bound = 1e250
-    min_bound = 1e-10
+    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):

GPy/examples/classification.py

@@ -9,8 +9,8 @@ import pylab as pb
 import numpy as np
 import GPy
 
-default_seed=10000
+default_seed = 10000
-def crescent_data(seed=default_seed): #FIXME
+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 = GPy.models.sparse_GP(data['X'], likelihood=likelihood, kernel=kernel, Z=Z, normalize_X=False)
-    m.set('len',2.)
+    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)
+    sample = np.random.randint(0, data['X'].shape[0], inducing)
-    Z = data['X'][sample,:]
+    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()
 

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)
@@ -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):
+def mrd_simulation(optimize=True, plot_sim=False):
-    D1, D2, D3, N, M, Q = 150, 250, 300, 700, 3, 7
+    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():

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.
+    d = -gradnew # Initial search direction.
-    success = True  # Force calculation of directional derivs.
+    success = True # Force calculation of directional derivs.
-    nsuccess = 0  # nsuccess counts number of successes.
+    nsuccess = 0 # nsuccess counts number of successes.
-    beta = 1.0  # Initial scale parameter.
+    beta = 1.0 # Initial scale parameter.
-    betamin = 1.0e-15  # Lower bound on scale.
+    betamin = 1.0e-15 # Lower bound on scale.
-    betamax = 1.0e100  # Upper bound on scale.
+    betamax = 1.0e100 # Upper bound on scale.
     status = "Not converged"
 
     flog = [fold]
@@ -107,12 +110,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()
 

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,53 @@ 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)/10)
+        elif self.learning_rate_adaptation == 'semi_pesky':
+            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
+
+            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()
-        self.model.grads = np.zeros_like(self.x_opt)
+        self.grads = []
 
         X, Y = self.model.X.copy(), self.model.likelihood.Y.copy()
 
@@ -235,6 +287,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 +325,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
@@ -303,6 +357,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()

GPy/kern/bias.py

@@ -55,7 +55,8 @@ 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

GPy/kern/kern.py

@@ -315,31 +315,20 @@ 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':
+            # TODO psi1 this must be faster/better/precached/more nice
-                pass
+            tmp1 = np.zeros((mu.shape[0], Z.shape[0]))
-            # rbf X bias
+            p1.psi1(Z, mu, S, tmp1)
-            elif p1.name == 'bias' and p2.name == 'rbf':
+            tmp2 = np.zeros((mu.shape[0], Z.shape[0]))
-                target += p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :])
+            p2.psi1(Z, mu, S, tmp2)
-            elif p2.name == 'bias' and p1.name == 'rbf':
+
-                target += p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :])
+            prod = np.multiply(tmp1, tmp2)
-            # linear X bias
+            crossterms += prod[:,:,None] + prod[:, None, :]
-            elif p1.name == 'bias' and p2.name == 'linear':
+            
-                tmp = np.zeros((mu.shape[0], Z.shape[0]))
+        target += crossterms
-                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):
@@ -348,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
+            tmp = np.zeros((mu.shape[0], Z.shape[0]))
-            if p1.name == 'white' or p2.name == 'white':
+            p1.psi1(Z, mu, S, tmp)
-                pass
+            p2.dpsi1_dtheta((tmp[:,None,:]*dL_dpsi2).sum(1)*2., Z, mu, S, target[ps2])
-            # 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"
 
         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):
+        for p1, p2 in itertools.permutations(self.parts, 2):
-            # white doesn;t combine with anything
+            if p1.name == 'linear' and p2.name == 'linear':
-            if p1.name == 'white' or p2.name == 'white':
+                raise NotImplementedError("We don't handle linear/linear cross-terms")
-                pass
+            tmp = np.zeros((mu.shape[0], Z.shape[0]))
-            # rbf X bias
+            p1.psi1(Z, mu, S, tmp)
-            elif p1.name == 'bias' and p2.name == 'rbf':
+            tmp2 = np.zeros_like(target)
-                p2.dpsi1_dX(dL_dpsi2.sum(1).T * p1.variance, Z, mu, S, target)
+            p2.dpsi1_dZ((tmp[:,None,:]*dL_dpsi2).sum(1).T, Z, mu, S, tmp2)
-            elif p2.name == 'bias' and p1.name == 'rbf':
+            target += tmp2
-                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"
 
-        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]))
@@ -420,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):
+        for p1, p2 in itertools.permutations(self.parts, 2):
-            # white doesn;t combine with anything
+            if p1.name == 'linear' and p2.name == 'linear':
-            if p1.name == 'white' or p2.name == 'white':
+                raise NotImplementedError("We don't handle linear/linear cross-terms")
-                pass
+
-            # rbf X bias
+            tmp = np.zeros((mu.shape[0], Z.shape[0]))
-            elif p1.name == 'bias' and p2.name == 'rbf':
+            p1.psi1(Z, mu, S, tmp)
-                p2.dpsi1_dmuS(dL_dpsi2.sum(1).T * p1.variance * 2., Z, mu, S, target_mu, target_S)
+            p2.dpsi1_dmuS((tmp[:,None,:]*dL_dpsi2).sum(1).T*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"
 
         return target_mu, target_S
 

GPy/kern/kernpart.py

@@ -54,5 +54,3 @@ class kernpart(object):
         raise NotImplementedError
     def dK_dX(self,X,X2,target):
         raise NotImplementedError
-
-

GPy/kern/white.py

@@ -18,6 +18,7 @@ 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
-

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:

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)
+    :type likelihood: GPy.likelihood.(Gaussian | EP | Laplace)
-    :param kernel : the kernel/covariance function. See link kernels
+    :param kernel : the kernel (covariance function). See link kernels
-    :type kernel: a GPy kernel
+    :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,49 +61,29 @@ class sparse_GP(GP):
             self.psi2 = None
 
     def _computations(self):
-#         sf = self.scale_factor
-#         sf2 = sf ** 2
 
         # factor Kmm
         self.Lm = jitchol(self.Kmm)
 
         # The rather complex computations of self.A
-        if self.likelihood.is_heteroscedastic:
+        if self.has_uncertain_inputs:
-            assert self.likelihood.D == 1  # TODO: what if the likelihood is heterscedatic and there are multiple independent outputs?
+            if self.likelihood.is_heteroscedastic:
-            if self.has_uncertain_inputs:
+                psi2_beta = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0)
-#                 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)
             else:
-#                 tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)) / sf)
+                psi2_beta = self.psi2.sum(0) * self.likelihood.precision
-                tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)))
+            evals, evecs = linalg.eigh(psi2_beta)
-                tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
+            clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable
-                self.A = tdot(tmp)
+            tmp = evecs * np.sqrt(clipped_evals)
         else:
-            if self.has_uncertain_inputs:
+            if self.likelihood.is_heteroscedastic:
-#                 psi2_beta_scaled = (self.psi2 * (self.likelihood.precision / sf2)).sum(0)
+                tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)))
-                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)
             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)
+        tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
-                self.A = tdot(tmp)
+        self.A = tdot(tmp)
+
 
         # factor B
-#         self.B = np.eye(self.M) / sf2 + self.A
         self.B = np.eye(self.M) + self.A
         self.LB = jitchol(self.B)
 
@@ -121,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)
@@ -132,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.flatten()[:, None, None] * dL_dpsi2_beta[None, :, :]
@@ -158,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)))
 
@@ -166,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
@@ -185,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):
@@ -205,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:

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