From 58b206a24541fe55183676d11b225116a61c9bb7 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Sun, 19 May 2013 19:16:25 +0100 Subject: [PATCH 1/3] enabled DSYR on feeble machines i.e. where numpy is compiled without proper blas linkage --- GPy/util/linalg.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index fa7de8c1..abfb1900 100644 --- a/GPy/util/linalg.py +++ b/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 From 3f22e61d2d09f17684cfc35a8b5d62ca9f273ad3 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Sun, 19 May 2013 19:17:55 +0100 Subject: [PATCH 2/3] tidied the computation of A in sparse_GP --- GPy/models/sparse_GP.py | 41 ++++++++++++----------------------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 2aafab16..ef7b445d 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -70,39 +70,22 @@ class sparse_GP(GP): self.Lm = jitchol(self.Kmm) # The rather complex computations of self.A - if self.likelihood.is_heteroscedastic: - assert self.likelihood.D == 1 # TODO: what if the likelihood is heterscedatic and there are multiple independent outputs? - if self.has_uncertain_inputs: -# 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) + if self.has_uncertain_inputs: + if self.likelihood.is_heteroscedastic: + psi2_beta = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0) else: -# tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)) / sf) - tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N))) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) - self.A = tdot(tmp) + psi2_beta = self.psi2.sum(0) * self.likelihood.precision + evals, evecs = linalg.eigh(psi2_beta) + clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable + tmp = evecs * np.sqrt(clipped_evals) else: - if self.has_uncertain_inputs: -# psi2_beta_scaled = (self.psi2 * (self.likelihood.precision / sf2)).sum(0) - 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) + if self.likelihood.is_heteroscedastic: + tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N))) 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) - self.A = tdot(tmp) + tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) + self.A = tdot(tmp) + # factor B # self.B = np.eye(self.M) / sf2 + self.A From 025f31a43fb51423732122f3e1c0121e85016bf9 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 20 May 2013 09:38:08 +0100 Subject: [PATCH 3/3] more tidying in sparse_GP --- GPy/models/sparse_GP.py | 29 +++++------------------------ 1 file changed, 5 insertions(+), 24 deletions(-) diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index ef7b445d..aeff83c2 100644 --- a/GPy/models/sparse_GP.py +++ b/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) - :param kernel : the kernel/covariance function. See link kernels - :type kernel: a GPy kernel + :type likelihood: GPy.likelihood.(Gaussian | EP | Laplace) + :param kernel : the kernel (covariance function). See link kernels + :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,8 +61,6 @@ class sparse_GP(GP): self.psi2 = None def _computations(self): -# sf = self.scale_factor -# sf2 = sf ** 2 # factor Kmm self.Lm = jitchol(self.Kmm) @@ -88,7 +84,6 @@ class sparse_GP(GP): # factor B -# self.B = np.eye(self.M) / sf2 + self.A self.B = np.eye(self.M) + self.A self.LB = jitchol(self.B) @@ -104,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) @@ -115,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[:, None, None] * dL_dpsi2_beta[None, :, :] @@ -141,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))) @@ -149,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 @@ -168,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): @@ -188,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: