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

Factored out lapack into utils so we can check version and give deprecation warnings

Browse source
This commit is contained in:
Alan Saul 2013-06-06 14:04:14 +01:00
parent b142b68876
commit e587618339
6 changed files with 73 additions and 38 deletions
Download patch file Download diff file Expand all files Collapse all files

16
GPy/core/fitc.py
View file

@ -3,10 +3,10 @@
import numpy as np
import pylab as pb
from ..util.linalg import mdot, jitchol, chol_inv, tdot, symmetrify,pdinv
from ..util.linalg import mdot, jitchol, chol_inv, tdot, symmetrify, pdinv, dtrtrs
from ..util.plot import gpplot
from .. import kern
from scipy import stats, linalg
from scipy import stats
from sparse_gp import SparseGP
class FITC(SparseGP):
@ -50,7 +50,7 @@ class FITC(SparseGP):
def _computations(self):
#factor Kmm
self.Lm = jitchol(self.Kmm)
self.Lmi,info = linalg.lapack.flapack.dtrtrs(self.Lm,np.eye(self.num_inducing),lower=1)
self.Lmi,info = dtrtrs(self.Lm,np.eye(self.num_inducing),lower=1)
Lmipsi1 = np.dot(self.Lmi,self.psi1)
self.Qnn = np.dot(Lmipsi1.T,Lmipsi1).copy()
self.Diag0 = self.psi0 - np.diag(self.Qnn)
@ -59,7 +59,7 @@ class FITC(SparseGP):
# The rather complex computations of self.A
tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.num_data)))
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
tmp, _ = dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
self.A = tdot(tmp)
# factor B
@ -68,8 +68,8 @@ class FITC(SparseGP):
self.LBi = chol_inv(self.LB)
self.psi1V = np.dot(self.psi1, self.V_star)
Lmi_psi1V, info = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0)
self._LBi_Lmi_psi1V, _ = linalg.lapack.flapack.dtrtrs(self.LB, np.asfortranarray(Lmi_psi1V), lower=1, trans=0)
Lmi_psi1V, info = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0)
self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(Lmi_psi1V), lower=1, trans=0)
Kmmipsi1 = np.dot(self.Lmi.T,Lmipsi1)
b_psi1_Ki = self.beta_star * Kmmipsi1.T
@ -188,7 +188,7 @@ class FITC(SparseGP):
self.P = Iplus_Dprod_i[:,None] * self.psi1.T
self.RPT0 = np.dot(self.Lmi,self.psi1)
self.L = np.linalg.cholesky(np.eye(self.num_inducing) + np.dot(self.RPT0,((1. - Iplus_Dprod_i)/self.Diag0)[:,None]*self.RPT0.T))
self.R,info = linalg.lapack.flapack.dtrtrs(self.L,self.Lmi,lower=1)
self.R,info = dtrtrs(self.L,self.Lmi,lower=1)
self.RPT = np.dot(self.R,self.P.T)
self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT)
self.w = self.Diag * self.likelihood.v_tilde
@ -210,7 +210,7 @@ class FITC(SparseGP):
# = I - [RPT0] * (U*U.T)^-1 * [RPT0].T
# = I - V.T * V
U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn)
V,info = linalg.lapack.flapack.dtrtrs(U,self.RPT0.T,lower=1)
V,info = dtrtrs(U,self.RPT0.T,lower=1)
C = np.eye(self.num_inducing) - np.dot(V.T,V)
mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:])
#self.C = C

13
GPy/core/gp.py
View file

@ -3,10 +3,9 @@
import numpy as np
from scipy import linalg
import pylab as pb
from .. import kern
from ..util.linalg import pdinv, mdot, tdot
from ..util.linalg import pdinv, mdot, tdot, dpotrs, dtrtrs
#from ..util.plot import gpplot, Tango
from ..likelihoods import EP
from gp_base import GPBase
@ -44,13 +43,13 @@ class GP(GPBase):
# the gradient of the likelihood wrt the covariance matrix
if self.likelihood.YYT is None:
#alpha = np.dot(self.Ki, self.likelihood.Y)
alpha,_ = linalg.lapack.flapack.dpotrs(self.L, self.likelihood.Y,lower=1)
alpha,_ = dpotrs(self.L, self.likelihood.Y,lower=1)
self.dL_dK = 0.5 * (tdot(alpha) - self.output_dim * self.Ki)
else:
#tmp = mdot(self.Ki, self.likelihood.YYT, self.Ki)
tmp, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(self.likelihood.YYT), lower=1)
tmp, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(tmp.T), lower=1)
tmp, _ = dpotrs(self.L, np.asfortranarray(self.likelihood.YYT), lower=1)
tmp, _ = dpotrs(self.L, np.asfortranarray(tmp.T), lower=1)
self.dL_dK = 0.5 * (tmp - self.output_dim * self.Ki)
def _get_params(self):
@ -76,7 +75,7 @@ class GP(GPBase):
Computes the model fit using YYT if it's available
"""
if self.likelihood.YYT is None:
tmp, _ = linalg.lapack.flapack.dtrtrs(self.L, np.asfortranarray(self.likelihood.Y), lower=1)
tmp, _ = dtrtrs(self.L, np.asfortranarray(self.likelihood.Y), lower=1)
return -0.5 * np.sum(np.square(tmp))
#return -0.5 * np.sum(np.square(np.dot(self.Li, self.likelihood.Y)))
else:
@ -108,7 +107,7 @@ class GP(GPBase):
"""
Kx = self.kern.K(_Xnew,self.X,which_parts=which_parts).T
#KiKx = np.dot(self.Ki, Kx)
KiKx, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(Kx), lower=1)
KiKx, _ = dpotrs(self.L, np.asfortranarray(Kx), lower=1)
mu = np.dot(KiKx.T, self.likelihood.Y)
if full_cov:
Kxx = self.kern.K(_Xnew, which_parts=which_parts)

14
GPy/core/sparse_gp.py
View file

@ -3,7 +3,7 @@
import numpy as np
import pylab as pb
from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides, chol_inv
from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides, chol_inv, dtrtrs, dpotrs, dpotri
from scipy import linalg
from ..likelihoods import Gaussian
from gp_base import GPBase
@ -80,7 +80,7 @@ class SparseGP(GPBase):
tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.num_data)))
else:
tmp = self.psi1 * (np.sqrt(self.likelihood.precision))
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
tmp, _ = dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
self.A = tdot(tmp)
@ -92,10 +92,10 @@ class SparseGP(GPBase):
self.psi1V = np.dot(self.psi1, self.likelihood.V)
# back substutue C into psi1V
tmp, info1 = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0)
self._LBi_Lmi_psi1V, _ = linalg.lapack.flapack.dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
tmp, info2 = linalg.lapack.flapack.dpotrs(self.LB, tmp, lower=1)
self.Cpsi1V, info3 = linalg.lapack.flapack.dtrtrs(self.Lm, tmp, lower=1, trans=1)
tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0)
self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
tmp, info2 = dpotrs(self.LB, tmp, lower=1)
self.Cpsi1V, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1)
# Compute dL_dKmm
tmp = tdot(self._LBi_Lmi_psi1V)
@ -219,7 +219,7 @@ class SparseGP(GPBase):
def _raw_predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False):
"""Internal helper function for making predictions, does not account for normalization"""
Bi, _ = linalg.lapack.flapack.dpotri(self.LB, lower=0) # WTH? this lower switch should be 1, but that doesn't work!
Bi, _ = dpotri(self.LB, lower=0) # WTH? this lower switch should be 1, but that doesn't work!
symmetrify(Bi)
Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi)

14
GPy/likelihoods/ep.py
View file

@ -1,6 +1,6 @@
import numpy as np
from scipy import stats, linalg
from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot
from scipy import stats
from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot,dtrtrs
from likelihood import likelihood
class EP(likelihood):
@ -124,7 +124,7 @@ class EP(likelihood):
Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K
B = np.eye(self.N) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K
L = jitchol(B)
V,info = linalg.lapack.flapack.dtrtrs(L,Sroot_tilde_K,lower=1)
V,info = dtrtrs(L,Sroot_tilde_K,lower=1)
Sigma = K - np.dot(V.T,V)
mu = np.dot(Sigma,self.v_tilde)
epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N
@ -210,7 +210,7 @@ class EP(likelihood):
DSYR(LLT,Kmn[:,i].copy(),Delta_tau) #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
L = jitchol(LLT)
#cholUpdate(L,Kmn[:,i]*np.sqrt(Delta_tau))
V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1)
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
@ -218,8 +218,8 @@ class EP(likelihood):
#Sigma recomputation with Cholesky decompositon
LLT = LLT0 + np.dot(Kmn*self.tau_tilde[None,:],Kmn.T)
L = jitchol(LLT)
V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1)
V2,info = linalg.lapack.flapack.dtrtrs(L.T,V,lower=0)
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,self.v_tilde)
mu = np.dot(V2.T,Knmv_tilde)
@ -322,7 +322,7 @@ class EP(likelihood):
P = Iplus_Dprod_i[:,None] * P0
safe_diag = np.where(Diag0 < self.tau_tilde, self.tau_tilde/(1.+Diag0*self.tau_tilde), (1. - Iplus_Dprod_i)/Diag0)
L = jitchol(np.eye(num_inducing) + np.dot(RPT0,safe_diag[:,None]*RPT0.T))
R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1)
R,info = dtrtrs(L,R0,lower=1)
RPT = np.dot(R,P.T)
Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1)
self.w = Diag * self.v_tilde

2
GPy/testing/examples_tests.py
View file

@ -39,7 +39,7 @@ def model_instance(model):
#assert isinstance(model, GPy.core.model)
return isinstance(model, GPy.core.Model)
#@nottest
@nottest
def test_models():
examples_path = os.path.dirname(GPy.examples.__file__)
# Load modules

52
GPy/util/linalg.py
View file

@ -11,6 +11,12 @@ import types
import ctypes
from ctypes import byref, c_char, c_int, c_double # TODO
# import scipy.lib.lapack
import scipy
if np.all(np.float64((scipy.__version__).split('.')[:2]) >= np.array([0, 10])):
import scipy.linalg.lapack as lapack
else:
import scipy.linalg.lapack.flapack as lapack
try:
_blaslib = ctypes.cdll.LoadLibrary(np.core._dotblas.__file__) # @UndefinedVariable
@ -18,6 +24,36 @@ try:
except:
_blas_available = False
def dtrtrs(A, B, lower=0, trans=0, unitdiag=0, overwrite_b=0):
"""Wrapper for lapack dtrtrs function
:param A: Matrix A
:param B: Matrix B
:param lower: is matrix lower (true) or upper (false)
:returns:
"""
return lapack.dtrtrs(A, B, lower=lower, trans=trans, unitdiag=unitdiag, overwrite_b=overwrite_b)
def dpotrs(A, B, overwrite_b=0, lower=0):
"""Wrapper for lapack dpotrs function
:param A: Matrix A
:param B: Matrix B
:param lower: is matrix lower (true) or upper (false)
:returns:
"""
return lapack.dpotrs(A, B, overwrite_b=overwrite_b, lower=lower)
def dpotri(A, B, overwrite_b=0, lower=0):
"""Wrapper for lapack dpotri function
:param A: Matrix A
:param B: Matrix B
:param lower: is matrix lower (true) or upper (false)
:returns:
"""
return lapack.dpotri(A, B, overwrite_b=overwrite_b, lower=lower)
def trace_dot(a, b):
"""
efficiently compute the trace of the matrix product of a and b
@ -56,7 +92,7 @@ def _mdot_r(a, b):
def jitchol(A, maxtries=5):
A = np.asfortranarray(A)
L, info = linalg.lapack.flapack.dpotrf(A, lower=1)
L, info = lapack.dpotrf(A, lower=1)
if info == 0:
return L
else:
@ -117,7 +153,7 @@ def pdinv(A, *args):
L = jitchol(A, *args)
logdet = 2.*np.sum(np.log(np.diag(L)))
Li = chol_inv(L)
Ai, _ = linalg.lapack.flapack.dpotri(L)
Ai, _ = lapack.dpotri(L)
# Ai = np.tril(Ai) + np.tril(Ai,-1).T
symmetrify(Ai)
@ -133,7 +169,7 @@ def chol_inv(L):
"""
return linalg.lapack.flapack.dtrtri(L, lower=True)[0]
return lapack.dtrtri(L, lower=True)[0]
def multiple_pdinv(A):
@ -150,7 +186,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 = [linalg.lapack.flapack.dpotri(L[0], True)[0] for L in chols]
invs = [lapack.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)
@ -351,9 +387,9 @@ 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, _ = linalg.lapack.flapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=1)
return linalg.lapack.flapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=1)[0].T
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
else:
tmp, _ = linalg.lapack.flapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=0)
return linalg.lapack.flapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=0)[0].T
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