mvpa.clfs.kernel

1 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 2 # vi: set ft=python sts=4 ts=4 sw=4 et: 3 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 4 # 5 # Copyright (c) 2008 Emanuele Olivetti <emanuele@relativita.com> 6 # See COPYING file distributed along with the PyMVPA package for the 7 # copyright and license terms. 8 # 9 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 10 """Kernels for Gaussian Process Regression and Classification.""" 11 12 13 _DEV__DOC__ = """ 14 Make use of Parameter Collections to keep parameters of the 15 kernels. Then we would get a uniform .reset() functionality. Now reset 16 is provided just for parts which are failing in the unittests, but 17 there is many more places where they are not reset properly if 18 classifier gets trained on some new data of different dimensionality 19 """ 20 21 __docformat__ = 'restructuredtext' 22 23 24 import numpy as N 25 26 from mvpa.misc.exceptions import InvalidHyperparameterError 27 from mvpa.clfs.distance import squared_euclidean_distance 28 29 if __debug__: 30 from mvpa.base import debug, warning 31 32

33 -class Kernel(object):

34 """Kernel function base class. 35 36 """ 37

38 - def __init__(self):

39 pass

40

41 - def __repr__(self):

42 return "Kernel()"

43

44 - def compute(self, data1, data2=None):

45 raise NotImplementedError

46

47 - def reset(self):

48 """Resets the kernel dropping internal variables to the original values""" 49 pass

50

51 - def compute_gradient(self,alphaalphaTK):

52 raise NotImplementedError

53

54 - def compute_lml_gradient(self,alphaalphaT_Kinv,data):

55 raise NotImplementedError

56

57 - def compute_lml_gradient_logscale(self,alphaalphaT_Kinv,data):

58 raise NotImplementedError

59 60 pass

61

62 -class KernelConstant(Kernel):

63 """The constant kernel class. 64 """

65 - def __init__(self, sigma_0=1.0, **kwargs):

66 """Initialize the constant kernel instance. 67 68 :Parameters: 69 sigma_0 : float 70 standard deviation of the Gaussian prior probability 71 N(0,sigma_0**2) of the intercept of the constant regression. 72 (Defaults to 1.0) 73 """ 74 # init base class first 75 Kernel.__init__(self, **kwargs) 76 77 self.sigma_0 = sigma_0 78 self.kernel_matrix = None

79

80 - def __repr__(self):

81 return "%s(sigma_0=%s)" % (self.__class__.__name__, str(self.sigma_0))

82

83 - def compute(self, data1, data2=None):

84 """Compute kernel matrix. 85 86 :Parameters: 87 data1 : numpy.ndarray 88 data 89 data2 : numpy.ndarray 90 data 91 (Defaults to None) 92 """ 93 if data2 is None: 94 data2 = data1 95 pass 96 self.kernel_matrix = \ 97 (self.sigma_0 ** 2) * N.ones((data1.shape[0], data2.shape[0])) 98 return self.kernel_matrix

99

100 - def set_hyperparameters(self, hyperparameter):

101 if hyperparameter < 0: 102 raise InvalidHyperparameterError() 103 self.sigma_0 = hyperparameter 104 return

105

106 - def compute_lml_gradient(self,alphaalphaT_Kinv,data):

107 K_grad_sigma_0 = 2*self.sigma_0 108 # self.lml_gradient = 0.5*(N.trace(N.dot(alphaalphaT_Kinv,K_grad_sigma_0*N.ones(alphaalphaT_Kinv.shape))) 109 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 110 # Fastest when B is a constant: B*A.sum() 111 self.lml_gradient = 0.5*N.array(K_grad_sigma_0*alphaalphaT_Kinv.sum()) 112 return self.lml_gradient

113

114 - def compute_lml_gradient_logscale(self,alphaalphaT_Kinv,data):

115 K_grad_sigma_0 = 2*self.sigma_0**2 116 self.lml_gradient = 0.5*N.array(K_grad_sigma_0*alphaalphaT_Kinv.sum()) 117 return self.lml_gradient

118 119 pass

120 121

122 -class KernelLinear(Kernel):

123 """The linear kernel class. 124 """

125 - def __init__(self, Sigma_p=None, sigma_0=1.0, **kwargs):

126 """Initialize the linear kernel instance. 127 128 :Parameters: 129 Sigma_p : numpy.ndarray 130 Covariance matrix of the Gaussian prior probability N(0,Sigma_p) 131 on the weights of the linear regression. 132 (Defaults to None) 133 sigma_0 : float 134 the standard deviation of the Gaussian prior N(0,sigma_0**2) 135 of the intercept of the linear regression. 136 (Deafults to 1.0) 137 """ 138 # init base class first 139 Kernel.__init__(self, **kwargs) 140 141 # TODO: figure out cleaner way... probably by using KernelParameters ;-) 142 self.Sigma_p = Sigma_p 143 self.sigma_0 = sigma_0 144 self.kernel_matrix = None

145 146

147 - def __repr__(self):

148 return "%s(Sigma_p=%s, sigma_0=%s)" \ 149 % (self.__class__.__name__, str(self.Sigma_p), str(self.sigma_0))

150 151

152 - def reset(self):

153 super(KernelLinear, self).reset() 154 self._Sigma_p = self._Sigma_p_orig

155 156

157 - def compute(self, data1, data2=None):

158 """Compute kernel matrix. 159 Set Sigma_p to correct dimensions and default value if necessary. 160 161 :Parameters: 162 data1 : numpy.ndarray 163 data 164 data2 : numpy.ndarray 165 data 166 (Defaults to None) 167 """ 168 if data2 is None: 169 data2 = data1 170 pass 171 172 # it is better to use separate lines of computation, to don't 173 # incure computation cost without need (otherwise 174 # N.dot(self.Sigma_p, data2.T) can take forever for relatively 175 # large number of features) 176 177 Sigma_p = self.Sigma_p # local binding 178 #if scalar - scale second term appropriately 179 if N.isscalar(Sigma_p): 180 if Sigma_p == 1.0: 181 data2_sc = data2.T 182 else: 183 data2_sc = Sigma_p * data2.T 184 185 # if vector use it as diagonal matrix -- ie scale each row by 186 # the given value 187 elif len(Sigma_p.shape) == 1 and \ 188 Sigma_p.shape[0] == data1.shape[1]: 189 # which due to numpy broadcasting is the same as product 190 # with scalar above 191 data2_sc = (Sigma_p * data1).T 192 193 # if it is a full matrix -- full-featured and lengthy 194 # matrix product 195 else: 196 data2_sc = N.dot(Sigma_p, data2.T) 197 pass 198 199 # XXX if Sigma_p is changed a warning should be issued! 200 # XXX other cases of incorrect Sigma_p could be catched 201 self.kernel_matrix = N.dot(data1, data2_sc) + self.sigma_0 ** 2 202 return self.kernel_matrix

203

204 - def set_hyperparameters(self, hyperparameter):

205 # XXX in the next line we assume that the values we want to 206 # assign to Sigma_p are a constant or a vector (the diagonal 207 # of Sigma_p actually). This is a limitation since these 208 # values could be in general an hermitian matrix (i.e., a 209 # covariance matrix)... but how to tell ModelSelector/OpenOpt 210 # to proved just "hermitian" set of values? So for now we skip 211 # the general case, which seems not to useful indeed. 212 if N.any(hyperparameter < 0): 213 raise InvalidHyperparameterError() 214 self.sigma_0 = N.array(hyperparameter[0]) 215 self._Sigma_p = N.diagflat(hyperparameter[1:]) 216 return

217

218 - def compute_lml_gradient(self,alphaalphaT_Kinv,data):

219 def lml_grad(K_grad_i): 220 # return N.trace(N.dot(alphaalphaT_Kinv,K_grad_i)) 221 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 222 return (alphaalphaT_Kinv*(K_grad_i.T)).sum()

223 lml_gradient = [] 224 lml_gradient.append(2*self.sigma_0*alphaalphaT_Kinv.sum()) 225 for i in range(self.Sigma_p.shape[0]): 226 # Note that Sigma_p is not squared in compute() so it 227 # disappears in the partial derivative: 228 K_grad_i = N.multiply.outer(data[:,i],data[:,i]) 229 lml_gradient.append(lml_grad(K_grad_i)) 230 pass 231 self.lml_gradient = lml_gradient = 0.5*N.array(lml_gradient) 232 return lml_gradient

233

234 - def compute_lml_gradient_logscale(self,alphaalphaT_Kinv,data):

235 def lml_grad(K_grad_i): 236 # return N.trace(N.dot(alphaalphaT_Kinv,K_grad_i)) 237 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 238 return (alphaalphaT_Kinv*(K_grad_i.T)).sum()

239 lml_gradient = [] 240 lml_gradient.append(2*self.sigma_0**2*alphaalphaT_Kinv.sum()) 241 Sigma_p = self.Sigma_p # local binding 242 for i in range(Sigma_p.shape[0]): 243 # Note that Sigma_p is not squared in compute() so it 244 # disappears in the partial derivative: 245 K_grad_log_i = Sigma_p[i,i]*N.multiply.outer(data[:,i],data[:,i]) 246 lml_gradient.append(lml_grad(K_grad_log_i)) 247 pass 248 self.lml_gradient = lml_gradient = 0.5*N.array(lml_gradient) 249 return lml_gradient 250 251

252 - def _setSigma_p(self, v):

253 """Set Sigma_p value and store _orig for reset 254 """ 255 if (v is None): 256 # if Sigma_p is not set use a scalar 1.0 257 v = 1.0 258 self._Sigma_p_orig = self._Sigma_p = v

259 260 Sigma_p = property(fget=lambda x:x._Sigma_p, 261 fset=_setSigma_p) 262 263 pass 264 265

266 -class KernelExponential(Kernel):

267 """The Exponential kernel class. 268 269 Note that it can handle a length scale for each dimension for 270 Automtic Relevance Determination. 271 272 """

273 - def __init__(self, length_scale=1.0, sigma_f = 1.0, **kwargs):

274 """Initialize an Exponential kernel instance. 275 276 :Parameters: 277 length_scale : float OR numpy.ndarray 278 the characteristic length-scale (or length-scales) of the 279 phenomenon under investigation. 280 (Defaults to 1.0) 281 sigma_f : float 282 Signal standard deviation. 283 (Defaults to 1.0) 284 """ 285 # init base class first 286 Kernel.__init__(self, **kwargs) 287 288 self.length_scale = length_scale 289 self.sigma_f = sigma_f 290 self.kernel_matrix = None

291 292

293 - def __repr__(self):

294 return "%s(length_scale=%s, sigma_f=%s)" \ 295 % (self.__class__.__name__, str(self.length_scale), str(self.sigma_f))

296

297 - def compute(self, data1, data2=None):

298 """Compute kernel matrix. 299 300 :Parameters: 301 data1 : numpy.ndarray 302 data 303 data2 : numpy.ndarray 304 data 305 (Defaults to None) 306 """ 307 # XXX the following computation can be (maybe) made more 308 # efficient since length_scale is squared and then 309 # square-rooted uselessly. 310 # Weighted euclidean distance matrix: 311 self.wdm = N.sqrt(squared_euclidean_distance( 312 data1, data2, weight=(self.length_scale**-2))) 313 self.kernel_matrix = \ 314 self.sigma_f**2 * N.exp(-self.wdm) 315 return self.kernel_matrix

316

317 - def gradient(self, data1, data2):

318 """Compute gradient of the kernel matrix. A must for fast 319 model selection with high-dimensional data. 320 """ 321 raise NotImplementedError

322

323 - def set_hyperparameters(self, hyperparameter):

324 """Set hyperaparmeters from a vector. 325 326 Used by model selection. 327 """ 328 if N.any(hyperparameter < 0): 329 raise InvalidHyperparameterError() 330 self.sigma_f = hyperparameter[0] 331 self.length_scale = hyperparameter[1:] 332 return

333

334 - def compute_lml_gradient(self,alphaalphaT_Kinv,data):

335 """Compute grandient of the kernel and return the portion of 336 log marginal likelihood gradient due to the kernel. 337 Shorter formula. Allows vector of lengthscales (ARD) 338 BUT THIS LAST OPTION SEEMS NOT TO WORK FOR (CURRENTLY) 339 UNKNOWN REASONS. 340 """ 341 self.lml_gradient = [] 342 def lml_grad(K_grad_i): 343 # return N.trace(N.dot(alphaalphaT_Kinv,K_grad_i)) 344 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 345 return (alphaalphaT_Kinv*(K_grad_i.T)).sum()

346 grad_sigma_f = 2.0/self.sigma_f*self.kernel_matrix 347 self.lml_gradient.append(lml_grad(grad_sigma_f)) 348 if N.isscalar(self.length_scale) or self.length_scale.size==1: 349 # use the same length_scale for all dimensions: 350 K_grad_l = self.wdm*self.kernel_matrix*(self.length_scale**-1) 351 self.lml_gradient.append(lml_grad(K_grad_l)) 352 else: 353 # use one length_scale for each dimension: 354 for i in range(self.length_scale.size): 355 K_grad_i = (self.length_scale[i]**-3)*(self.wdm**-1)*self.kernel_matrix*N.subtract.outer(data[:,i],data[:,i])**2 356 self.lml_gradient.append(lml_grad(K_grad_i)) 357 pass 358 pass 359 self.lml_gradient = 0.5*N.array(self.lml_gradient) 360 return self.lml_gradient

361

362 - def compute_lml_gradient_logscale(self,alphaalphaT_Kinv,data):

363 """Compute grandient of the kernel and return the portion of 364 log marginal likelihood gradient due to the kernel. 365 Shorter formula. Allows vector of lengthscales (ARD). 366 BUT THIS LAST OPTION SEEMS NOT TO WORK FOR (CURRENTLY) 367 UNKNOWN REASONS. 368 """ 369 self.lml_gradient = [] 370 def lml_grad(K_grad_i): 371 # return N.trace(N.dot(alphaalphaT_Kinv,K_grad_i)) 372 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 373 return (alphaalphaT_Kinv*(K_grad_i.T)).sum()

374 grad_log_sigma_f = 2.0*self.kernel_matrix 375 self.lml_gradient.append(lml_grad(grad_log_sigma_f)) 376 if N.isscalar(self.length_scale) or self.length_scale.size==1: 377 # use the same length_scale for all dimensions: 378 K_grad_l = self.wdm*self.kernel_matrix 379 self.lml_gradient.append(lml_grad(K_grad_l)) 380 else: 381 # use one length_scale for each dimension: 382 for i in range(self.length_scale.size): 383 K_grad_i = (self.length_scale[i]**-2)*(self.wdm**-1)*self.kernel_matrix*N.subtract.outer(data[:,i],data[:,i])**2 384 self.lml_gradient.append(lml_grad(K_grad_i)) 385 pass 386 pass 387 self.lml_gradient = 0.5*N.array(self.lml_gradient) 388 return self.lml_gradient 389 390 pass 391 392

393 -class KernelSquaredExponential(Kernel):

394 """The Squared Exponential kernel class. 395 396 Note that it can handle a length scale for each dimension for 397 Automtic Relevance Determination. 398 399 """

400 - def __init__(self, length_scale=1.0, sigma_f=1.0, **kwargs):

401 """Initialize a Squared Exponential kernel instance. 402 403 :Parameters: 404 length_scale : float OR numpy.ndarray 405 the characteristic length-scale (or length-scales) of the 406 phenomenon under investigation. 407 (Defaults to 1.0) 408 sigma_f : float 409 Signal standard deviation. 410 (Defaults to 1.0) 411 """ 412 # init base class first 413 Kernel.__init__(self, **kwargs) 414 415 self.length_scale = length_scale 416 self.sigma_f = sigma_f 417 self.kernel_matrix = None

418 419

420 - def reset(self):

421 super(KernelSquaredExponential, self).reset() 422 self._length_scale = self._length_scale_orig

423 424

425 - def __repr__(self):

426 return "%s(length_scale=%s, sigma_f=%s)" \ 427 % (self.__class__.__name__, str(self.length_scale), str(self.sigma_f))

428

429 - def compute(self, data1, data2=None):

430 """Compute kernel matrix. 431 432 :Parameters: 433 data1 : numpy.ndarray 434 data 435 data2 : numpy.ndarray 436 data 437 (Defaults to None) 438 """ 439 # weighted squared euclidean distance matrix: 440 self.wdm2 = squared_euclidean_distance(data1, data2, weight=(self.length_scale**-2)) 441 self.kernel_matrix = self.sigma_f**2 * N.exp(-0.5*self.wdm2) 442 # XXX EO: old implementation: 443 # self.kernel_matrix = \ 444 # self.sigma_f * N.exp(-squared_euclidean_distance( 445 # data1, data2, weight=0.5 / (self.length_scale ** 2))) 446 return self.kernel_matrix

447

448 - def set_hyperparameters(self, hyperparameter):

449 """Set hyperaparmeters from a vector. 450 451 Used by model selection. 452 """ 453 if N.any(hyperparameter < 0): 454 raise InvalidHyperparameterError() 455 self.sigma_f = hyperparameter[0] 456 self._length_scale = hyperparameter[1:] 457 return

458

459 - def compute_lml_gradient(self,alphaalphaT_Kinv,data):

460 """Compute grandient of the kernel and return the portion of 461 log marginal likelihood gradient due to the kernel. 462 Shorter formula. Allows vector of lengthscales (ARD). 463 """ 464 self.lml_gradient = [] 465 def lml_grad(K_grad_i): 466 # return N.trace(N.dot(alphaalphaT_Kinv,K_grad_i)) 467 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 468 return (alphaalphaT_Kinv*(K_grad_i.T)).sum()

469 grad_sigma_f = 2.0/self.sigma_f*self.kernel_matrix 470 self.lml_gradient.append(lml_grad(grad_sigma_f)) 471 if N.isscalar(self.length_scale) or self.length_scale.size==1: 472 # use the same length_scale for all dimensions: 473 K_grad_l = self.wdm2*self.kernel_matrix*(1.0/self.length_scale) 474 self.lml_gradient.append(lml_grad(K_grad_l)) 475 else: 476 # use one length_scale for each dimension: 477 for i in range(self.length_scale.size): 478 K_grad_i = 1.0/(self.length_scale[i]**3)*self.kernel_matrix*N.subtract.outer(data[:,i],data[:,i])**2 479 self.lml_gradient.append(lml_grad(K_grad_i)) 480 pass 481 pass 482 self.lml_gradient = 0.5*N.array(self.lml_gradient) 483 return self.lml_gradient

484

485 - def compute_lml_gradient_logscale(self,alphaalphaT_Kinv,data):

486 """Compute grandient of the kernel and return the portion of 487 log marginal likelihood gradient due to the kernel. 488 Hyperparameters are in log scale which is sometimes more 489 stable. Shorter formula. Allows vector of lengthscales (ARD). 490 """ 491 self.lml_gradient = [] 492 def lml_grad(K_grad_i): 493 # return N.trace(N.dot(alphaalphaT_Kinv,K_grad_i)) 494 # Faster formula: N.trace(N.dot(A,B)) = (A*(B.T)).sum() 495 return (alphaalphaT_Kinv*(K_grad_i.T)).sum()

496 K_grad_log_sigma_f = 2.0*self.kernel_matrix 497 self.lml_gradient.append(lml_grad(K_grad_log_sigma_f)) 498 if N.isscalar(self.length_scale) or self.length_scale.size==1: 499 # use the same length_scale for all dimensions: 500 K_grad_log_l = self.wdm2*self.kernel_matrix 501 self.lml_gradient.append(lml_grad(K_grad_log_l)) 502 else: 503 # use one length_scale for each dimension: 504 for i in range(self.length_scale.size): 505 K_grad_log_l_i = 1.0/(self.length_scale[i]**2)*self.kernel_matrix*N.subtract.outer(data[:,i],data[:,i])**2 506 self.lml_gradient.append(lml_grad(K_grad_log_l_i)) 507 pass 508 pass 509 self.lml_gradient = 0.5*N.array(self.lml_gradient) 510 return self.lml_gradient 511

512 - def _setlength_scale(self, v):

513 """Set value of length_scale and its _orig 514 """ 515 self._length_scale = self._length_scale_orig = v

516 517 length_scale = property(fget=lambda x:x._length_scale, 518 fset=_setlength_scale) 519 pass 520

521 -class KernelMatern_3_2(Kernel):

522 """The Matern kernel class for the case ni=3/2 or ni=5/2. 523 524 Note that it can handle a length scale for each dimension for 525 Automtic Relevance Determination. 526 527 """

528 - def __init__(self, length_scale=1.0, sigma_f=1.0, numerator=3.0, **kwargs):

529 """Initialize a Squared Exponential kernel instance. 530 531 :Parameters: 532 length_scale : float OR numpy.ndarray 533 the characteristic length-scale (or length-scales) of the 534 phenomenon under investigation. 535 (Defaults to 1.0) 536 sigma_f : float 537 Signal standard deviation. 538 (Defaults to 1.0) 539 numerator: float 540 the numerator of parameter ni of Matern covariance functions. 541 Currently only numerator=3.0 and numerator=5.0 are implemented. 542 (Defaults to 3.0) 543 """ 544 # init base class first 545 Kernel.__init__(self, **kwargs) 546 547 self.length_scale = length_scale 548 self.sigma_f = sigma_f 549 self.kernel_matrix = None 550 if numerator == 3.0 or numerator == 5.0: 551 self.numerator = numerator 552 else: 553 raise NotImplementedError

554

555 - def __repr__(self):

556 return "%s(length_scale=%s, ni=%d/2)" \ 557 % (self.__class__.__name__, str(self.length_scale), self.numerator)

558

559 - def compute(self, data1, data2=None):

560 """Compute kernel matrix. 561 562 :Parameters: 563 data1 : numpy.ndarray 564 data 565 data2 : numpy.ndarray 566 data 567 (Defaults to None) 568 """ 569 tmp = squared_euclidean_distance( 570 data1, data2, weight=0.5 / (self.length_scale ** 2)) 571 if self.numerator == 3.0: 572 tmp = N.sqrt(tmp) 573 self.kernel_matrix = \ 574 self.sigma_f**2 * (1.0 + N.sqrt(3.0) * tmp) \ 575 * N.exp(-N.sqrt(3.0) * tmp) 576 elif self.numerator == 5.0: 577 tmp2 = N.sqrt(tmp) 578 self.kernel_matrix = \ 579 self.sigma_f**2 * (1.0 + N.sqrt(5.0) * tmp2 + 5.0 / 3.0 * tmp) \ 580 * N.exp(-N.sqrt(5.0) * tmp2) 581 return self.kernel_matrix

582

583 - def gradient(self, data1, data2):

584 """Compute gradient of the kernel matrix. A must for fast 585 model selection with high-dimensional data. 586 """ 587 # TODO SOON 588 # grad = ... 589 # return grad 590 raise NotImplementedError

591

592 - def set_hyperparameters(self, hyperparameter):

593 """Set hyperaparmeters from a vector. 594 595 Used by model selection. 596 Note: 'numerator' is not considered as an hyperparameter. 597 """ 598 if N.any(hyperparameter < 0): 599 raise InvalidHyperparameterError() 600 self.sigma_f = hyperparameter[0] 601 self.length_scale = hyperparameter[1:] 602 return

603 604 pass

605 606

607 -class KernelMatern_5_2(KernelMatern_3_2):

608 """The Matern kernel class for the case ni=5/2. 609 610 This kernel is just KernelMatern_3_2(numerator=5.0). 611 """

612 - def __init__(self, **kwargs):

613 """Initialize a Squared Exponential kernel instance. 614 615 :Parameters: 616 length_scale : float OR numpy.ndarray 617 the characteristic length-scale (or length-scales) of the 618 phenomenon under investigation. 619 (Defaults to 1.0) 620 """ 621 KernelMatern_3_2.__init__(self, numerator=5.0, **kwargs) 622 pass

623 624

625 -class KernelRationalQuadratic(Kernel):

626 """The Rational Quadratic (RQ) kernel class. 627 628 Note that it can handle a length scale for each dimension for 629 Automtic Relevance Determination. 630 631 """

632 - def __init__(self, length_scale=1.0, sigma_f=1.0, alpha=0.5, **kwargs):

633 """Initialize a Squared Exponential kernel instance. 634 635 :Parameters: 636 length_scale : float OR numpy.ndarray 637 the characteristic length-scale (or length-scales) of the 638 phenomenon under investigation. 639 (Defaults to 1.0) 640 sigma_f : float 641 Signal standard deviation. 642 (Defaults to 1.0) 643 alpha: float 644 The parameter of the RQ functions family. 645 (Defaults to 2.0) 646 """ 647 # init base class first 648 Kernel.__init__(self, **kwargs) 649 650 self.length_scale = length_scale 651 self.sigma_f = sigma_f 652 self.kernel_matrix = None 653 self.alpha = alpha

654

655 - def __repr__(self):

656 return "%s(length_scale=%s, alpha=%f)" \ 657 % (self.__class__.__name__, str(self.length_scale), self.alpha)

658

659 - def compute(self, data1, data2=None):

660 """Compute kernel matrix. 661 662 :Parameters: 663 data1 : numpy.ndarray 664 data 665 data2 : numpy.ndarray 666 data 667 (Defaults to None) 668 """ 669 tmp = squared_euclidean_distance( 670 data1, data2, weight=1.0 / (self.length_scale ** 2)) 671 self.kernel_matrix = \ 672 self.sigma_f**2 * (1.0 + tmp / (2.0 * self.alpha)) ** -self.alpha 673 return self.kernel_matrix

674

675 - def gradient(self, data1, data2):

676 """Compute gradient of the kernel matrix. A must for fast 677 model selection with high-dimensional data. 678 """ 679 # TODO SOON 680 # grad = ... 681 # return grad 682 raise NotImplementedError

683

684 - def set_hyperparameters(self, hyperparameter):

685 """Set hyperaparmeters from a vector. 686 687 Used by model selection. 688 Note: 'alpha' is not considered as an hyperparameter. 689 """ 690 if N.any(hyperparameter < 0): 691 raise InvalidHyperparameterError() 692 self.sigma_f = hyperparameter[0] 693 self.length_scale = hyperparameter[1:] 694 return

695 696 pass

697 698 699 # dictionary of avalable kernels with names as keys: 700 kernel_dictionary = {'constant': KernelConstant, 701 'linear': KernelLinear, 702 'exponential': KernelExponential, 703 'squared exponential': KernelSquaredExponential, 704 'Matern ni=3/2': KernelMatern_3_2, 705 'Matern ni=5/2': KernelMatern_5_2, 706 'rational quadratic': KernelRationalQuadratic} 707

Source Code for Module mvpa.clfs.kernel