mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-10 20:42:39 +02:00
Merge branch 'params' of github.com:SheffieldML/GPy into params
This commit is contained in:
James Hensman
commit
6d2e462b5e
1 changed files with 61 additions and 360 deletions
|
|
@ -3,390 +3,91 @@ from scipy import stats
|
|||
from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot,dtrtrs
|
||||
from likelihood import likelihood
|
||||
|
||||
class EP(likelihood):
|
||||
def __init__(self,data,noise_model):
|
||||
"""
|
||||
Expectation Propagation
|
||||
|
||||
:param data: data to model
|
||||
:type data: numpy array
|
||||
:param noise_model: noise distribution
|
||||
:type noise_model: A GPy noise model
|
||||
|
||||
"""
|
||||
self.noise_model = noise_model
|
||||
self.data = data
|
||||
self.num_data, self.output_dim = self.data.shape
|
||||
self.is_heteroscedastic = True
|
||||
self.num_params = 0
|
||||
|
||||
#Initial values - Likelihood approximation parameters:
|
||||
#p(y|f) = t(f|tau_tilde,v_tilde)
|
||||
self.tau_tilde = np.zeros(self.num_data)
|
||||
self.v_tilde = np.zeros(self.num_data)
|
||||
|
||||
#initial values for the GP variables
|
||||
self.Y = np.zeros((self.num_data,1))
|
||||
self.covariance_matrix = np.eye(self.num_data)
|
||||
self.precision = np.ones(self.num_data)[:,None]
|
||||
self.Z = 0
|
||||
self.YYT = None
|
||||
self.V = self.precision * self.Y
|
||||
self.VVT_factor = self.V
|
||||
self.trYYT = 0.
|
||||
|
||||
super(EP, self).__init__()
|
||||
|
||||
def restart(self):
|
||||
self.tau_tilde = np.zeros(self.num_data)
|
||||
self.v_tilde = np.zeros(self.num_data)
|
||||
self.Y = np.zeros((self.num_data,1))
|
||||
self.covariance_matrix = np.eye(self.num_data)
|
||||
self.precision = np.ones(self.num_data)[:,None]
|
||||
self.Z = 0
|
||||
self.YYT = None
|
||||
self.V = self.precision * self.Y
|
||||
self.VVT_factor = self.V
|
||||
self.trYYT = 0.
|
||||
|
||||
def predictive_values(self,mu,var,full_cov,**noise_args):
|
||||
if full_cov:
|
||||
raise NotImplementedError, "Cannot make correlated predictions with an EP likelihood"
|
||||
return self.noise_model.predictive_values(mu,var,**noise_args)
|
||||
|
||||
def log_predictive_density(self, y_test, mu_star, var_star):
|
||||
"""
|
||||
Calculation of the log predictive density
|
||||
|
||||
.. math:
|
||||
p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})
|
||||
|
||||
:param y_test: test observations (y_{*})
|
||||
:type y_test: (Nx1) array
|
||||
:param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*})
|
||||
:type mu_star: (Nx1) array
|
||||
:param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*})
|
||||
:type var_star: (Nx1) array
|
||||
"""
|
||||
return self.noise_model.log_predictive_density(y_test, mu_star, var_star)
|
||||
|
||||
def _get_params(self):
|
||||
#return np.zeros(0)
|
||||
return self.noise_model._get_params()
|
||||
|
||||
def _get_param_names(self):
|
||||
#return []
|
||||
return self.noise_model._get_param_names()
|
||||
|
||||
def _set_params(self,p):
|
||||
#pass # TODO: the EP likelihood might want to take some parameters...
|
||||
self.noise_model._set_params(p)
|
||||
|
||||
def _gradients(self,partial):
|
||||
#return np.zeros(0) # TODO: the EP likelihood might want to take some parameters...
|
||||
return self.noise_model._gradients(partial)
|
||||
|
||||
def _compute_GP_variables(self):
|
||||
#Variables to be called from GP
|
||||
mu_tilde = self.v_tilde/self.tau_tilde #When calling EP, this variable is used instead of Y in the GP model
|
||||
sigma_sum = 1./self.tau_ + 1./self.tau_tilde
|
||||
mu_diff_2 = (self.v_/self.tau_ - mu_tilde)**2
|
||||
self.Z = np.sum(np.log(self.Z_hat)) + 0.5*np.sum(np.log(sigma_sum)) + 0.5*np.sum(mu_diff_2/sigma_sum) #Normalization constant, aka Z_ep
|
||||
self.Z += 0.5*self.num_data*np.log(2*np.pi)
|
||||
|
||||
self.Y = mu_tilde[:,None]
|
||||
self.YYT = np.dot(self.Y,self.Y.T)
|
||||
self.covariance_matrix = np.diag(1./self.tau_tilde)
|
||||
self.precision = self.tau_tilde[:,None]
|
||||
self.V = self.precision * self.Y
|
||||
self.VVT_factor = self.V
|
||||
self.trYYT = np.trace(self.YYT)
|
||||
|
||||
def fit_full(self, K, epsilon=1e-3,power_ep=[1.,1.]):
|
||||
class EP(object):
|
||||
def __init__(self, epsilon=1e-6, eta=1., delta=1.):
|
||||
"""
|
||||
The expectation-propagation algorithm.
|
||||
For nomenclature see Rasmussen & Williams 2006.
|
||||
|
||||
:param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float)
|
||||
:type epsilon: float
|
||||
:param power_ep: Power EP parameters
|
||||
:type power_ep: list of floats
|
||||
|
||||
:param eta: Power EP thing TODO: Ricardo: what, exactly?
|
||||
:type eta: float64
|
||||
:param delta: Power EP thing TODO: Ricardo: what, exactly?
|
||||
:type delta: float64
|
||||
"""
|
||||
self.epsilon = epsilon
|
||||
self.eta, self.delta = power_ep
|
||||
self.epsilon, self.eta, self.delta = epsilon, eta, delta
|
||||
self.reset()
|
||||
|
||||
def reset(self):
|
||||
self.old_mutilde, self.old_vtilde = None, None
|
||||
|
||||
def inference(self, kern, X, likelihood, Y, Y_metadata=None):
|
||||
|
||||
K = kern.K(X)
|
||||
|
||||
mu_tilde, tau_tilde = self.expectation_propagation()
|
||||
|
||||
|
||||
def expectation_propagation(self, K, Y, Y_metadata, likelihood)
|
||||
|
||||
num_data, data_dim = Y.shape
|
||||
assert data_dim == 1, "This EP methods only works for 1D outputs"
|
||||
|
||||
|
||||
#Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
|
||||
mu = np.zeros(self.num_data)
|
||||
Sigma = K.copy()
|
||||
|
||||
"""
|
||||
Initial values - Cavity distribution parameters:
|
||||
q_(f|mu_,sigma2_) = Product{q_i(f|mu_i,sigma2_i)}
|
||||
sigma_ = 1./tau_
|
||||
mu_ = v_/tau_
|
||||
"""
|
||||
self.tau_ = np.empty(self.num_data,dtype=float)
|
||||
self.v_ = np.empty(self.num_data,dtype=float)
|
||||
|
||||
#Initial values - Marginal moments
|
||||
z = np.empty(self.num_data,dtype=float)
|
||||
self.Z_hat = np.empty(self.num_data,dtype=float)
|
||||
phi = np.empty(self.num_data,dtype=float)
|
||||
mu_hat = np.empty(self.num_data,dtype=float)
|
||||
sigma2_hat = np.empty(self.num_data,dtype=float)
|
||||
Z_hat = np.empty(num_data,dtype=np.float64)
|
||||
mu_hat = np.empty(num_data,dtype=np.float64)
|
||||
sigma2_hat = np.empty(num_data,dtype=np.float64)
|
||||
|
||||
#initial values - Gaussian factors
|
||||
if self.old_mutilde is None:
|
||||
tau_tilde, mu_tilde, v_tilde = np.zeros((3, num_data, num_data))
|
||||
else:
|
||||
assert old_mutilde.size == num_data, "data size mis-match: did you change the data? try resetting!"
|
||||
mu_tilde, v_tilde = self.old_mutilde, self.old_vtilde
|
||||
tau_tilde = v_tilde/mu_tilde
|
||||
|
||||
#Approximation
|
||||
epsilon_np1 = self.epsilon + 1.
|
||||
epsilon_np2 = self.epsilon + 1.
|
||||
self.iterations = 0
|
||||
self.np1 = [self.tau_tilde.copy()]
|
||||
self.np2 = [self.v_tilde.copy()]
|
||||
while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
|
||||
update_order = np.random.permutation(self.num_data)
|
||||
iterations = 0
|
||||
while (epsilon_np1 > self.epsilon) or (epsilon_np2 > self.epsilon):
|
||||
update_order = np.random.permutation(num_data)
|
||||
for i in update_order:
|
||||
#Cavity distribution parameters
|
||||
self.tau_[i] = 1./Sigma[i,i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = mu[i]/Sigma[i,i] - self.eta*self.v_tilde[i]
|
||||
tau_cav = 1./Sigma[i,i] - self.eta*tau_tilde[i]
|
||||
v_cav = mu[i]/Sigma[i,i] - self.eta*v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self.data[i],self.tau_[i],self.v_[i])
|
||||
Z_hat[i], mu_hat[i], sigma2_hat[i] = likelihood.moments_match(Y[i], tau_cav, v_cav, Y_metadata=(None if Y_metadata is None else Y_metadata[i]))
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i])
|
||||
self.tau_tilde[i] += Delta_tau
|
||||
self.v_tilde[i] += Delta_v
|
||||
delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma[i,i])
|
||||
delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma[i,i])
|
||||
tau_tilde[i] += delta_tau
|
||||
v_tilde[i] += delta_v
|
||||
#Posterior distribution parameters update
|
||||
DSYR(Sigma,Sigma[:,i].copy(), -float(Delta_tau/(1.+ Delta_tau*Sigma[i,i])))
|
||||
mu = np.dot(Sigma,self.v_tilde)
|
||||
self.iterations += 1
|
||||
#Sigma recomptutation with Cholesky decompositon
|
||||
Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K
|
||||
B = np.eye(self.num_data) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K
|
||||
DSYR(Sigma, Sigma[:,i].copy(), -Delta_tau/(1.+ Delta_tau*Sigma[i,i]))
|
||||
mu = np.dot(Sigma, v_tilde)
|
||||
iterations += 1
|
||||
|
||||
#(re) compute Sigma and mu using full Cholesky decompy
|
||||
tau_tilde_root = np.sqrt(tau_tilde)
|
||||
Sroot_tilde_K = tau_tilde_root[:,None] * K
|
||||
B = np.eye(num_data) + Sroot_tilde_K * tau_tilde_root[None,:]
|
||||
L = jitchol(B)
|
||||
V,info = dtrtrs(L,Sroot_tilde_K,lower=1)
|
||||
V, _ = 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.num_data
|
||||
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.num_data
|
||||
self.np1.append(self.tau_tilde.copy())
|
||||
self.np2.append(self.v_tilde.copy())
|
||||
mu = np.dot(Sigma,v_tilde)
|
||||
|
||||
return self._compute_GP_variables()
|
||||
#monitor convergence
|
||||
epsilon_np1 = np.mean(np.square(tau_tilde-tau_tilde_old))
|
||||
epsilon_np2 = np.mean(np.square(v_tilde-v_tilde_old))
|
||||
tau_tilde_old = tau_tilde.copy()
|
||||
v_tilde_old = v_tilde.copy()
|
||||
|
||||
def fit_DTC(self, Kmm, Kmn, epsilon=1e-3,power_ep=[1.,1.]):
|
||||
"""
|
||||
The expectation-propagation algorithm with sparse pseudo-input.
|
||||
For nomenclature see ... 2013.
|
||||
return mu, Sigma, mu_tilde, tau_tilde
|
||||
|
||||
:param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float)
|
||||
:type epsilon: float
|
||||
:param power_ep: Power EP parameters
|
||||
:type power_ep: list of floats
|
||||
|
||||
"""
|
||||
self.epsilon = epsilon
|
||||
self.eta, self.delta = power_ep
|
||||
|
||||
num_inducing = Kmm.shape[0]
|
||||
|
||||
#TODO: this doesn't work with uncertain inputs!
|
||||
|
||||
"""
|
||||
Prior approximation parameters:
|
||||
q(f|X) = int_{df}{N(f|KfuKuu_invu,diag(Kff-Qff)*N(u|0,Kuu)} = N(f|0,Sigma0)
|
||||
Sigma0 = Qnn = Knm*Kmmi*Kmn
|
||||
"""
|
||||
KmnKnm = np.dot(Kmn,Kmn.T)
|
||||
Lm = jitchol(Kmm)
|
||||
Lmi = chol_inv(Lm)
|
||||
Kmmi = np.dot(Lmi.T,Lmi)
|
||||
KmmiKmn = np.dot(Kmmi,Kmn)
|
||||
Qnn_diag = np.sum(Kmn*KmmiKmn,-2)
|
||||
LLT0 = Kmm.copy()
|
||||
|
||||
#Kmmi, Lm, Lmi, Kmm_logdet = pdinv(Kmm)
|
||||
#KmnKnm = np.dot(Kmn, Kmn.T)
|
||||
#KmmiKmn = np.dot(Kmmi,Kmn)
|
||||
#Qnn_diag = np.sum(Kmn*KmmiKmn,-2)
|
||||
#LLT0 = Kmm.copy()
|
||||
|
||||
"""
|
||||
Posterior approximation: q(f|y) = N(f| mu, Sigma)
|
||||
Sigma = Diag + P*R.T*R*P.T + K
|
||||
mu = w + P*Gamma
|
||||
"""
|
||||
mu = np.zeros(self.num_data)
|
||||
LLT = Kmm.copy()
|
||||
Sigma_diag = Qnn_diag.copy()
|
||||
|
||||
"""
|
||||
Initial values - Cavity distribution parameters:
|
||||
q_(g|mu_,sigma2_) = Product{q_i(g|mu_i,sigma2_i)}
|
||||
sigma_ = 1./tau_
|
||||
mu_ = v_/tau_
|
||||
"""
|
||||
self.tau_ = np.empty(self.num_data,dtype=float)
|
||||
self.v_ = np.empty(self.num_data,dtype=float)
|
||||
|
||||
#Initial values - Marginal moments
|
||||
z = np.empty(self.num_data,dtype=float)
|
||||
self.Z_hat = np.empty(self.num_data,dtype=float)
|
||||
phi = np.empty(self.num_data,dtype=float)
|
||||
mu_hat = np.empty(self.num_data,dtype=float)
|
||||
sigma2_hat = np.empty(self.num_data,dtype=float)
|
||||
|
||||
#Approximation
|
||||
epsilon_np1 = 1
|
||||
epsilon_np2 = 1
|
||||
self.iterations = 0
|
||||
np1 = [self.tau_tilde.copy()]
|
||||
np2 = [self.v_tilde.copy()]
|
||||
while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
|
||||
update_order = np.random.permutation(self.num_data)
|
||||
for i in update_order:
|
||||
#Cavity distribution parameters
|
||||
self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self.data[i],self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
|
||||
self.tau_tilde[i] += Delta_tau
|
||||
self.v_tilde[i] += Delta_v
|
||||
#Posterior distribution parameters update
|
||||
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 = 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
|
||||
self.iterations += 1
|
||||
#Sigma recomputation with Cholesky decompositon
|
||||
LLT = LLT0 + np.dot(Kmn*self.tau_tilde[None,:],Kmn.T)
|
||||
L = jitchol(LLT)
|
||||
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)
|
||||
epsilon_np1 = sum((self.tau_tilde-np1[-1])**2)/self.num_data
|
||||
epsilon_np2 = sum((self.v_tilde-np2[-1])**2)/self.num_data
|
||||
np1.append(self.tau_tilde.copy())
|
||||
np2.append(self.v_tilde.copy())
|
||||
|
||||
self._compute_GP_variables()
|
||||
|
||||
def fit_FITC(self, Kmm, Kmn, Knn_diag, epsilon=1e-3,power_ep=[1.,1.]):
|
||||
"""
|
||||
The expectation-propagation algorithm with sparse pseudo-input.
|
||||
For nomenclature see Naish-Guzman and Holden, 2008.
|
||||
|
||||
:param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float)
|
||||
:type epsilon: float
|
||||
:param power_ep: Power EP parameters
|
||||
:type power_ep: list of floats
|
||||
"""
|
||||
self.epsilon = epsilon
|
||||
self.eta, self.delta = power_ep
|
||||
|
||||
num_inducing = Kmm.shape[0]
|
||||
|
||||
"""
|
||||
Prior approximation parameters:
|
||||
q(f|X) = int_{df}{N(f|KfuKuu_invu,diag(Kff-Qff)*N(u|0,Kuu)} = N(f|0,Sigma0)
|
||||
Sigma0 = diag(Knn-Qnn) + Qnn, Qnn = Knm*Kmmi*Kmn
|
||||
"""
|
||||
Lm = jitchol(Kmm)
|
||||
Lmi = chol_inv(Lm)
|
||||
Kmmi = np.dot(Lmi.T,Lmi)
|
||||
P0 = Kmn.T
|
||||
KmnKnm = np.dot(P0.T, P0)
|
||||
KmmiKmn = np.dot(Kmmi,P0.T)
|
||||
Qnn_diag = np.sum(P0.T*KmmiKmn,-2)
|
||||
Diag0 = Knn_diag - Qnn_diag
|
||||
R0 = jitchol(Kmmi).T
|
||||
|
||||
"""
|
||||
Posterior approximation: q(f|y) = N(f| mu, Sigma)
|
||||
Sigma = Diag + P*R.T*R*P.T + K
|
||||
mu = w + P*Gamma
|
||||
"""
|
||||
self.w = np.zeros(self.num_data)
|
||||
self.Gamma = np.zeros(num_inducing)
|
||||
mu = np.zeros(self.num_data)
|
||||
P = P0.copy()
|
||||
R = R0.copy()
|
||||
Diag = Diag0.copy()
|
||||
Sigma_diag = Knn_diag
|
||||
RPT0 = np.dot(R0,P0.T)
|
||||
|
||||
"""
|
||||
Initial values - Cavity distribution parameters:
|
||||
q_(g|mu_,sigma2_) = Product{q_i(g|mu_i,sigma2_i)}
|
||||
sigma_ = 1./tau_
|
||||
mu_ = v_/tau_
|
||||
"""
|
||||
self.tau_ = np.empty(self.num_data,dtype=float)
|
||||
self.v_ = np.empty(self.num_data,dtype=float)
|
||||
|
||||
#Initial values - Marginal moments
|
||||
z = np.empty(self.num_data,dtype=float)
|
||||
self.Z_hat = np.empty(self.num_data,dtype=float)
|
||||
phi = np.empty(self.num_data,dtype=float)
|
||||
mu_hat = np.empty(self.num_data,dtype=float)
|
||||
sigma2_hat = np.empty(self.num_data,dtype=float)
|
||||
|
||||
#Approximation
|
||||
epsilon_np1 = 1
|
||||
epsilon_np2 = 1
|
||||
self.iterations = 0
|
||||
self.np1 = [self.tau_tilde.copy()]
|
||||
self.np2 = [self.v_tilde.copy()]
|
||||
while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
|
||||
update_order = np.random.permutation(self.num_data)
|
||||
for i in update_order:
|
||||
#Cavity distribution parameters
|
||||
self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
|
||||
self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
|
||||
#Marginal moments
|
||||
self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self.data[i],self.tau_[i],self.v_[i])
|
||||
#Site parameters update
|
||||
Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
|
||||
Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
|
||||
self.tau_tilde[i] += Delta_tau
|
||||
self.v_tilde[i] += Delta_v
|
||||
#Posterior distribution parameters update
|
||||
dtd1 = Delta_tau*Diag[i] + 1.
|
||||
dii = Diag[i]
|
||||
Diag[i] = dii - (Delta_tau * dii**2.)/dtd1
|
||||
pi_ = P[i,:].reshape(1,num_inducing)
|
||||
P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_
|
||||
Rp_i = np.dot(R,pi_.T)
|
||||
RTR = np.dot(R.T,np.dot(np.eye(num_inducing) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
|
||||
R = jitchol(RTR).T
|
||||
self.w[i] += (Delta_v - Delta_tau*self.w[i])*dii/dtd1
|
||||
self.Gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T)
|
||||
RPT = np.dot(R,P.T)
|
||||
Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1)
|
||||
mu = self.w + np.dot(P,self.Gamma)
|
||||
self.iterations += 1
|
||||
#Sigma recomptutation with Cholesky decompositon
|
||||
Iplus_Dprod_i = 1./(1.+ Diag0 * self.tau_tilde)
|
||||
Diag = Diag0 * Iplus_Dprod_i
|
||||
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 = 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
|
||||
self.Gamma = np.dot(R.T, np.dot(RPT,self.v_tilde))
|
||||
mu = self.w + np.dot(P,self.Gamma)
|
||||
epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.num_data
|
||||
epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.num_data
|
||||
self.np1.append(self.tau_tilde.copy())
|
||||
self.np2.append(self.v_tilde.copy())
|
||||
|
||||
return self._compute_GP_variables()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue