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

GPy/inference/latent_function_inference/__init__.py

@@ -64,8 +64,7 @@ class InferenceMethodList(LatentFunctionInference, list):
 from .exact_gaussian_inference import ExactGaussianInference
 from .laplace import Laplace,LaplaceBlock
 from GPy.inference.latent_function_inference.var_dtc import VarDTC
-from .expectation_propagation import EP
-from .expectation_propagation_dtc import EPDTC
+from .expectation_propagation import EP, EPDTC
 from .dtc import DTC
 from .fitc import FITC
 from .var_dtc_parallel import VarDTC_minibatch

GPy/inference/latent_function_inference/dtc.py

@@ -28,8 +28,8 @@ class DTC(LatentFunctionInference):
         num_data, output_dim = Y.shape
 
         #make sure the noise is not hetero
-        beta = 1./likelihood.gaussian_variance(Y_metadata)
-        if beta.size > 1:
+        precision = 1./likelihood.gaussian_variance(Y_metadata)
+        if precision.size > 1:
             raise NotImplementedError("no hetero noise with this implementation of DTC")
 
         Kmm = kern.K(Z)
@@ -42,7 +42,7 @@ class DTC(LatentFunctionInference):
         Kmmi, L, Li, _ = pdinv(Kmm)
 
         # Compute A
-        LiUTbeta = np.dot(Li, U.T)*np.sqrt(beta)
+        LiUTbeta = np.dot(Li, U.T)*np.sqrt(precision)
         A = tdot(LiUTbeta) + np.eye(num_inducing)
 
         # factor A
@@ -50,7 +50,7 @@ class DTC(LatentFunctionInference):
 
         # back substutue to get b, P, v
         tmp, _ = dtrtrs(L, Uy, lower=1)
-        b, _ = dtrtrs(LA, tmp*beta, lower=1)
+        b, _ = dtrtrs(LA, tmp*precision, lower=1)
         tmp, _ = dtrtrs(LA, b, lower=1, trans=1)
         v, _ = dtrtrs(L, tmp, lower=1, trans=1)
         tmp, _ = dtrtrs(LA, Li, lower=1, trans=0)
@@ -59,8 +59,8 @@ class DTC(LatentFunctionInference):
         #compute log marginal
         log_marginal = -0.5*num_data*output_dim*np.log(2*np.pi) + \
                        -np.sum(np.log(np.diag(LA)))*output_dim + \
-                       0.5*num_data*output_dim*np.log(beta) + \
-                       -0.5*beta*np.sum(np.square(Y)) + \
+                       0.5*num_data*output_dim*np.log(precision) + \
+                       -0.5*precision*np.sum(np.square(Y)) + \
                        0.5*np.sum(np.square(b))
 
         # Compute dL_dKmm
@@ -70,11 +70,11 @@ class DTC(LatentFunctionInference):
         # Compute dL_dU
         vY = np.dot(v.reshape(-1,1),Y.T)
         dL_dU = vY - np.dot(vvT_P, U.T)
-        dL_dU *= beta
+        dL_dU *= precision
 
         #compute dL_dR
         Uv = np.dot(U, v)
-        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1))*beta**2
+        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./precision + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1))*precision**2
 
         dL_dthetaL = likelihood.exact_inference_gradients(dL_dR)
 
@@ -97,8 +97,8 @@ class vDTC(object):
         num_data, output_dim = Y.shape
 
         #make sure the noise is not hetero
-        beta = 1./likelihood.gaussian_variance(Y_metadata)
-        if beta.size > 1:
+        precision = 1./likelihood.gaussian_variance(Y_metadata)
+        if precision.size > 1:
             raise NotImplementedError("no hetero noise with this implementation of DTC")
 
         Kmm = kern.K(Z)
@@ -111,9 +111,9 @@ class vDTC(object):
         Kmmi, L, Li, _ = pdinv(Kmm)
 
         # Compute A
-        LiUTbeta = np.dot(Li, U.T)*np.sqrt(beta)
+        LiUTbeta = np.dot(Li, U.T)*np.sqrt(precision)
         A_ = tdot(LiUTbeta)
-        trace_term = -0.5*(np.sum(Knn)*beta - np.trace(A_))
+        trace_term = -0.5*(np.sum(Knn)*precision - np.trace(A_))
         A = A_ + np.eye(num_inducing)
 
         # factor A
@@ -121,7 +121,7 @@ class vDTC(object):
 
         # back substutue to get b, P, v
         tmp, _ = dtrtrs(L, Uy, lower=1)
-        b, _ = dtrtrs(LA, tmp*beta, lower=1)
+        b, _ = dtrtrs(LA, tmp*precision, lower=1)
         tmp, _ = dtrtrs(LA, b, lower=1, trans=1)
         v, _ = dtrtrs(L, tmp, lower=1, trans=1)
         tmp, _ = dtrtrs(LA, Li, lower=1, trans=0)
@@ -131,8 +131,8 @@ class vDTC(object):
         #compute log marginal
         log_marginal = -0.5*num_data*output_dim*np.log(2*np.pi) + \
                        -np.sum(np.log(np.diag(LA)))*output_dim + \
-                       0.5*num_data*output_dim*np.log(beta) + \
-                       -0.5*beta*np.sum(np.square(Y)) + \
+                       0.5*num_data*output_dim*np.log(precision) + \
+                       -0.5*precision*np.sum(np.square(Y)) + \
                        0.5*np.sum(np.square(b)) + \
                        trace_term
 
@@ -145,15 +145,15 @@ class vDTC(object):
         vY = np.dot(v.reshape(-1,1),Y.T)
         #dL_dU = vY - np.dot(vvT_P, U.T)
         dL_dU = vY - np.dot(vvT_P - Kmmi, U.T)
-        dL_dU *= beta
+        dL_dU *= precision
 
         #compute dL_dR
         Uv = np.dot(U, v)
-        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) )*beta**2
-        dL_dR -=beta*trace_term/num_data
+        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./precision + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) )*precision**2
+        dL_dR -=precision*trace_term/num_data
 
         dL_dthetaL = likelihood.exact_inference_gradients(dL_dR)
-        grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn) + -0.5*beta, 'dL_dKnm':dL_dU.T, 'dL_dthetaL':dL_dthetaL}
+        grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn) + -0.5*precision, 'dL_dKnm':dL_dU.T, 'dL_dthetaL':dL_dthetaL}
 
         #construct a posterior object
         post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L)

GPy/inference/latent_function_inference/exact_gaussian_inference.py

@@ -22,21 +22,7 @@ class ExactGaussianInference(LatentFunctionInference):
     def __init__(self):
         pass#self._YYTfactor_cache = caching.cache()
 
-    def get_YYTfactor(self, Y):
-        """
-        find a matrix L which satisfies LL^T = YY^T.
-
-        Note that L may have fewer columns than Y, else L=Y.
-        """
-        N, D = Y.shape
-        if (N>D):
-            return Y
-        else:
-            #if Y in self.cache, return self.Cache[Y], else store Y in cache and return L.
-            #print "WARNING: N>D of Y, we need caching of L, such that L*L^T = Y, returning Y still!"
-            return Y
-
-    def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None):
+    def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None, K=None, precision=None):
         """
         Returns a Posterior class containing essential quantities of the posterior
         """
@@ -46,13 +32,17 @@ class ExactGaussianInference(LatentFunctionInference):
         else:
             m = mean_function.f(X)
 
+        if precision is None:
+            precision = likelihood.gaussian_variance(Y_metadata)
 
-        YYT_factor = self.get_YYTfactor(Y-m)
+        YYT_factor = Y-m
 
-        K = kern.K(X)
+        if K is None:
+            K = kern.K(X)
 
         Ky = K.copy()
-        diag.add(Ky, likelihood.gaussian_variance(Y_metadata)+1e-8)
+        diag.add(Ky, precision+1e-8)
+
         Wi, LW, LWi, W_logdet = pdinv(Ky)
 
         alpha, _ = dpotrs(LW, YYT_factor, lower=1)

GPy/inference/latent_function_inference/expectation_propagation.py

@@ -1,12 +1,14 @@
 # Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 import numpy as np
-from ...util.linalg import pdinv,jitchol,DSYR,tdot,dtrtrs, dpotrs
-from .posterior import Posterior
-from . import LatentFunctionInference
+from ...util.linalg import jitchol, DSYR, dtrtrs, dtrtri
+from ...core.parameterization.observable_array import ObsAr
+from . import ExactGaussianInference, VarDTC
+from ...util import diag
+
 log_2_pi = np.log(2*np.pi)
 
-class EP(LatentFunctionInference):
+class EPBase(object):
     def __init__(self, epsilon=1e-6, eta=1., delta=1.):
         """
         The expectation-propagation algorithm.
@@ -19,6 +21,7 @@ class EP(LatentFunctionInference):
         :param delta: damping EP updates factor.
         :type delta: float64
         """
+        super(EPBase, self).__init__()
         self.epsilon, self.eta, self.delta = epsilon, eta, delta
         self.reset()
 
@@ -33,32 +36,22 @@ class EP(LatentFunctionInference):
         # TODO: update approximation in the end as well? Maybe even with a switch?
         pass
 
-    def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None, Z=None):
-        assert mean_function is None, "inference with a mean function not implemented"
+class EP(EPBase, ExactGaussianInference):
+    def inference(self, kern, X, likelihood, Y, mean_function=None, Y_metadata=None, precision=None, K=None):
         num_data, output_dim = Y.shape
         assert output_dim ==1, "ep in 1D only (for now!)"
 
-        K = kern.K(X)
+        if K is None:
+            K = kern.K(X)
 
         if self._ep_approximation is None:
-
             #if we don't yet have the results of runnign EP, run EP and store the computed factors in self._ep_approximation
             mu, Sigma, mu_tilde, tau_tilde, Z_hat = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
         else:
             #if we've already run EP, just use the existing approximation stored in self._ep_approximation
             mu, Sigma, mu_tilde, tau_tilde, Z_hat = self._ep_approximation
 
-        Wi, LW, LWi, W_logdet = pdinv(K + np.diag(1./tau_tilde))
-
-        alpha, _ = dpotrs(LW, mu_tilde, lower=1)
-
-        log_marginal =  0.5*(-num_data * log_2_pi - W_logdet - np.sum(alpha * mu_tilde)) # TODO: add log Z_hat??
-
-        dL_dK = 0.5 * (tdot(alpha[:,None]) - Wi)
-
-        dL_dthetaL = np.zeros(likelihood.size)#TODO: derivatives of the likelihood parameters
-
-        return Posterior(woodbury_inv=Wi, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL}
+        return super(EP, self).inference(kern, X, likelihood, mu_tilde[:,None], mean_function=mean_function, Y_metadata=Y_metadata, precision=1./tau_tilde, K=K)
 
     def expectation_propagation(self, K, Y, likelihood, Y_metadata):
 
@@ -69,6 +62,7 @@ class EP(LatentFunctionInference):
         #Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
         mu = np.zeros(num_data)
         Sigma = K.copy()
+        diag.add(Sigma, 1e-7)
 
         #Initial values - Marginal moments
         Z_hat = np.empty(num_data,dtype=np.float64)
@@ -79,14 +73,14 @@ class EP(LatentFunctionInference):
         if self.old_mutilde is None:
             tau_tilde, mu_tilde, v_tilde = np.zeros((3, num_data))
         else:
-            assert old_mutilde.size == num_data, "data size mis-match: did you change the data? try resetting!"
+            assert self.old_mutilde.size == num_data, "data size mis-match: did you change the data? try resetting!"
             mu_tilde, v_tilde = self.old_mutilde, self.old_vtilde
             tau_tilde = v_tilde/mu_tilde
 
         #Approximation
         tau_diff = self.epsilon + 1.
         v_diff = self.epsilon + 1.
-       	iterations = 0
+        iterations = 0
         while (tau_diff > self.epsilon) or (v_diff > self.epsilon):
             update_order = np.random.permutation(num_data)
             for i in update_order:
@@ -124,3 +118,120 @@ class EP(LatentFunctionInference):
 
         mu_tilde = v_tilde/tau_tilde
         return mu, Sigma, mu_tilde, tau_tilde, Z_hat
+
+class EPDTC(EPBase, VarDTC):
+    def inference(self, kern, X, Z, likelihood, Y, mean_function=None, Y_metadata=None, Lm=None, dL_dKmm=None, psi0=None, psi1=None, psi2=None):
+        assert Y.shape[1]==1, "ep in 1D only (for now!)"
+
+        Kmm = kern.K(Z)
+        if psi1 is None:
+            try:
+                Kmn = kern.K(Z, X)
+            except TypeError:
+                Kmn = kern.psi1(Z, X).T
+        else:
+            Kmn = psi1.T
+
+        if self._ep_approximation is None:
+            mu, Sigma, mu_tilde, tau_tilde, Z_hat = self._ep_approximation = self.expectation_propagation(Kmm, Kmn, Y, likelihood, Y_metadata)
+        else:
+            mu, Sigma, mu_tilde, tau_tilde, Z_hat = self._ep_approximation
+
+        return super(EPDTC, self).inference(kern, X, Z, likelihood, mu_tilde,
+                                            mean_function=mean_function,
+                                            Y_metadata=Y_metadata,
+                                            precision=tau_tilde,
+                                            Lm=Lm, dL_dKmm=dL_dKmm,
+                                            psi0=psi0, psi1=psi1, psi2=psi2)
+
+    def expectation_propagation(self, Kmm, Kmn, Y, likelihood, Y_metadata):
+        num_data, output_dim = Y.shape
+        assert output_dim == 1, "This EP methods only works for 1D outputs"
+
+        LLT0 = Kmm.copy()
+        #diag.add(LLT0, 1e-8)
+
+        Lm = jitchol(LLT0)
+        Lmi = dtrtri(Lm)
+        Kmmi = np.dot(Lmi.T,Lmi)
+        KmmiKmn = np.dot(Kmmi,Kmn)
+        Qnn_diag = np.sum(Kmn*KmmiKmn,-2)
+
+        #Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
+        mu = np.zeros(num_data)
+        LLT = Kmm.copy() #Sigma = K.copy()
+        Sigma_diag = Qnn_diag.copy() + 1e-8
+
+        #Initial values - Marginal moments
+        Z_hat = np.zeros(num_data,dtype=np.float64)
+        mu_hat = np.zeros(num_data,dtype=np.float64)
+        sigma2_hat = np.zeros(num_data,dtype=np.float64)
+
+        #initial values - Gaussian factors
+        if self.old_mutilde is None:
+            tau_tilde, mu_tilde, v_tilde = np.zeros((3, num_data))
+        else:
+            assert self.old_mutilde.size == num_data, "data size mis-match: did you change the data? try resetting!"
+            mu_tilde, v_tilde = self.old_mutilde, self.old_vtilde
+            tau_tilde = v_tilde/mu_tilde
+
+        #Approximation
+        tau_diff = self.epsilon + 1.
+        v_diff = self.epsilon + 1.
+        iterations = 0
+        tau_tilde_old = 0.
+        v_tilde_old = 0.
+        update_order = np.random.permutation(num_data)
+
+        while (tau_diff > self.epsilon) or (v_diff > self.epsilon):
+            for i in update_order:
+                #Cavity distribution parameters
+                tau_cav = 1./Sigma_diag[i] - self.eta*tau_tilde[i]
+                v_cav = mu[i]/Sigma_diag[i] - self.eta*v_tilde[i]
+                #Marginal moments
+                Z_hat[i], mu_hat[i], sigma2_hat[i] = likelihood.moments_match_ep(Y[i], tau_cav, v_cav)#, Y_metadata=None)#=(None if Y_metadata is None else Y_metadata[i]))
+                #Site parameters update
+                delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
+                delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
+                tau_tilde[i] += delta_tau
+                v_tilde[i] += delta_v
+                #Posterior distribution parameters update
+
+                #DSYR(Sigma, Sigma[:,i].copy(), -delta_tau/(1.+ delta_tau*Sigma[i,i]))
+                DSYR(LLT,Kmn[:,i].copy(),delta_tau)
+                L = jitchol(LLT+np.eye(LLT.shape[0])*1e-7)
+
+                V,info = dtrtrs(L,Kmn,lower=1)
+                Sigma_diag = np.sum(V*V,-2)
+                si = np.sum(V.T*V[:,i],-1)
+                mu += (delta_v-delta_tau*mu[i])*si
+                #mu = np.dot(Sigma, v_tilde)
+
+            #(re) compute Sigma and mu using full Cholesky decompy
+            LLT = LLT0 + np.dot(Kmn*tau_tilde[None,:],Kmn.T)
+            #diag.add(LLT, 1e-8)
+            L = jitchol(LLT)
+            V, _ = dtrtrs(L,Kmn,lower=1)
+            V2, _ = dtrtrs(L.T,V,lower=0)
+            #Sigma_diag = np.sum(V*V,-2)
+            #Knmv_tilde = np.dot(Kmn,v_tilde)
+            #mu = np.dot(V2.T,Knmv_tilde)
+            Sigma = np.dot(V2.T,V2)
+            mu = np.dot(Sigma,v_tilde)
+
+            #monitor convergence
+            #if iterations>0:
+            tau_diff = np.mean(np.square(tau_tilde-tau_tilde_old))
+            v_diff = np.mean(np.square(v_tilde-v_tilde_old))
+
+            tau_tilde_old = tau_tilde.copy()
+            v_tilde_old = v_tilde.copy()
+
+            # Only to while loop once:?
+            tau_diff = 0
+            v_diff = 0
+            iterations += 1
+
+        mu_tilde = v_tilde/tau_tilde
+        return mu, Sigma, ObsAr(mu_tilde[:,None]), tau_tilde, Z_hat
+

GPy/inference/latent_function_inference/expectation_propagation_dtc.py

@@ -1,352 +0,0 @@
-# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-import numpy as np
-from ...util import diag
-from ...util.linalg import mdot, jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify, DSYR
-from ...core.parameterization.variational import VariationalPosterior
-from . import LatentFunctionInference
-from .posterior import Posterior
-log_2_pi = np.log(2*np.pi)
-
-class EPDTC(LatentFunctionInference):
-    const_jitter = 1e-6
-    def __init__(self, epsilon=1e-6, eta=1., delta=1., limit=1):
-        from ...util.caching import Cacher
-        self.limit = limit
-        self.get_trYYT = Cacher(self._get_trYYT, limit)
-        self.get_YYTfactor = Cacher(self._get_YYTfactor, limit)
-
-        self.epsilon, self.eta, self.delta = epsilon, eta, delta
-        self.reset()
-
-    def set_limit(self, limit):
-        self.get_trYYT.limit = limit
-        self.get_YYTfactor.limit = limit
-
-    def on_optimization_start(self):
-        self._ep_approximation = None
-
-    def on_optimization_end(self):
-        # TODO: update approximation in the end as well? Maybe even with a switch?
-        pass
-
-    def _get_trYYT(self, Y):
-        return np.sum(np.square(Y))
-
-    def __getstate__(self):
-        # has to be overridden, as Cacher objects cannot be pickled.
-        return self.limit
-
-    def __setstate__(self, state):
-        # has to be overridden, as Cacher objects cannot be pickled.
-        self.limit = state
-        from ...util.caching import Cacher
-        self.get_trYYT = Cacher(self._get_trYYT, self.limit)
-        self.get_YYTfactor = Cacher(self._get_YYTfactor, self.limit)
-
-    def _get_YYTfactor(self, Y):
-        """
-        find a matrix L which satisfies LLT = YYT.
-
-        Note that L may have fewer columns than Y.
-        """
-        N, D = Y.shape
-        if (N>=D):
-            return Y
-        else:
-            return jitchol(tdot(Y))
-
-    def get_VVTfactor(self, Y, prec):
-        return Y * prec # TODO chache this, and make it effective
-
-    def reset(self):
-        self.old_mutilde, self.old_vtilde = None, None
-        self._ep_approximation = None
-
-    def inference(self, kern, X, Z, likelihood, Y, mean_function=None, Y_metadata=None):
-        assert mean_function is None, "inference with a mean function not implemented"
-        num_data, output_dim = Y.shape
-        assert output_dim ==1, "ep in 1D only (for now!)"
-
-        Kmm = kern.K(Z)
-        Kmn = kern.K(Z,X)
-
-        if self._ep_approximation is None:
-            mu, Sigma, mu_tilde, tau_tilde, Z_hat = self._ep_approximation = self.expectation_propagation(Kmm, Kmn, Y, likelihood, Y_metadata)
-        else:
-            mu, Sigma, mu_tilde, tau_tilde, Z_hat = self._ep_approximation
-
-
-        if isinstance(X, VariationalPosterior):
-            uncertain_inputs = True
-            psi0 = kern.psi0(Z, X)
-            psi1 = Kmn.T#kern.psi1(Z, X)
-            psi2 = kern.psi2(Z, X)
-        else:
-            uncertain_inputs = False
-            psi0 = kern.Kdiag(X)
-            psi1 = Kmn.T#kern.K(X, Z)
-            psi2 = None
-
-        #see whether we're using variational uncertain inputs
-
-        _, output_dim = Y.shape
-
-        #see whether we've got a different noise variance for each datum
-        #beta = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6)
-        beta = tau_tilde
-        VVT_factor = beta[:,None]*mu_tilde[:,None]
-        trYYT = self.get_trYYT(mu_tilde[:,None])
-
-        # do the inference:
-        het_noise = beta.size > 1
-        num_inducing = Z.shape[0]
-        num_data = Y.shape[0]
-        # kernel computations, using BGPLVM notation
-
-        Kmm = kern.K(Z).copy()
-        diag.add(Kmm, self.const_jitter)
-        Lm = jitchol(Kmm)
-
-        # The rather complex computations of A
-        if uncertain_inputs:
-            if het_noise:
-                psi2_beta = psi2 * (beta.flatten().reshape(num_data, 1, 1)).sum(0)
-            else:
-                psi2_beta = psi2.sum(0) * beta
-            LmInv = dtrtri(Lm)
-            A = LmInv.dot(psi2_beta.dot(LmInv.T))
-        else:
-            if het_noise:
-                tmp = psi1 * (np.sqrt(beta.reshape(num_data, 1)))
-            else:
-                tmp = psi1 * (np.sqrt(beta))
-            tmp, _ = dtrtrs(Lm, tmp.T, lower=1)
-            A = tdot(tmp) #print A.sum()
-
-        # factor B
-        B = np.eye(num_inducing) + A
-        LB = jitchol(B)
-        psi1Vf = np.dot(psi1.T, VVT_factor)
-        # back substutue C into psi1Vf
-        tmp, _ = dtrtrs(Lm, psi1Vf, lower=1, trans=0)
-        _LBi_Lmi_psi1Vf, _ = dtrtrs(LB, tmp, lower=1, trans=0)
-        tmp, _ = dtrtrs(LB, _LBi_Lmi_psi1Vf, lower=1, trans=1)
-        Cpsi1Vf, _ = dtrtrs(Lm, tmp, lower=1, trans=1)
-
-        # data fit and derivative of L w.r.t. Kmm
-        delit = tdot(_LBi_Lmi_psi1Vf)
-        data_fit = np.trace(delit)
-        DBi_plus_BiPBi = backsub_both_sides(LB, output_dim * np.eye(num_inducing) + delit)
-        delit = -0.5 * DBi_plus_BiPBi
-        delit += -0.5 * B * output_dim
-        delit += output_dim * np.eye(num_inducing)
-        # Compute dL_dKmm
-        dL_dKmm = backsub_both_sides(Lm, delit)
-
-        # derivatives of L w.r.t. psi
-        dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm,
-            VVT_factor, Cpsi1Vf, DBi_plus_BiPBi,
-            psi1, het_noise, uncertain_inputs)
-
-        # log marginal likelihood
-        log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise,
-            psi0, A, LB, trYYT, data_fit, VVT_factor)
-
-        #put the gradients in the right places
-        dL_dR = _compute_dL_dR(likelihood,
-            het_noise, uncertain_inputs, LB,
-            _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A,
-            psi0, psi1, beta,
-            data_fit, num_data, output_dim, trYYT, mu_tilde[:,None])
-
-        dL_dthetaL = 0#likelihood.exact_inference_gradients(dL_dR,Y_metadata)
-
-        if uncertain_inputs:
-            grad_dict = {'dL_dKmm': dL_dKmm,
-                         'dL_dpsi0':dL_dpsi0,
-                         'dL_dpsi1':dL_dpsi1,
-                         'dL_dpsi2':dL_dpsi2,
-                         'dL_dthetaL':dL_dthetaL}
-        else:
-            grad_dict = {'dL_dKmm': dL_dKmm,
-                         'dL_dKdiag':dL_dpsi0,
-                         'dL_dKnm':dL_dpsi1,
-                         'dL_dthetaL':dL_dthetaL}
-
-        #get sufficient things for posterior prediction
-        #TODO: do we really want to do this in  the loop?
-        if VVT_factor.shape[1] == Y.shape[1]:
-            woodbury_vector = Cpsi1Vf # == Cpsi1V
-        else:
-            print('foobar')
-            psi1V = np.dot(mu_tilde[:,None].T*beta, psi1).T
-            tmp, _ = dtrtrs(Lm, psi1V, lower=1, trans=0)
-            tmp, _ = dpotrs(LB, tmp, lower=1)
-            woodbury_vector, _ = dtrtrs(Lm, tmp, lower=1, trans=1)
-        Bi, _ = dpotri(LB, lower=1)
-        symmetrify(Bi)
-        Bi = -dpotri(LB, lower=1)[0]
-        diag.add(Bi, 1)
-
-        woodbury_inv = backsub_both_sides(Lm, Bi)
-
-        #construct a posterior object
-        post = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, K=Kmm, mean=None, cov=None, K_chol=Lm)
-        return post, log_marginal, grad_dict
-
-
-
-
-
-    def expectation_propagation(self, Kmm, Kmn, Y, likelihood, Y_metadata):
-
-        num_data, data_dim = Y.shape
-        assert data_dim == 1, "This EP methods only works for 1D outputs"
-
-        KmnKnm = np.dot(Kmn,Kmn.T)
-        Lm = jitchol(Kmm)
-        Lmi = dtrtrs(Lm,np.eye(Lm.shape[0]))[0] #chol_inv(Lm)
-        Kmmi = np.dot(Lmi.T,Lmi)
-        KmmiKmn = np.dot(Kmmi,Kmn)
-        Qnn_diag = np.sum(Kmn*KmmiKmn,-2)
-        LLT0 = Kmm.copy()
-
-        #Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
-        mu = np.zeros(num_data)
-        LLT = Kmm.copy() #Sigma = K.copy()
-        Sigma_diag = Qnn_diag.copy()
-
-        #Initial values - Marginal moments
-        Z_hat = np.empty(num_data,dtype=np.float64)
-        mu_hat = np.empty(num_data,dtype=np.float64)
-        sigma2_hat = np.empty(num_data,dtype=np.float64)
-
-        #initial values - Gaussian factors
-        if self.old_mutilde is None:
-            tau_tilde, mu_tilde, v_tilde = np.zeros((3, num_data))
-        else:
-            assert old_mutilde.size == num_data, "data size mis-match: did you change the data? try resetting!"
-            mu_tilde, v_tilde = self.old_mutilde, self.old_vtilde
-            tau_tilde = v_tilde/mu_tilde
-
-        #Approximation
-        tau_diff = self.epsilon + 1.
-        v_diff = self.epsilon + 1.
-       	iterations = 0
-        while (tau_diff > self.epsilon) or (v_diff > self.epsilon):
-            update_order = np.random.permutation(num_data)
-            for i in update_order:
-                #Cavity distribution parameters
-                tau_cav = 1./Sigma_diag[i] - self.eta*tau_tilde[i]
-                v_cav = mu[i]/Sigma_diag[i] - self.eta*v_tilde[i]
-                #Marginal moments
-                Z_hat[i], mu_hat[i], sigma2_hat[i] = likelihood.moments_match_ep(Y[i], tau_cav, v_cav)#, Y_metadata=None)#=(None if Y_metadata is None else Y_metadata[i]))
-                #Site parameters update
-                delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
-                delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
-                tau_tilde[i] += delta_tau
-                v_tilde[i] += delta_v
-                #Posterior distribution parameters update
-
-                #DSYR(Sigma, Sigma[:,i].copy(), -delta_tau/(1.+ delta_tau*Sigma[i,i]))
-                DSYR(LLT,Kmn[:,i].copy(),delta_tau)
-                L = jitchol(LLT)
-
-                V,info = dtrtrs(L,Kmn,lower=1)
-                Sigma_diag = np.sum(V*V,-2)
-                si = np.sum(V.T*V[:,i],-1)
-                mu += (delta_v-delta_tau*mu[i])*si
-                #mu = np.dot(Sigma, v_tilde)
-
-            #(re) compute Sigma and mu using full Cholesky decompy
-            LLT = LLT0 + np.dot(Kmn*tau_tilde[None,:],Kmn.T)
-            L = jitchol(LLT)
-            V,info = dtrtrs(L,Kmn,lower=1)
-            V2,info = dtrtrs(L.T,V,lower=0)
-            #Sigma_diag = np.sum(V*V,-2)
-            #Knmv_tilde = np.dot(Kmn,v_tilde)
-            #mu = np.dot(V2.T,Knmv_tilde)
-            Sigma = np.dot(V2.T,V2)
-            mu = np.dot(Sigma,v_tilde)
-
-            #monitor convergence
-            if iterations>0:
-                tau_diff = np.mean(np.square(tau_tilde-tau_tilde_old))
-                v_diff = np.mean(np.square(v_tilde-v_tilde_old))
-            tau_tilde_old = tau_tilde.copy()
-            v_tilde_old = v_tilde.copy()
-
-            tau_diff = 0
-            v_diff = 0
-            iterations += 1
-
-        mu_tilde = v_tilde/tau_tilde
-        return mu, Sigma, mu_tilde, tau_tilde, Z_hat
-
-def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs):
-    dL_dpsi0 = -0.5 * output_dim * (beta[:,None] * np.ones([num_data, 1])).flatten()
-    dL_dpsi1 = np.dot(VVT_factor, Cpsi1Vf.T)
-    dL_dpsi2_beta = 0.5 * backsub_both_sides(Lm, output_dim * np.eye(num_inducing) - DBi_plus_BiPBi)
-    if het_noise:
-        if uncertain_inputs:
-            dL_dpsi2 = beta[:, None, None] * dL_dpsi2_beta[None, :, :]
-        else:
-            dL_dpsi1 += 2.*np.dot(dL_dpsi2_beta, (psi1 * beta.reshape(num_data, 1)).T).T
-            dL_dpsi2 = None
-    else:
-        dL_dpsi2 = beta * dL_dpsi2_beta
-        if uncertain_inputs:
-            # repeat for each of the N psi_2 matrices
-            dL_dpsi2 = np.repeat(dL_dpsi2[None, :, :], num_data, axis=0)
-        else:
-            # subsume back into psi1 (==Kmn)
-            dL_dpsi1 += 2.*np.dot(psi1, dL_dpsi2)
-            dL_dpsi2 = None
-
-    return dL_dpsi0, dL_dpsi1, dL_dpsi2
-
-
-def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, psi0, psi1, beta, data_fit, num_data, output_dim, trYYT, Y):
-    # the partial derivative vector for the likelihood
-    if likelihood.size == 0:
-        # save computation here.
-        dL_dR = None
-    elif het_noise:
-        if uncertain_inputs:
-            raise NotImplementedError("heteroscedatic derivates with uncertain inputs not implemented")
-        else:
-            #from ...util.linalg import chol_inv
-            #LBi = chol_inv(LB)
-            LBi, _ = dtrtrs(LB,np.eye(LB.shape[0]))
-
-            Lmi_psi1, nil = dtrtrs(Lm, psi1.T, lower=1, trans=0)
-            _LBi_Lmi_psi1, _ = dtrtrs(LB, Lmi_psi1, lower=1, trans=0)
-
-            dL_dR = -0.5 * beta + 0.5 * (beta*Y)**2
-            dL_dR += 0.5 * output_dim * (psi0 - np.sum(Lmi_psi1**2,0))[:,None] * beta**2
-
-            dL_dR += 0.5*np.sum(mdot(LBi.T,LBi,Lmi_psi1)*Lmi_psi1,0)[:,None]*beta**2
-
-            dL_dR += -np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T * Y * beta**2
-            dL_dR += 0.5*np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T**2 * beta**2
-    else:
-        # likelihood is not heteroscedatic
-        dL_dR = -0.5 * num_data * output_dim * beta + 0.5 * trYYT * beta ** 2
-        dL_dR += 0.5 * output_dim * (psi0.sum() * beta ** 2 - np.trace(A) * beta)
-        dL_dR += beta * (0.5 * np.sum(A * DBi_plus_BiPBi) - data_fit)
-    return dL_dR
-
-def _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit,Y):
-    #compute log marginal likelihood
-    if het_noise:
-        lik_1 = -0.5 * num_data * output_dim * np.log(2. * np.pi) + 0.5 * np.sum(np.log(beta)) - 0.5 * np.sum(beta * np.square(Y).sum(axis=-1))
-        lik_2 = -0.5 * output_dim * (np.sum(beta.flatten() * psi0) - np.trace(A))
-    else:
-        lik_1 = -0.5 * num_data * output_dim * (np.log(2. * np.pi) - np.log(beta)) - 0.5 * beta * trYYT
-        lik_2 = -0.5 * output_dim * (np.sum(beta * psi0) - np.trace(A))
-    lik_3 = -output_dim * (np.sum(np.log(np.diag(LB))))
-    lik_4 = 0.5 * data_fit
-    log_marginal = lik_1 + lik_2 + lik_3 + lik_4
-    return log_marginal

GPy/inference/latent_function_inference/var_dtc.py

@@ -64,31 +64,30 @@ class VarDTC(LatentFunctionInference):
     def get_VVTfactor(self, Y, prec):
         return Y * prec # TODO chache this, and make it effective
 
-    def inference(self, kern, X, Z, likelihood, Y, Y_metadata=None, Lm=None, dL_dKmm=None, psi0=None, psi1=None, psi2=None):
+    def inference(self, kern, X, Z, likelihood, Y, Y_metadata=None, mean_function=None, precision=None, Lm=None, dL_dKmm=None, psi0=None, psi1=None, psi2=None):
+        assert mean_function is None, "inference with a mean function not implemented"
+
+        num_data, output_dim = Y.shape
+        num_inducing = Z.shape[0]
 
-        _, output_dim = Y.shape
         uncertain_inputs = isinstance(X, VariationalPosterior)
 
-        #see whether we've got a different noise variance for each datum
-        beta = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6)
-        # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
-        #self.YYTfactor = self.get_YYTfactor(Y)
-        #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta)
-        het_noise = beta.size > 1
-        if beta.ndim == 1:
-            beta = beta[:, None]
-        VVT_factor = beta*Y
-        #VVT_factor = beta*Y
+        if precision is None:
+            #assume Gaussian likelihood
+            precision = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), self.const_jitter)
+
+        if precision.ndim == 1:
+            precision = precision[:, None]
+        het_noise = precision.size > 1
+
+        VVT_factor = precision*Y
+        #VVT_factor = precision*Y
         trYYT = self.get_trYYT(Y)
 
-        # do the inference:
-        num_inducing = Z.shape[0]
-        num_data = Y.shape[0]
         # kernel computations, using BGPLVM notation
-
-        Kmm = kern.K(Z).copy()
-        diag.add(Kmm, self.const_jitter)
         if Lm is None:
+            Kmm = kern.K(Z).copy()
+            diag.add(Kmm, self.const_jitter)
             Lm = jitchol(Kmm)
 
         # The rather complex computations of A, and the psi stats
@@ -99,15 +98,16 @@ class VarDTC(LatentFunctionInference):
                 psi1 = kern.psi1(Z, X)
             if het_noise:
                 if psi2 is None:
-                    assert len(psi2.shape) == 3  # Need to have not summed out N
-                    #FIXME: Need testing
-                    psi2_beta = np.sum([psi2[X[i:i+1,:], :, :] * beta_i for i,beta_i in enumerate(beta)],0)
+                    psi2_beta = (kern.psi2n(Z, X) * precision[:, :, None]).sum(0)
                 else:
-                    psi2_beta = np.sum([kern.psi2(Z,X[i:i+1,:]) * beta_i for i,beta_i in enumerate(beta)],0)
+                    psi2_beta = (psi2 * precision[:, :, None]).sum(0)
             else:
                 if psi2 is None:
-                    psi2 = kern.psi2(Z,X)
-                psi2_beta =  psi2 * beta
+                    psi2_beta = kern.psi2(Z,X) * precision
+                elif psi2.ndim == 3:
+                    psi2_beta = psi2.sum(0) * precision
+                else:
+                    psi2_beta = psi2 * precision
             LmInv = dtrtri(Lm)
             A = LmInv.dot(psi2_beta.dot(LmInv.T))
         else:
@@ -116,9 +116,9 @@ class VarDTC(LatentFunctionInference):
             if psi1 is None:
                 psi1 = kern.K(X, Z)
             if het_noise:
-                tmp = psi1 * (np.sqrt(beta))
+                tmp = psi1 * (np.sqrt(precision))
             else:
-                tmp = psi1 * (np.sqrt(beta))
+                tmp = psi1 * (np.sqrt(precision))
             tmp, _ = dtrtrs(Lm, tmp.T, lower=1)
             A = tdot(tmp) #print A.sum()
 
@@ -144,19 +144,19 @@ class VarDTC(LatentFunctionInference):
             dL_dKmm = backsub_both_sides(Lm, delit)
 
         # derivatives of L w.r.t. psi
-        dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm,
+        dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, precision, Lm,
             VVT_factor, Cpsi1Vf, DBi_plus_BiPBi,
             psi1, het_noise, uncertain_inputs)
 
         # log marginal likelihood
-        log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise,
+        log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, precision, het_noise,
             psi0, A, LB, trYYT, data_fit, Y)
 
         #noise derivatives
         dL_dR = _compute_dL_dR(likelihood,
             het_noise, uncertain_inputs, LB,
             _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A,
-            psi0, psi1, beta,
+            psi0, psi1, precision,
             data_fit, num_data, output_dim, trYYT, Y, VVT_factor)
 
         dL_dthetaL = likelihood.exact_inference_gradients(dL_dR,Y_metadata)
@@ -181,7 +181,7 @@ class VarDTC(LatentFunctionInference):
         else:
             print('foobar')
             import ipdb; ipdb.set_trace()
-            psi1V = np.dot(Y.T*beta, psi1).T
+            psi1V = np.dot(Y.T*precision, psi1).T
             tmp, _ = dtrtrs(Lm, psi1V, lower=1, trans=0)
             tmp, _ = dpotrs(LB, tmp, lower=1)
             woodbury_vector, _ = dtrtrs(Lm, tmp, lower=1, trans=1)

GPy/inference/optimization/optimization.py

@@ -228,13 +228,35 @@ class opt_SCG(Optimizer):
         self.f_opt = self.trace[-1]
         self.funct_eval = opt_result[2]
         self.status = opt_result[3]
+        
+class Opt_Adadelta(Optimizer):
+    def __init__(self, step_rate=0.1, decay=0.9, momentum=0, *args, **kwargs):
+        Optimizer.__init__(self, *args, **kwargs)
+        self.opt_name = "Adadelta (climin)"
+        self.step_rate=step_rate
+        self.decay = decay
+        self.momentum = momentum
+
+    def opt(self, f_fp=None, f=None, fp=None):
+        assert not fp is None
+        
+        import climin
+                    
+        opt = climin.adadelta.Adadelta(self.x_init, fp, step_rate=self.step_rate, decay=self.decay, momentum=self.momentum)
+        
+        for info in opt:
+            if info['n_iter']>=self.max_iters:
+                self.x_opt =  opt.wrt
+                self.status = 'maximum number of function evaluations exceeded '
+                break
 
 def get_optimizer(f_min):
 
     optimizers = {'fmin_tnc': opt_tnc,
           'simplex': opt_simplex,
           'lbfgsb': opt_lbfgsb,
-          'scg': opt_SCG}
+          'scg': opt_SCG,
+          'adadelta':Opt_Adadelta}
 
     if rasm_available:
         optimizers['rasmussen'] = opt_rasm