diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index bbbea2fe..09e84c12 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -2,7 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np -from ..util.linalg import mdot, tdot, symmetrify, backsub_both_sides, chol_inv, dtrtrs, dpotrs, dpotri +from ..util.linalg import mdot, tdot, symmetrify, backsub_both_sides, dtrtrs, dpotrs, dpotri from gp import GP from parameterization.param import Param from ..inference.latent_function_inference import varDTC diff --git a/GPy/inference/latent_function_inference/posterior.py b/GPy/inference/latent_function_inference/posterior.py index 09ac96e8..f28bf9d1 100644 --- a/GPy/inference/latent_function_inference/posterior.py +++ b/GPy/inference/latent_function_inference/posterior.py @@ -115,7 +115,7 @@ class Posterior(object): @property def woodbury_inv(self): if self._woodbury_inv is None: - self._woodbury_inv, _ = dpotri(self.woodbury_chol, lower=0) + self._woodbury_inv, _ = dpotri(self.woodbury_chol, lower=1) #self._woodbury_inv, _ = dpotrs(self.woodbury_chol, np.eye(self.woodbury_chol.shape[0]), lower=1) symmetrify(self._woodbury_inv) return self._woodbury_inv diff --git a/GPy/inference/latent_function_inference/varDTC.py b/GPy/inference/latent_function_inference/varDTC.py index d7f770c8..90ca7a6a 100644 --- a/GPy/inference/latent_function_inference/varDTC.py +++ b/GPy/inference/latent_function_inference/varDTC.py @@ -91,7 +91,7 @@ class VarDTC(object): tmp = psi1 * (np.sqrt(beta.reshape(num_data, 1))) else: tmp = psi1 * (np.sqrt(beta)) - tmp, _ = dtrtrs(Lm, np.asfortranarray(tmp.T), lower=1) + tmp, _ = dtrtrs(Lm, tmp.T, lower=1) A = tdot(tmp) # factor B @@ -100,14 +100,16 @@ class VarDTC(object): LB = jitchol(B) # 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) + #self.YYTfactor = self.get_YYTfactor(Y) + #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) + VVT_factor = beta*param_to_array(Y) trYYT = self.get_trYYT(Y) + psi1Vf = np.dot(psi1.T, VVT_factor) # back substutue C into psi1Vf - tmp, info1 = dtrtrs(Lm, np.asfortranarray(psi1Vf), lower=1, trans=0) - _LBi_Lmi_psi1Vf, _ = dtrtrs(LB, np.asfortranarray(tmp), lower=1, trans=0) + tmp, info1 = dtrtrs(Lm, psi1Vf, lower=1, trans=0) + _LBi_Lmi_psi1Vf, _ = dtrtrs(LB, tmp, lower=1, trans=0) tmp, info2 = dtrtrs(LB, _LBi_Lmi_psi1Vf, lower=1, trans=1) Cpsi1Vf, info3 = dtrtrs(Lm, tmp, lower=1, trans=1) @@ -152,8 +154,8 @@ class VarDTC(object): raise NotImplementedError, "heteroscedatic derivates with uncertain inputs not implemented" else: LBi = chol_inv(LB) - Lmi_psi1, nil = dtrtrs(Lm, np.asfortranarray(psi1.T), lower=1, trans=0) - _LBi_Lmi_psi1, _ = dtrtrs(LB, np.asfortranarray(Lmi_psi1), lower=1, trans=0) + Lmi_psi1, nil = dtrtrs(Lm, psi1.T, lower=1, trans=0) + _LBi_Lmi_psi1, _ = dtrtrs(LB, Lmi_psi1, lower=1, trans=0) partial_for_likelihood = -0.5 * beta + 0.5 * likelihood.V**2 partial_for_likelihood += 0.5 * output_dim * (psi0 - np.sum(Lmi_psi1**2,0))[:,None] * beta**2 @@ -195,12 +197,14 @@ class VarDTC(object): if VVT_factor.shape[1] == Y.shape[1]: woodbury_vector = Cpsi1Vf # == Cpsi1V else: + print 'foobar' psi1V = np.dot(Y.T*beta, psi1).T - tmp, _ = dtrtrs(Lm, np.asfortranarray(psi1V), lower=1, trans=0) + 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=0) + Bi, _ = dpotri(LB, lower=1) symmetrify(Bi) + Bi = dpotri(LB, lower=1)[0] woodbury_inv = backsub_both_sides(Lm, np.eye(num_inducing) - Bi) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 44f3700d..6bcddc3e 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -17,7 +17,8 @@ import os from config import * if np.all(np.float64((scipy.__version__).split('.')[:2]) >= np.array([0, 12])): - import scipy.linalg.lapack as lapack + #import scipy.linalg.lapack.clapack as lapack + from scipy.linalg import lapack else: from scipy.linalg.lapack import flapack as lapack @@ -41,42 +42,103 @@ else: _blas_available = False warnings.warn("warning: caught this exception:" + str(e)) -def dtrtri(L, lower=0): +def force_F_ordered_symmetric(A): """ - Wrapper for lapack dtrtrs function + return a F ordered version of A, assuming A is symmetric + """ + if A.flags['F_CONTIGUOUS']: + return A + if A.flags['C_CONTIGUOUS']: + return A.T + else: + return np.asfortranarray(A) + +def force_F_ordered(A): + """ + return a F ordered version of A, assuming A is triangular + """ + if A.flags['F_CONTIGUOUS']: + return A + print "why are your arrays not F order?" + return np.asfortranarray(A) + +def jitchol(A, maxtries=5): + A = force_F_ordered_symmetric(A) + L, info = lapack.dpotrf(A, lower=1) + if info == 0: + return L + else: + if maxtries==0: + raise linalg.LinAlgError, "not positive definite, even with jitter." + diagA = np.diag(A) + if np.any(diagA <= 0.): + raise linalg.LinAlgError, "not pd: non-positive diagonal elements" + jitter = diagA.mean() * 1e-6 + + return jitchol(A+np.eye(A.shape[0])*jitter, maxtries-1) + +#def jitchol(A, maxtries=5): +# A = np.ascontiguousarray(A) +# L, info = lapack.dpotrf(A, lower=1) +# if info == 0: +# return L +# else: +# diagA = np.diag(A) +# if np.any(diagA <= 0.): +# raise linalg.LinAlgError, "not pd: non-positive diagonal elements" +# jitter = diagA.mean() * 1e-6 +# while maxtries > 0 and np.isfinite(jitter): +# print 'Warning: adding jitter of {:.10e}'.format(jitter) +# try: +# return linalg.cholesky(A + np.eye(A.shape[0]).T * jitter, lower=True) +# except: +# jitter *= 10 +# finally: +# maxtries -= 1 +# raise linalg.LinAlgError, "not positive definite, even with jitter." +# + + + + +def dtrtri(L, lower=1): + """ + Wrapper for lapack dtrtri function Inverse of L :param L: Triangular Matrix L :param lower: is matrix lower (true) or upper (false) :returns: Li, info """ + L = force_F_ordered(L) return lapack.dtrtri(L, lower=lower) -def dtrtrs(A, B, lower=0, trans=0, unitdiag=0): +def dtrtrs(A, B, lower=1, trans=0, unitdiag=0): """ Wrapper for lapack dtrtrs function - :param A: Matrix A + :param A: Matrix A(triangular) :param B: Matrix B :param lower: is matrix lower (true) or upper (false) :returns: """ + A = force_F_ordered(A) + #Note: B does not seem to need to be F ordered! return lapack.dtrtrs(A, B, lower=lower, trans=trans, unitdiag=unitdiag) -def dpotrs(A, B, lower=0): +def dpotrs(A, B, lower=1): """ Wrapper for lapack dpotrs function - :param A: Matrix A :param B: Matrix B :param lower: is matrix lower (true) or upper (false) :returns: - """ + A = force_F_ordered(A) return lapack.dpotrs(A, B, lower=lower) -def dpotri(A, lower=0): +def dpotri(A, lower=1): """ Wrapper for lapack dpotri function @@ -85,7 +147,12 @@ def dpotri(A, lower=0): :returns: A inverse """ - return lapack.dpotri(A, lower=lower) + assert lower==1, "scipy linalg behaviour is very weird. please use lower, fortran ordered arrays" + + A = force_F_ordered(A) + R, info = lapack.dpotri(A, lower=0) + symmetrify(R) + return R, info def pddet(A): """ @@ -133,60 +200,8 @@ def _mdot_r(a, b): b = b[0] return np.dot(a, b) -def jitchol(A, maxtries=5): - A = np.asfortranarray(A) - L, info = lapack.dpotrf(A, lower=1) - if info == 0: - return L - else: - diagA = np.diag(A) - if np.any(diagA <= 0.): - raise linalg.LinAlgError, "not pd: non-positive diagonal elements" - jitter = diagA.mean() * 1e-6 - while maxtries > 0 and np.isfinite(jitter): - print 'Warning: adding jitter of {:.10e}'.format(jitter) - try: - return linalg.cholesky(A + np.eye(A.shape[0]).T * jitter, lower=True) - except: - jitter *= 10 - finally: - maxtries -= 1 - raise linalg.LinAlgError, "not positive definite, even with jitter." - - - -def jitchol_old(A, maxtries=5): - """ - :param A: An almost pd square matrix - - :rval L: the Cholesky decomposition of A - - .. note: - - Adds jitter to K, to enforce positive-definiteness - if stuff breaks, please check: - np.allclose(sp.linalg.cholesky(XXT, lower = True), np.triu(sp.linalg.cho_factor(XXT)[0]).T) - - """ - try: - return linalg.cholesky(A, lower=True) - except linalg.LinAlgError: - diagA = np.diag(A) - if np.any(diagA < 0.): - 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), - try: - return linalg.cholesky(A + np.eye(A.shape[0]).T * jitter, lower=True) - except: - jitter *= 10 - - raise linalg.LinAlgError, "not positive definite, even with jitter." - def pdinv(A, *args): """ - :param A: A DxD pd numpy array :rval Ai: the inverse of A @@ -201,7 +216,7 @@ def pdinv(A, *args): """ L = jitchol(A, *args) logdet = 2.*np.sum(np.log(np.diag(L))) - Li = chol_inv(L) + Li = dtrtri(L) Ai, _ = lapack.dpotri(L) # Ai = np.tril(Ai) + np.tril(Ai,-1).T symmetrify(Ai) @@ -209,7 +224,7 @@ def pdinv(A, *args): return Ai, L, Li, logdet -def chol_inv(L): +def dtrtri(L): """ Inverts a Cholesky lower triangular matrix @@ -218,7 +233,8 @@ def chol_inv(L): """ - return lapack.dtrtri(L, lower=True)[0] + L = force_F_ordered(L) + return lapack.dtrtri(L, lower=1)[0] def multiple_pdinv(A): @@ -234,7 +250,7 @@ def multiple_pdinv(A): N = A.shape[-1] chols = [jitchol(A[:, :, i]) for i in range(N)] halflogdets = [np.sum(np.log(np.diag(L[0]))) for L in chols] - invs = [lapack.dpotri(L[0], True)[0] for L in chols] + invs = [dpotri(L[0], True)[0] for L in chols] invs = [np.triu(I) + np.triu(I, 1).T for I in invs] return np.dstack(invs), np.array(halflogdets) @@ -425,7 +441,8 @@ def tdot_blas(mat, out=None): symmetrify(out, upper=True) - return out + + return np.ascontiguousarray(out) def tdot(*args, **kwargs): if _blas_available: @@ -557,24 +574,24 @@ def cholupdate(L, x): def backsub_both_sides(L, X, transpose='left'): """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky""" if transpose == 'left': - tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=1) - return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=1)[0].T + tmp, _ = dtrtrs(L, X, lower=1, trans=1) + return dtrtrs(L, tmp.T, lower=1, trans=1)[0].T else: - tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=0) - return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=0)[0].T + tmp, _ = dtrtrs(L, X, lower=1, trans=0) + return dtrtrs(L, tmp.T, lower=1, trans=0)[0].T def PCA(Y, input_dim): """ -Principal component analysis: maximum likelihood solution by SVD + Principal component analysis: maximum likelihood solution by SVD -:param Y: NxD np.array of data -:param input_dim: int, dimension of projection + :param Y: NxD np.array of data + :param input_dim: int, dimension of projection -:rval X: - Nxinput_dim np.array of dimensionality reduced data -:rval W: - input_dimxD mapping from X to Y + :rval X: - Nxinput_dim np.array of dimensionality reduced data + :rval W: - input_dimxD mapping from X to Y -""" + """ if not np.allclose(Y.mean(axis=0), 0.0): print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"