Fixed merge conflict

GPy/inference/latent_function_inference/__init__.py

@@ -26,3 +26,4 @@ etc.
 from exact_gaussian_inference import ExactGaussianInference
 from laplace import LaplaceInference
 expectation_propagation = 'foo' # TODO
+from dtc import DTC

GPy/inference/latent_function_inference/dtc.py

@@ -0,0 +1,96 @@
+# Copyright (c) 2012, James Hensman
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+from posterior import Posterior
+from ...util.linalg import jitchol, tdot, dtrtrs, dpotri, pdinv
+import numpy as np
+log_2_pi = np.log(2*np.pi)
+
+class DTC(object):
+    """
+    An object for inference when the likelihood is Gaussian, but we want to do sparse inference.
+
+    The function self.inference returns a Posterior object, which summarizes
+    the posterior.
+
+    NB. It's not recommended to use this function! It's here for historical purposes. 
+
+    """
+    def __init__(self):
+        self.const_jitter = 1e-6
+
+    def inference(self, kern, X, X_variance, Z, likelihood, Y):
+        assert X_variance is None, "cannot use X_variance with DTC. Try varDTC."
+
+        num_inducing, _ = Z.shape
+        num_data, output_dim = Y.shape
+
+        #make sure the noise is not hetero
+        beta = 1./np.squeeze(likelihood.variance)
+        if beta.size <1:
+            raise NotImplementedError, "no hetero noise with this implementatino of DTC"
+
+        Kmm = kern.K(Z)
+        Knn = kern.Kdiag(X)
+        Knm = kern.K(X, Z)
+        U = Knm
+        Uy = np.dot(U.T,Y)
+
+        #factor Kmm 
+        Kmmi, L, Li, _ = pdinv(Kmm)
+
+        # Compute A
+        LiUT, _ = dtrtrs(L, U.T*np.sqrt(beta), lower=1)
+        A_I = tdot(LiUT)
+        A = A_I + np.eye(num_inducing)
+
+        # factor A
+        LA = jitchol(A)
+
+        # back substutue to get b, P, v
+        tmp, _ = dtrtrs(L, Uy, lower=1)
+        b, _ = dtrtrs(LA, tmp*beta, 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)
+        P = tdot(tmp.T)
+
+        #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*np.sum(np.square(b))
+
+        # Compute dL_dKmm
+        tmp, _ = dtrtrs(L, A_I, lower=1, trans=1)
+        dL_dK, _ = dtrtrs(L, tmp.T, lower=1, trans=0)
+        tmp, _ = dtrtrs(LA, tmp.T, lower=1, trans=1)
+        dL_dK -= tdot(tmp.T)
+        dL_dK *= output_dim
+        dL_dK -= tdot(v)
+        dL_dK /=2.
+
+        # Compute dL_dU
+        vvT_P = tdot(v.reshape(-1,1)) + P
+        vY = np.dot(v.reshape(-1,1),Y.T)
+        dL_dU = vY + np.dot(vvT_P, U.T)
+        dL_dU *= beta
+
+        #compute dL_dR
+        Uv = np.dot(U, v)
+        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - beta * np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1)
+                )*beta**2
+
+        grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn), 'dL_dKnm':dL_dU.T}
+
+        #update gradients
+        kern.update_gradients_sparse(X=X, Z=Z, **grad_dict)
+        likelihood.update_gradients(dL_dR)
+
+        #construct a posterior object
+        post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L)
+
+        return post, log_marginal, grad_dict
+
+

GPy/inference/latent_function_inference/varDTC.py

@@ -2,9 +2,8 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 from posterior import Posterior
-from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dpotri, symmetrify
+from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify
 import numpy as np
-from ...util.linalg import dtrtri
 from ...util.caching import Cacher
 from ...util.misc import param_to_array
 log_2_pi = np.log(2*np.pi)
@@ -85,7 +84,7 @@ class VarDTC(object):
                 tmp = tmp.T
             # no backsubstitution because of bound explosion on tr(A) if not...
             LmInv, _ = dtrtri(Lm, lower=1)
-            A = LmInv.T.dot(psi2_beta.dot(LmInv))
+            A = LmInv.dot(psi2_beta.dot(LmInv.T))
             #print A.sum()
         else:
             if het_noise:
@@ -97,6 +96,7 @@ class VarDTC(object):
 
         # factor B
         B = np.eye(num_inducing) + A
+        self.A = A
         LB = jitchol(B)
 
         # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!

GPy/inference/optimization/scg.py

@@ -69,8 +69,8 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=np.inf, display=True,
     success = True # Force calculation of directional derivs.
     nsuccess = 0 # nsuccess counts number of successes.
     beta = 1.0 # Initial scale parameter.
-    betamin = 1.0e-60 # Lower bound on scale.
-    betamax = 1.0e50 # Upper bound on scale.
+    betamin = 1.0e-15 # Lower bound on scale.
+    betamax = 1.0e15 # Upper bound on scale.
     status = "Not converged"
 
     flog = [fold]