mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-06-11 15:15:15 +02:00
Broken it by getting rid of squeeze, but now working on making it faster using proper vector multiplciation for diagonals
This commit is contained in:
Alan Saul
parent
20227fb2ac
commit
f9857e08c0
4 changed files with 69 additions and 70 deletions
|
|
@ -37,9 +37,10 @@ def timing():
|
|||
|
||||
def debug_student_t_noise_approx():
|
||||
plot = False
|
||||
real_var = 0.4
|
||||
real_var = 0.1
|
||||
#Start a function, any function
|
||||
X = np.linspace(0.0, 10.0, 100)[:, None]
|
||||
#X = np.array([0.5])[:, None]
|
||||
Y = np.sin(X) + np.random.randn(*X.shape)*real_var
|
||||
|
||||
X_full = np.linspace(0.0, 10.0, 500)[:, None]
|
||||
|
|
@ -52,7 +53,7 @@ def debug_student_t_noise_approx():
|
|||
real_sd = np.sqrt(real_var)
|
||||
print "Real noise: ", real_sd
|
||||
|
||||
initial_var_guess = 1
|
||||
initial_var_guess = 0.02
|
||||
#t_rv = t(deg_free, loc=0, scale=real_var)
|
||||
#noise = t_rvrvs(size=Y.shape)
|
||||
#Y += noise
|
||||
|
|
@ -91,12 +92,14 @@ def debug_student_t_noise_approx():
|
|||
#m.constrain_positive('rbf')
|
||||
#m.constrain_fixed('rbf_v', 1.0898)
|
||||
#m.constrain_fixed('rbf_l', 1.8651)
|
||||
m.constrain_positive('t_noi')
|
||||
#m.constrain_fixed('t_noise_variance', real_sd)
|
||||
m.constrain_positive('rbf')
|
||||
m.constrain_fixed('t_noi', real_sd)
|
||||
m.ensure_default_constraints()
|
||||
m.update_likelihood_approximation()
|
||||
m.optimize(messages=True)
|
||||
print(m)
|
||||
return m
|
||||
#return m
|
||||
#m.optimize('lbfgsb', messages=True, callback=m._update_params_callback)
|
||||
if plot:
|
||||
plt.suptitle('Student-t likelihood')
|
||||
|
|
@ -104,6 +107,7 @@ def debug_student_t_noise_approx():
|
|||
m.plot()
|
||||
plt.plot(X_full, Y_full)
|
||||
plt.ylim(-2.5, 2.5)
|
||||
return m
|
||||
|
||||
#print "Clean student t, ncg"
|
||||
#t_distribution = GPy.likelihoods.likelihood_functions.student_t(deg_free, sigma=edited_real_sd)
|
||||
|
|
|
|||
|
|
@ -53,7 +53,7 @@ class Laplace(likelihood):
|
|||
|
||||
def predictive_values(self, mu, var, full_cov):
|
||||
if full_cov:
|
||||
raise NotImplementedError("Cannot make correlated predictions with an EP likelihood")
|
||||
raise NotImplementedError("Cannot make correlated predictions with an Laplace likelihood")
|
||||
return self.likelihood_function.predictive_values(mu, var)
|
||||
|
||||
def _get_params(self):
|
||||
|
|
@ -63,42 +63,28 @@ class Laplace(likelihood):
|
|||
return self.likelihood_function._get_param_names()
|
||||
|
||||
def _set_params(self, p):
|
||||
#print "Setting laplace param with: ", p
|
||||
return self.likelihood_function._set_params(p)
|
||||
|
||||
def both_gradients(self, dL_d_K_Sigma, dK_dthetaK):
|
||||
"""
|
||||
Find the gradients of the marginal likelihood w.r.t both thetaK and thetaL
|
||||
|
||||
dL_dthetaK differs from that of normal likelihoods as it has additional terms coming from
|
||||
changes to y_tilde and changes to Sigma_tilde when the kernel parameters are adjusted
|
||||
|
||||
Similar terms arise when finding the gradients with respect to changes in the liklihood
|
||||
parameters
|
||||
"""
|
||||
return (self._Kgradients(dL_d_K_Sigma, dK_dthetaK), self._gradients(dL_d_K_Sigma))
|
||||
|
||||
def _shared_gradients_components(self):
|
||||
#FIXME: Careful of side effects! And make sure W and K are up to date!
|
||||
Ki, _, _, _ = pdinv(self.K)
|
||||
d3lik_d3fhat = self.likelihood_function.d3lik_d3f(self.data, self.f_hat)
|
||||
#dL_dfhat = -0.5*np.diag(self.Ki_W_i)*d3lik_d3fhat
|
||||
dL_dfhat = -0.5*(np.diag(self.Ki_W_i)*d3lik_d3fhat)[:, None]
|
||||
Wi_K_i = mdot(self.W_12, self.Bi, self.W_12) #same as rasms R
|
||||
I_KW_i = np.eye(self.N) - np.dot(self.K, Wi_K_i)
|
||||
return dL_dfhat, Ki, I_KW_i, Wi_K_i
|
||||
return dL_dfhat, I_KW_i, Wi_K_i
|
||||
|
||||
def _Kgradients(self, dL_d_K_Sigma, dK_dthetaK):
|
||||
def _Kgradients(self, dK_dthetaK):
|
||||
"""
|
||||
Gradients with respect to prior kernel parameters
|
||||
"""
|
||||
dL_dfhat, Ki, I_KW_i, Wi_K_i = self._shared_gradients_components()
|
||||
dL_dfhat, I_KW_i, Wi_K_i = self._shared_gradients_components()
|
||||
dlp = self.likelihood_function.dlik_df(self.data, self.f_hat)[:, None]
|
||||
|
||||
dL_dthetaK = np.zeros(dK_dthetaK.shape)
|
||||
for thetaK_i, dK_dthetaK_i in enumerate(dK_dthetaK):
|
||||
#Explicit
|
||||
dL_dthetaK[thetaK_i] = 0.5*mdot(self.f_hat.T, Ki, dK_dthetaK_i, Ki, self.f_hat) - 0.5*np.trace(Wi_K_i*dK_dthetaK_i)
|
||||
f_Ki_dK_dtheta_Ki_f = mdot(self.Ki_f.T, dK_dthetaK_i, self.Ki_f)
|
||||
dL_dthetaK[thetaK_i] = 0.5*f_Ki_dK_dtheta_Ki_f - 0.5*np.trace(Wi_K_i*dK_dthetaK_i)
|
||||
#Implicit
|
||||
df_hat_dthetaK = mdot(I_KW_i, dK_dthetaK_i, dlp)
|
||||
dL_dthetaK[thetaK_i] += np.dot(dL_dfhat.T, df_hat_dthetaK)
|
||||
|
|
@ -109,11 +95,12 @@ class Laplace(likelihood):
|
|||
"""
|
||||
Gradients with respect to likelihood parameters
|
||||
"""
|
||||
dL_dfhat, Ki, I_KW_i, Wi_K_i = self._shared_gradients_components()
|
||||
return np.zeros(1)
|
||||
#return np.zeros(0)
|
||||
dL_dfhat, I_KW_i, Wi_K_i = self._shared_gradients_components()
|
||||
dlik_dthetaL, dlik_grad_dthetaL, dlik_hess_dthetaL = self.likelihood_function._gradients(self.data, self.f_hat)
|
||||
|
||||
num_params = len(dlik_dthetaL)
|
||||
#Ki_W_i = np.diag(inv(Ki + self.W))[:, None]
|
||||
dL_dthetaL = np.zeros((1, num_params)) # make space for one derivative for each likelihood parameter
|
||||
for thetaL_i in range(num_params):
|
||||
#Explicit
|
||||
|
|
@ -123,7 +110,6 @@ class Laplace(likelihood):
|
|||
dL_dthetaL[thetaL_i] = np.sum(dlik_dthetaL[thetaL_i]) - 0.5*np.dot(np.diag(self.Ki_W_i), dlik_hess_dthetaL[thetaL_i])
|
||||
#Implicit
|
||||
df_hat_dthetaL = mdot(I_KW_i, self.K, dlik_grad_dthetaL[thetaL_i])
|
||||
import ipdb; ipdb.set_trace() ### XXX BREAKPOINT
|
||||
dL_dthetaL[thetaL_i] += np.dot(dL_dfhat.T, df_hat_dthetaL)
|
||||
|
||||
return np.squeeze(dL_dthetaL) #should be array of length *params-being optimized*, for student t just optimising 1 parameter, this is (1,)
|
||||
|
|
@ -230,10 +216,8 @@ class Laplace(likelihood):
|
|||
self._compute_likelihood_variables()
|
||||
|
||||
def _compute_likelihood_variables(self):
|
||||
#At this point get the hessian matrix
|
||||
#print "Data: ", self.data
|
||||
#print "fhat: ", self.f_hat
|
||||
self.W = -np.diag(self.likelihood_function.d2lik_d2f(self.data, self.f_hat, extra_data=self.extra_data))
|
||||
#At this point get the hessian matrix (or vector as W is diagonal)
|
||||
self.W = -self.likelihood_function.d2lik_d2f(self.data, self.f_hat, extra_data=self.extra_data)
|
||||
|
||||
if not self.likelihood_function.log_concave:
|
||||
self.W[self.W < 0] = 1e-6 # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
|
||||
|
|
@ -273,7 +257,8 @@ class Laplace(likelihood):
|
|||
"""
|
||||
#W is diagnoal so its sqrt is just the sqrt of the diagonal elements
|
||||
W_12 = np.sqrt(W)
|
||||
B = np.eye(K.shape[0]) + np.dot(W_12, np.dot(K, W_12))
|
||||
assert np.all(W_12.T*K*W_12 == np.dot(np.diagflat(W_12), np.dot(K, np.diagflat(W_12)))) # FIXME Take this out when you've done multiinput
|
||||
B = np.eye(K.shape[0]) + W_12.T*K*W_12
|
||||
L = jitchol(B)
|
||||
return (B, L, W_12)
|
||||
|
||||
|
|
@ -330,7 +315,7 @@ class Laplace(likelihood):
|
|||
i = 0
|
||||
while difference > epsilon and i < MAX_ITER and rs < MAX_RESTART:
|
||||
#f_old = f.copy()
|
||||
W = -np.diag(self.likelihood_function.d2lik_d2f(self.data, f, extra_data=self.extra_data))
|
||||
W = -self.likelihood_function.d2lik_d2f(self.data, f, extra_data=self.extra_data)
|
||||
if not self.likelihood_function.log_concave:
|
||||
W[W < 0] = 1e-6 # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur
|
||||
# If the likelihood is non-log-concave. We wan't to say that there is a negative variance
|
||||
|
|
@ -339,7 +324,7 @@ class Laplace(likelihood):
|
|||
B, L, W_12 = self._compute_B_statistics(K, W)
|
||||
|
||||
W_f = np.dot(W, f)
|
||||
grad = self.likelihood_function.dlik_df(self.data, f, extra_data=self.extra_data)[:, None]
|
||||
grad = self.likelihood_function.dlik_df(self.data, f, extra_data=self.extra_data)
|
||||
#Find K_i_f
|
||||
b = W_f + grad
|
||||
|
||||
|
|
|
|||
|
|
@ -191,8 +191,8 @@ class student_t(likelihood_function):
|
|||
:returns: float(likelihood evaluated for this point)
|
||||
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
|
||||
e = y - f
|
||||
|
|
@ -207,7 +207,7 @@ class student_t(likelihood_function):
|
|||
"""
|
||||
Gradient of the link function at y, given f w.r.t f
|
||||
|
||||
$$\frac{dp(y_{i}|f_{i})}{df} = \frac{-(v+1)(f_{i}-y_{i})}{(f_{i}-y_{i})^{2} + \sigma^{2}v}$$
|
||||
$$\frac{dp(y_{i}|f_{i})}{df} = \frac{(v+1)(y_{i}-f_{i})}{(y_{i}-f_{i})^{2} + \sigma^{2}v}$$
|
||||
|
||||
:y: data
|
||||
:f: latent variables f
|
||||
|
|
@ -215,51 +215,52 @@ class student_t(likelihood_function):
|
|||
:returns: gradient of likelihood evaluated at points
|
||||
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
e = y - f
|
||||
grad = -((self.v + 1) * e) / (self.v * (self.sigma**2) + (e**2))
|
||||
return np.squeeze(grad)
|
||||
grad = ((self.v + 1) * e) / (self.v * (self.sigma**2) + (e**2))
|
||||
return grad
|
||||
|
||||
def d2lik_d2f(self, y, f, extra_data=None):
|
||||
"""
|
||||
Hessian at this point (if we are only looking at the link function not the prior) the hessian will be 0 unless i == j
|
||||
i.e. second derivative link_function at y given f f_j w.r.t f and f_j
|
||||
|
||||
Will return diagonal of hessian, since every where else it is 0
|
||||
Will return diagonal of hessian, since every where else it is 0, as the likelihood factorizes over cases
|
||||
(the distribution for y_{i} depends only on f_{i} not on f_{j!=i}
|
||||
|
||||
$$\frac{d^{2}p(y_{i}|f_{i})}{d^{3}f} = \frac{(v+1)((f_{i}-y_{i})^{2} - \sigma^{2}v)}{((f_{i}-y_{i})^{2} + \sigma^{2}v)^{2}}$$
|
||||
$$\frac{d^{2}p(y_{i}|f_{i})}{d^{3}f} = \frac{(v+1)((y_{i}-f_{i})^{2} - \sigma^{2}v)}{((y_{i}-f_{i})^{2} + \sigma^{2}v)^{2}}$$
|
||||
|
||||
:y: data
|
||||
:f: latent variables f
|
||||
:extra_data: extra_data which is not used in student t distribution
|
||||
:returns: array which is diagonal of covariance matrix (second derivative of likelihood evaluated at points)
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
|
||||
e = y - f
|
||||
hess = ((self.v + 1)*(e**2 - self.v*(self.sigma**2))) / ((((self.sigma**2)*self.v) + e**2)**2)
|
||||
return np.squeeze(hess)
|
||||
return hess
|
||||
|
||||
def d3lik_d3f(self, y, f, extra_data=None):
|
||||
"""
|
||||
Third order derivative link_function (log-likelihood ) at y given f f_j w.r.t f and f_j
|
||||
|
||||
$$\frac{d^{3}p(y_{i}|f_{i})}{d^{3}f} = \frac{-2(v+1)((f_{i} - y_{i})^3 - 3(f_{i} - y_{i}) \sigma^{2} v))}{((f_{i} - y_{i}) + \sigma^{2} v)^3}$$
|
||||
$$\frac{d^{3}p(y_{i}|f_{i})}{d^{3}f} = \frac{-2(v+1)((y_{i} - f_{i})^3 - 3(y_{i} - f_{i}) \sigma^{2} v))}{((y_{i} - f_{i}) + \sigma^{2} v)^3}$$
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
e = y - f
|
||||
d3lik_d3f = ( -(2*(self.v + 1)*(e**3 - e*3*self.v*(self.sigma**2))) /
|
||||
d3lik_d3f = ( (2*(self.v + 1)*(-e)*(e**2 - 3*self.v*(self.sigma**2))) /
|
||||
((e**2 + (self.sigma**2)*self.v)**3)
|
||||
)
|
||||
return np.squeeze(d3lik_d3f)
|
||||
return d3lik_d3f
|
||||
|
||||
def link_dstd(self, y, f, extra_data=None):
|
||||
def lik_dstd(self, y, f, extra_data=None):
|
||||
"""
|
||||
Gradient of the likelihood (lik) w.r.t sigma parameter (standard deviation)
|
||||
|
||||
|
|
@ -268,48 +269,48 @@ class student_t(likelihood_function):
|
|||
|
||||
$$\frac{dp(y_{i}|f_{i})}{d\sigma} = -\frac{1}{\sigma} + \frac{(1+v)(y_{i}-f_{i})^2}{\sigma^3 v(1 + \frac{1}{v}(\frac{(y_{i} - f_{i})}{\sigma^2})^2)}$$
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
e = y - f
|
||||
dlik_dsigma = ( (1/self.sigma) -
|
||||
((1+self.v)*(e**2))/((self.sigma**3)*self.v*(1 + (e**2) / ((self.sigma**2)*self.v) ) )
|
||||
dlik_dsigma = ( - (1/self.sigma) +
|
||||
((1+self.v)*(e**2))/((self.sigma**3)*self.v*(1 + ((e**2) / ((self.sigma**2)*self.v)) ) )
|
||||
)
|
||||
return np.squeeze(dlik_dsigma)
|
||||
return dlik_dsigma
|
||||
|
||||
def dlik_df_dstd(self, y, f, extra_data=None):
|
||||
"""
|
||||
Gradient of the dlik_df w.r.t sigma parameter (standard deviation)
|
||||
|
||||
$$\frac{d}{d\sigma}(\frac{dp(y_{i}|f_{i})}{df}) = \frac{2\sigma v(v + 1)(f-y)}{(f-y)^2 + \sigma^2 v)^2}$$
|
||||
$$\frac{d}{d\sigma}(\frac{dp(y_{i}|f_{i})}{df}) = \frac{-2\sigma v(v + 1)(y_{i}-f_{i})}{(y_{i}-f_{i})^2 + \sigma^2 v)^2}$$
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
e = y - f
|
||||
dlik_grad_dsigma = ((2*self.sigma*self.v*(self.v + 1)*e)
|
||||
dlik_grad_dsigma = ((-2*self.sigma*self.v*(self.v + 1)*e)
|
||||
/ ((self.v*(self.sigma**2) + e**2)**2)
|
||||
)
|
||||
return np.squeeze(dlik_grad_dsigma)
|
||||
return dlik_grad_dsigma
|
||||
|
||||
def d2lik_d2f_dstd(self, y, f, extra_data=None):
|
||||
"""
|
||||
Gradient of the hessian (d2lik_d2f) w.r.t sigma parameter (standard deviation)
|
||||
|
||||
$$\frac{d}{d\sigma}(\frac{d^{2}p(y_{i}|f_{i})}{d^{2}f}) = \frac{(v + 1)((f-y)^2 - \sigma^2 v)}{((f-y)^2 + \sigma^2 v)}$$
|
||||
$$\frac{d}{d\sigma}(\frac{d^{2}p(y_{i}|f_{i})}{d^{2}f}) = \frac{2\sigma v(v + 1)(\sigma^2 v - 3(y-f)^2)}{((y-f)^2 + \sigma^2 v)^3}$$
|
||||
"""
|
||||
y = np.squeeze(y)
|
||||
f = np.squeeze(f)
|
||||
#y = np.squeeze(y)
|
||||
#f = np.squeeze(f)
|
||||
assert y.shape == f.shape
|
||||
e = y - f
|
||||
dlik_hess_dsigma = ( ((self.v + 1)*(e**2 - (self.sigma**2)*self.v)) /
|
||||
((e**2 + (self.sigma**2)*self.v)**2)
|
||||
dlik_hess_dsigma = ( (2*self.sigma*self.v*(self.v + 1)*((self.sigma**2)*self.v - 3*(e**2))) /
|
||||
((e**2 + (self.sigma**2)*self.v)**3)
|
||||
)
|
||||
return dlik_hess_dsigma
|
||||
|
||||
def _gradients(self, y, f, extra_data=None):
|
||||
#must be listed in same order as 'get_param_names'
|
||||
derivs = ([self.link_dstd(y, f, extra_data=extra_data)],
|
||||
derivs = ([self.lik_dstd(y, f, extra_data=extra_data)],
|
||||
[self.dlik_df_dstd(y, f, extra_data=extra_data)],
|
||||
[self.d2lik_d2f_dstd(y, f, extra_data=extra_data)]
|
||||
) # lists as we might learn many parameters
|
||||
|
|
|
|||
|
|
@ -142,13 +142,22 @@ class GP(model):
|
|||
Note, we use the chain rule: dL_dtheta = dL_dK * d_K_dtheta
|
||||
"""
|
||||
dL_dthetaK = self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X)
|
||||
print "dL_dthetaK before: ",dL_dthetaK
|
||||
if isinstance(self.likelihood, Laplace):
|
||||
#Reapproximate incase it hasnt been done...
|
||||
if isinstance(self.likelihood, Laplace):
|
||||
self.likelihood.fit_full(self.kern.K(self.X))
|
||||
self.likelihood._set_params(self.likelihood._get_params())
|
||||
|
||||
#Need to pass in a matrix of ones to get access to raw dK_dthetaK values without being chained
|
||||
fake_dL_dKs = np.ones(self.dL_dK.shape) #FIXME: Check this is right...
|
||||
#fake_dL_dKs = np.eye(self.dL_dK.shape[0]) #FIXME: Check this is right...
|
||||
dK_dthetaK = self.kern.dK_dtheta(dL_dK=fake_dL_dKs, X=self.X)
|
||||
#THIS SHOULD NOT BE (1,num_k_params) matrix it should be (N,N,num_k_params)
|
||||
|
||||
dL_dthetaK = self.likelihood._Kgradients(dL_d_K_Sigma=self.dL_dK, dK_dthetaK=dK_dthetaK)
|
||||
dL_dthetaL = self.likelihood._gradients(partial=np.diag(self.dL_dK))
|
||||
dL_dthetaK = self.likelihood._Kgradients(dK_dthetaK=dK_dthetaK)
|
||||
dL_dthetaL = 0 # self.likelihood._gradients(partial=np.diag(self.dL_dK))
|
||||
print "dL_dthetaK after: ",dL_dthetaK
|
||||
#print "Stacked dL_dthetaK, dL_dthetaL: ", np.hstack((dL_dthetaK, dL_dthetaL))
|
||||
else:
|
||||
dL_dthetaL = self.likelihood._gradients(partial=np.diag(self.dL_dK))
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue