kernel.py 14.6 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)

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

179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    def kernel_hessian(self, X,Y):
        """ Jacobian of the kernel with respect to X
        """
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        if np.ndim(Y) == 1:
            Y = Y.reshape((1,-1))
        d = cdist(X / self.length_scale,
                  Y / self.length_scale, metric='sqeuclidean')
        u = 1/(2*self.length_scale**2)
        K = self.amplitude * np.exp(-0.5 * d)
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)

        N_X, Nf = X.shape
        N_Y = Y.shape[0]

        X_rep = np.tile(X, (N_Y,1,1)).swapaxes(0,1)
        Y_rep = np.tile(Y, (N_X,1,1))

        dd_dX = -2*(Y_rep-X_rep)  # shape: (N_X, N_Y, Nf)
        dd_dY = -dd_dX
        dd_dX_dd_dY = np.einsum('nmi,nmj->nmij',dd_dX, dd_dY)  # shape: (N_X, N_Y, Nf, Nf)
        if self.Nsplit_eta is not None:
            diag = np.diag(self.Nsplit_eta*[1]+(Nf-self.Nsplit_eta)*[self.eta**2])
        else:
            diag = np.identity(Nf)
        d2d_dXdY = np.tile(-2*diag, (N_X, N_Y,1,1))  # shape: (N_X, N_Y, Nf, Nf)
        hess = -u*K.reshape(N_X,N_Y,1,1) * (u*dd_dX_dd_dY - d2d_dXdY)

        return hess

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
213
214
    @property
    def theta(self):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
215
216
        """Returns the log-transformed hyperparameters of the kernel.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
217
        self._theta = np.array([self.amplitude, self.length_scale, self.noise, self.eta]) 
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
218
        return np.log(self._theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
219
220
221

    @theta.setter
    def theta(self, theta):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
222
223
224
225
226
        """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
227
228
        self.amplitude = self._theta[0]
        self.length_scale = self._theta[1]
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
229
230
231
232
233
234
235
236
237
238
239
240
        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
241
242
243
    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
244
        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
245
246
247
248
249
        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
250
        dK_dl = self.amplitude * d * np.exp(-0.5 * d)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
251
252
        return dK_dl

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
    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
276
    def kernel_hyperparameter_gradient(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
277
278
279
        """Calculates the derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
280
281
282
283
        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
284
285
286
    

class double_gauss_kernel(kernel):
287
    def __init__(self, amplitude=100., amplitude_bounds=(1e1,1e5),
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
288
289
                 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
290
291
                 weight=0.01, weight_bounds=None,
                 noise=1e-5, noise_bounds=None,
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
292
293
                 eta=1, eta_bounds=(0.1,10),
                 Nsplit_eta=None):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
294
295
296
297
        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
298
        self.noise = noise
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
299
300
        self.eta = eta
        self.Nsplit_eta = Nsplit_eta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
301
302
303
304

        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
305
306
        self.weight_bounds = (weight, weight)
        self.noise_bounds = (noise, noise)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
307
308
309
310
        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
311

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
312
        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
313

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
314
    def __call__(self, X, eval_gradient=False):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
315
        K = self.kernel(X, with_noise=True)
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
316
317
318
319
320
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K
321

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
322
323
324
    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
325
326
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
327
328
329
        if Y is None:
            Y = X
        elif np.ndim(Y) == 1:
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
330
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
331
332
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
333
334
335
336
        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
337
338
339
340
        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
341
342
        return K

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
343
344
345
346
    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
347
    def kernel_vector(self, x,Y):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
348
        K = self.kernel(x,Y).reshape(-1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
349
350
351
352
353
        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
354
355
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
356
357
358
359
360
361
362
363
364
365
366
        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
367
368
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
369
370
371
372
373
374
375
        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
376
377
        """Returns the log-transformed hyperparameters of the kernel.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
378
        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
379
        return np.log(self._theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
380
381
382

    @theta.setter
    def theta(self, theta):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
383
384
385
386
387
        """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
388
389
390
391
        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
392
        self.noise = self._theta[4]
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
393
394
395
396
397
398
399
400
401
402
        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
403
404
405
406
407
408

    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
409
        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
410
411
412
413
414
        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
415
        dK_dl1 = self.amplitude*d1 * np.exp(-0.5 * d1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
416
417
418
419
420
        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
421
        dK_dl2 = self.amplitude*self.weight*d2 * np.exp(-0.5 * d2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
422
423
424
425
426
        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
427
        dK_dl2 = self.amplitude*self.weight*np.exp(-0.5 * d2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
428
429
        return dK_dl2

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
430
    def dK_dn(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
431
        dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
432
433
        return dK_dn

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
    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
453
    def kernel_hyperparameter_gradient(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
454
455
456
        """Calculates the derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
457
        dK_deta = self.dK_deta(X)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
458
        X = self.apply_eta(X)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
459
        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])