apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-05 14:55:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

lots of F-ordering nonsense. Seems to work though

Browse source
This commit is contained in:
James Hensman 2014-02-14 11:38:34 +00:00
parent 80a734e153
commit 3375606039
4 changed files with 109 additions and 88 deletions
Download patch file Download diff file Expand all files Collapse all files

2
GPy/core/sparse_gp.py
View file

@ -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

2
GPy/inference/latent_function_inference/posterior.py
View file

@ -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

22
GPy/inference/latent_function_inference/varDTC.py
View file

@ -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)

171
GPy/util/linalg.py
View file

@ -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)"