kernel.py 13.5 KB
Newer Older
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
1
2
3
4
5
import numpy as np
from abc import ABC, abstractmethod

from scipy.spatial.distance import pdist, cdist, squareform

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
6
7
from time import time

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
8
9
10
11
12
class kernel(ABC):
    def __init__(self):
        self._theta = None

    @abstractmethod
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
13
    def __call__(self):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
14
15
        pass

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
16
17
18
    @abstractmethod
    def kernel(self):
        pass
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
19

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
20
21
22
    @abstractmethod
    def kernel_vector(self):
        pass
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

    @abstractmethod
    def kernel_jacobian(self):
        pass

    @abstractmethod
    def kernel_hyperparameter_gradient(self):
        pass

    @property
    def theta(self):
        return self._theta

    @theta.setter
    def theta(self, theta):
        self._theta = theta
    
    def numerical_jacobian(self,x,y, dx=1.e-5):
        if np.ndim(y) == 1:
            y = y.reshape((1,-1))
        nx = len(x)
        ny = y.shape[0]
        f_jac = np.zeros((ny,nx))
        for i in range(nx):
            x_up = np.copy(x)
            x_down = np.copy(x)
            x_up[i] += 0.5*dx
            x_down[i] -= 0.5*dx
            
            f_up = self.kernel(x_up,y)
            f_down = self.kernel(x_down,y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
54

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
55
56
57
58
            f_jac[:,i] = (f_up - f_down)/dx
        return f_jac

    def numerical_hyperparameter_gradient(self,X, dx=1.e-5):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
59
60
61
        """Calculates the numerical derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
62
63
64
65
66
67
68
69
70
71
72
        N_data = X.shape[0]
        theta = np.copy(self.theta)
        N_hyper = len(theta)
        dK_dTheta = np.zeros((N_hyper, N_data, N_data))
        for i in range(N_hyper):
            theta_up = np.copy(theta)
            theta_down = np.copy(theta)
            theta_up[i] += 0.5*dx
            theta_down[i] -= 0.5*dx
            
            self.theta = theta_up
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
73
            K_up = self(X, eval_gradient=False)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
74
            self.theta = theta_down
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
75
            K_down = self(X, eval_gradient=False)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
76

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
77
            dK_dTheta[i,:,:] = (K_up - K_down)/dx
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
78
        self.theta = theta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
79
80
81
82
        return dK_dTheta


class gauss_kernel(kernel):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
83
84
    def __init__(self, amplitude=100.0, amplitude_bounds=(1e0, 1e5),
                 length_scale=10.0, length_scale_bounds=(1e0, 1e3),
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
85
                 noise=1e-5, noise_bounds=None,
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
86
87
                 eta=1, eta_bounds=(0.1,10),
                 Nsplit_eta=None):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
88
89
        self.amplitude = amplitude
        self.length_scale = length_scale
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
90

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
91
92
93
        self.amplitude_bounds = amplitude_bounds
        self.length_scale_bounds = length_scale_bounds

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
94
95
        self.theta_bounds = np.array([amplitude_bounds, length_scale_bounds])

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
96
97
98
99
100
101
102
103
        self.amplitude = amplitude
        self.length_scale = length_scale
        self.noise = noise
        self.eta = eta
        self.Nsplit_eta = Nsplit_eta

        self.amplitude_bounds = amplitude_bounds
        self.length_scale_bounds = length_scale_bounds
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
104
        self.noise_bounds = (noise,noise)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
105
106
107
108
109
        if self.Nsplit_eta is None:
            self.eta_bounds = (eta,eta)
        else:
            self.eta_bounds = eta_bounds

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
110
        self.theta_bounds = np.log(np.array([amplitude_bounds, length_scale_bounds, self.noise_bounds, self.eta_bounds]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
111

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
112
113
114
115
116
117
118
119
120
121
122
123
    def __call__(self, X, eval_gradient=False):
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        d = cdist(X / self.length_scale,
                  X / self.length_scale, metric='sqeuclidean')
        K = self.amplitude * np.exp(-0.5 * d)
        
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
124
125
126
127
128
129
130
131
132
133
134
135

    def __call__(self, X, eval_gradient=False):
        K = self.kernel(X, with_noise=True)
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K

    def kernel(self, X, Y=None, with_noise=False):
        if with_noise:
            assert Y is None
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
136
137
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
138
139
140
        if Y is None:
            Y = X
        elif np.ndim(Y) == 1:
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
141
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
142
143
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
144
145
        d = cdist(X / self.length_scale,
                  Y / self.length_scale, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
146
147
148
149
        if with_noise:
            K = self.amplitude * (np.exp(-0.5 * d) + self.noise*np.eye(X.shape[0]))
        else:
            K = self.amplitude * np.exp(-0.5 * d)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
150
151
        return K

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
152
153
154
155
    def kernel_value(self, x,y):
        K = self.kernel(x,y)
        return np.asscalar(K)
    
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
156
    def kernel_vector(self, x,Y):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
157
        K = self.kernel(x,Y).reshape(-1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
158
159
160
        return K

    def kernel_jacobian(self, X,Y):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
161
162
163
164
        """ Jacobian of the kernel with respect to X
        """
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
165
166
167
168
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        if np.ndim(Y) == 1:
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
169
170
171
172
173
        d = cdist(X / self.length_scale,
                  Y / self.length_scale, metric='sqeuclidean')
        dK_dd = -self.amplitude * 1/(2*self.length_scale**2) * np.exp(-0.5 * d)
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
174
175
        dd_df = 2*(X - Y)

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
176
        dk_df = np.multiply(dK_dd.T, dd_df)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
177
178
179
180
        return dk_df

    @property
    def theta(self):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
181
182
        """Returns the log-transformed hyperparameters of the kernel.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
183
        self._theta = np.array([self.amplitude, self.length_scale, self.noise, self.eta]) 
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
184
        return np.log(self._theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
185
186
187

    @theta.setter
    def theta(self, theta):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
188
189
190
191
192
        """Sets the hyperparameters of the kernel.

        theta: log-transformed hyperparameters
        """
        self._theta = np.exp(theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
193
194
        self.amplitude = self._theta[0]
        self.length_scale = self._theta[1]
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
195
196
197
198
199
200
201
202
203
204
205
206
        self.noise = self._theta[2]
        self.eta = self._theta[3]

    def apply_eta(self, X):
        Xeta = np.copy(X)
        if self.Nsplit_eta is not None:
            if np.ndim(X) == 1:
                Xeta[self.Nsplit_eta:] *= self.eta
            else:
                Xeta[:,self.Nsplit_eta:] *= self.eta
        return Xeta

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
207
208
209
    def dK_da(self, X):
        d = cdist(X / self.length_scale,
                  X / self.length_scale, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
210
        dK_da = self.amplitude * (np.exp(-0.5 * d) + self.noise*np.eye(X.shape[0]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
211
212
213
214
215
        return dK_da
        
    def dK_dl(self, X):
        d = cdist(X / self.length_scale,
                  X / self.length_scale, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
216
        dK_dl = self.amplitude * d * np.exp(-0.5 * d)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
217
218
        return dK_dl

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
    def dK_dn(self, X):
        dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
        return dK_dn

    def dK_deta(self, X, dx=1e-5):
        N_data = X.shape[0]
        theta = np.copy(self.theta)
        dK_deta = np.zeros((N_data, N_data))

        theta_up = np.copy(theta)
        theta_down = np.copy(theta)
        theta_up[-1] += 0.5*dx
        theta_down[-1] -= 0.5*dx
        
        self.theta = theta_up
        K_up = self(X, eval_gradient=False)
        self.theta = theta_down
        K_down = self(X, eval_gradient=False)
        dK_dTheta = (K_up - K_down)/dx

        self.theta = theta
        return dK_dTheta

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
242
    def kernel_hyperparameter_gradient(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
243
244
245
        """Calculates the derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
246
247
248
249
        dK_deta = self.dK_deta(X)
        X = self.apply_eta(X)
        return np.array([self.dK_da(X), self.dK_dl(X), self.dK_dn(X), dK_deta])

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
250
251
252
    

class double_gauss_kernel(kernel):
253
    def __init__(self, amplitude=100., amplitude_bounds=(1e1,1e5),
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
254
255
                 length_scale1=10.0, length_scale1_bounds=(1e0, 1e3),
                 length_scale2=10.0, length_scale2_bounds=(1e0, 1e3),
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
256
257
                 weight=0.01, weight_bounds=None,
                 noise=1e-5, noise_bounds=None,
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
258
259
                 eta=1, eta_bounds=(0.1,10),
                 Nsplit_eta=None):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
260
261
262
263
        self.amplitude = amplitude
        self.length_scale1 = length_scale1
        self.length_scale2 = length_scale2
        self.weight = weight
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
264
        self.noise = noise
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
265
266
        self.eta = eta
        self.Nsplit_eta = Nsplit_eta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
267
268
269
270

        self.amplitude_bounds = amplitude_bounds
        self.length_scale1_bounds = length_scale1_bounds
        self.length_scale2_bounds = length_scale2_bounds
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
271
272
        self.weight_bounds = (weight, weight)
        self.noise_bounds = (noise, noise)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
273
274
275
276
        if self.Nsplit_eta is None:
            self.eta_bounds = (eta,eta)
        else:
            self.eta_bounds = eta_bounds
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
277

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
278
        self.theta_bounds = np.log(np.array([amplitude_bounds, length_scale1_bounds, length_scale2_bounds, self.weight_bounds, self.noise_bounds, self.eta_bounds]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
279

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
280
    def __call__(self, X, eval_gradient=False):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
281
        K = self.kernel(X, with_noise=True)
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
282
283
284
285
286
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K
287

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
288
289
290
    def kernel(self, X, Y=None, with_noise=False):
        if with_noise:
            assert Y is None
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
291
292
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
293
294
295
        if Y is None:
            Y = X
        elif np.ndim(Y) == 1:
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
296
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
297
298
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
299
300
301
302
        d1 = cdist(X / self.length_scale1,
                  Y / self.length_scale1, metric='sqeuclidean')
        d2 = cdist(X / self.length_scale2,
                  Y / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
303
304
305
306
        if with_noise:
            K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2) + self.noise*np.eye(X.shape[0]))
        else:
            K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
307
308
        return K

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
309
310
311
312
    def kernel_value(self, x,y):
        K = self.kernel(x,y)
        return np.asscalar(K)
    
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
313
    def kernel_vector(self, x,Y):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
314
        K = self.kernel(x,Y).reshape(-1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
315
316
317
318
319
        return K

    def kernel_jacobian(self, X,Y):
        """ Jacobian of the kernel with respect to X
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
320
321
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
322
323
324
325
326
327
328
329
330
331
332
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        if np.ndim(Y) == 1:
            Y = Y.reshape((1,-1))
        d1 = cdist(X / self.length_scale1,
                   Y / self.length_scale1, metric='sqeuclidean')
        d2 = cdist(X / self.length_scale2,
                   Y / self.length_scale2, metric='sqeuclidean')
        dK1_dd1 = -1/(2*self.length_scale1**2) * np.exp(-0.5 * d1)
        dK2_dd2 = -1/(2*self.length_scale2**2) * np.exp(-0.5 * d2)
        dK_dd = self.amplitude * (dK1_dd1 + self.weight*dK2_dd2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
333
334
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
335
336
337
338
339
340
341
        dd_df = 2*(X - Y)

        dk_df = np.multiply(dK_dd.T, dd_df)
        return dk_df

    @property
    def theta(self):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
342
343
        """Returns the log-transformed hyperparameters of the kernel.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
344
        self._theta = np.array([self.amplitude, self.length_scale1, self.length_scale2, self.weight, self.noise, self.eta])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
345
        return np.log(self._theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
346
347
348

    @theta.setter
    def theta(self, theta):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
349
350
351
352
353
        """Sets the hyperparameters of the kernel.

        theta: log-transformed hyperparameters
        """
        self._theta = np.exp(theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
354
355
356
357
        self.amplitude = self._theta[0]
        self.length_scale1 = self._theta[1]
        self.length_scale2 = self._theta[2]
        self.weight = self._theta[3]
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
358
        self.noise = self._theta[4]
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
359
360
361
362
363
364
365
366
367
368
        self.eta = self._theta[5]

    def apply_eta(self, X):
        Xeta = np.copy(X)
        if self.Nsplit_eta is not None:
            if np.ndim(X) == 1:
                Xeta[self.Nsplit_eta:] *= self.eta
            else:
                Xeta[:,self.Nsplit_eta:] *= self.eta
        return Xeta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
369
370
371
372
373
374

    def dK_da(self, X):
        d1 = cdist(X / self.length_scale1,
                   X / self.length_scale1, metric='sqeuclidean')
        d2 = cdist(X / self.length_scale2,
                   X / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
375
        dK_da = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2) + self.noise*np.eye(X.shape[0]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
376
377
378
379
380
        return dK_da
        
    def dK_dl1(self, X):
        d1 = cdist(X / self.length_scale1,
                   X / self.length_scale1, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
381
        dK_dl1 = self.amplitude*d1 * np.exp(-0.5 * d1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
382
383
384
385
386
        return dK_dl1

    def dK_dl2(self, X):
        d2 = cdist(X / self.length_scale2,
                   X / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
387
        dK_dl2 = self.amplitude*self.weight*d2 * np.exp(-0.5 * d2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
388
389
390
391
392
        return dK_dl2

    def dK_dw(self, X):
        d2 = cdist(X / self.length_scale2,
                   X / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
393
        dK_dl2 = self.amplitude*self.weight*np.exp(-0.5 * d2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
394
395
        return dK_dl2

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
396
    def dK_dn(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
397
        dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
398
399
        return dK_dn

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
    def dK_deta(self, X, dx=1e-5):
        N_data = X.shape[0]
        theta = np.copy(self.theta)
        dK_deta = np.zeros((N_data, N_data))

        theta_up = np.copy(theta)
        theta_down = np.copy(theta)
        theta_up[-1] += 0.5*dx
        theta_down[-1] -= 0.5*dx
        
        self.theta = theta_up
        K_up = self(X, eval_gradient=False)
        self.theta = theta_down
        K_down = self(X, eval_gradient=False)
        dK_dTheta = (K_up - K_down)/dx

        self.theta = theta
        return dK_dTheta

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
419
    def kernel_hyperparameter_gradient(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
420
421
422
        """Calculates the derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
423
        dK_deta = self.dK_deta(X)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
424
        X = self.apply_eta(X)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
425
        return np.array([self.dK_da(X), self.dK_dl1(X), self.dK_dl2(X), self.dK_dw(X), self.dK_dn(X), dK_deta])