Merge branch 'master' of github.com:SheffieldML/GPy

This commit is contained in:
James Hensman 2013-02-21 14:27:08 +00:00
commit a166187231
18 changed files with 347 additions and 86 deletions

View file

@ -333,8 +333,7 @@ class parameterised(object):
header_string = ["{h:^{col}}".format(h = header[i], col = cols[i]) for i in range(len(cols))] header_string = ["{h:^{col}}".format(h = header[i], col = cols[i]) for i in range(len(cols))]
header_string = map(lambda x: '|'.join(x), [header_string]) header_string = map(lambda x: '|'.join(x), [header_string])
separator = '-'*len(header_string[0]) separator = '-'*len(header_string[0])
param_string = ["{n:^{c0}}|{v:^{c1}}|{c:^{c2}}|{t:^{c3}}".format(n = names[i], v = values[i], c = constraints[i], t = ties[i], param_string = ["{n:^{c0}}|{v:^{c1}}|{c:^{c2}}|{t:^{c3}}".format(n = names[i], v = values[i], c = constraints[i], t = ties[i], c0 = cols[0], c1 = cols[1], c2 = cols[2], c3 = cols[3]) for i in range(len(values))]
c0 = cols[0], c1 = cols[1], c2 = cols[2], c3 = cols[3]) for i in range(len(values))]
return ('\n'.join([header_string[0], separator]+param_string)) + '\n' return ('\n'.join([header_string[0], separator]+param_string)) + '\n'

View file

@ -2,5 +2,5 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt) # Licensed under the BSD 3-clause license (see LICENSE.txt)
from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product_orthogonal from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product, product_orthogonal
from kern import kern from kern import kern

View file

@ -18,6 +18,7 @@ from Brownian import Brownian as Brownianpart
from periodic_exponential import periodic_exponential as periodic_exponentialpart from periodic_exponential import periodic_exponential as periodic_exponentialpart
from periodic_Matern32 import periodic_Matern32 as periodic_Matern32part from periodic_Matern32 import periodic_Matern32 as periodic_Matern32part
from periodic_Matern52 import periodic_Matern52 as periodic_Matern52part from periodic_Matern52 import periodic_Matern52 as periodic_Matern52part
from product import product as productpart
from product_orthogonal import product_orthogonal as product_orthogonalpart from product_orthogonal import product_orthogonal as product_orthogonalpart
#TODO these s=constructors are not as clean as we'd like. Tidy the code up #TODO these s=constructors are not as clean as we'd like. Tidy the code up
#using meta-classes to make the objects construct properly wthout them. #using meta-classes to make the objects construct properly wthout them.
@ -242,10 +243,21 @@ def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,
part = periodic_Matern52part(D,variance, lengthscale, period, n_freq, lower, upper) part = periodic_Matern52part(D,variance, lengthscale, period, n_freq, lower, upper)
return kern(D, [part]) return kern(D, [part])
def product(k1,k2):
"""
Construct a product kernel over D from two kernels over D
:param k1, k2: the kernels to multiply
:type k1, k2: kernpart
:rtype: kernel object
"""
part = productpart(k1,k2)
return kern(k1.D, [part])
def product_orthogonal(k1,k2): def product_orthogonal(k1,k2):
""" """
Construct a product kernel Construct a product kernel over D1 x D2 from a kernel over D1 and another over D2.
:param k1, k2: the kernels to multiply :param k1, k2: the kernels to multiply
:type k1, k2: kernpart :type k1, k2: kernpart
:rtype: kernel object :rtype: kernel object

View file

@ -8,6 +8,7 @@ from ..core.parameterised import parameterised
from kernpart import kernpart from kernpart import kernpart
import itertools import itertools
from product_orthogonal import product_orthogonal from product_orthogonal import product_orthogonal
from product import product
class kern(parameterised): class kern(parameterised):
def __init__(self,D,parts=[], input_slices=None): def __init__(self,D,parts=[], input_slices=None):
@ -149,24 +150,38 @@ class kern(parameterised):
""" """
Shortcut for `prod_orthogonal`. Note that `+` assumes that we sum 2 kernels defines on the same space whereas `*` assumes that the kernels are defined on different subspaces. Shortcut for `prod_orthogonal`. Note that `+` assumes that we sum 2 kernels defines on the same space whereas `*` assumes that the kernels are defined on different subspaces.
""" """
return self.prod_orthogonal(other) return self.prod(other)
def prod_orthogonal(self,other): def prod(self,other):
""" """
multiply two kernels. Both kernels are defined on separate spaces. Note that the constrains on the parameters of the kernels to multiply will be lost. multiply two kernels defined on the same spaces.
:param other: the other kernel to be added :param other: the other kernel to be added
:type other: GPy.kern :type other: GPy.kern
""" """
K1 = self.copy() K1 = self.copy()
K2 = other.copy() K2 = other.copy()
K1.unconstrain('')
K2.unconstrain('')
prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices] newkernparts = [product(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
K1.untie_everything()
K2.untie_everything()
D = K1.D + K2.D slices = []
for sl1, sl2 in itertools.product(K1.input_slices,K2.input_slices):
s1, s2 = [False]*K1.D, [False]*K2.D
s1[sl1], s2[sl2] = [True], [True]
slices += [s1+s2]
newkern = kern(K1.D, newkernparts, slices)
newkern._follow_constrains(K1,K2)
return newkern
def prod_orthogonal(self,other):
"""
multiply two kernels. Both kernels are defined on separate spaces.
:param other: the other kernel to be added
:type other: GPy.kern
"""
K1 = self.copy()
K2 = other.copy()
newkernparts = [product_orthogonal(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)] newkernparts = [product_orthogonal(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
@ -176,9 +191,14 @@ class kern(parameterised):
s1[sl1], s2[sl2] = [True], [True] s1[sl1], s2[sl2] = [True], [True]
slices += [s1+s2] slices += [s1+s2]
newkern = kern(D, newkernparts, slices) newkern = kern(K1.D + K2.D, newkernparts, slices)
newkern._follow_constrains(K1,K2)
# create the ties return newkern
def _follow_constrains(self,K1,K2):
# Build the array that allows to go from the initial indices of the param to the new ones
K1_param = [] K1_param = []
n = 0 n = 0
for k1 in K1.parts: for k1 in K1.parts:
@ -192,19 +212,40 @@ class kern(parameterised):
index_param = [] index_param = []
for p1 in K1_param: for p1 in K1_param:
for p2 in K2_param: for p2 in K2_param:
index_param += [0] + p1[1:] + p2[1:] index_param += p1 + p2
index_param = np.array(index_param) index_param = np.array(index_param)
# Get the ties and constrains of the kernels before the multiplication
prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices]
prev_constr_pos = np.append(K1.constrained_positive_indices, K1.Nparam + K2.constrained_positive_indices)
prev_constr_neg = np.append(K1.constrained_negative_indices, K1.Nparam + K2.constrained_negative_indices)
prev_constr_fix = K1.constrained_fixed_indices + [arr + K1.Nparam for arr in K2.constrained_fixed_indices]
prev_constr_fix_values = K1.constrained_fixed_values + K2.constrained_fixed_values
prev_constr_bou = K1.constrained_bounded_indices + [arr + K1.Nparam for arr in K2.constrained_bounded_indices]
prev_constr_bou_low = K1.constrained_bounded_lowers + K2.constrained_bounded_lowers
prev_constr_bou_upp = K1.constrained_bounded_uppers + K2.constrained_bounded_uppers
# follow the previous ties # follow the previous ties
for arr in prev_ties: for arr in prev_ties:
for j in arr: for j in arr:
index_param[np.where(index_param==j)[0]] = arr[0] index_param[np.where(index_param==j)[0]] = arr[0]
# tie # ties and constrains
for i in np.unique(index_param)[1:]: for i in range(K1.Nparam + K2.Nparam):
newkern.tie_param(np.where(index_param==i)[0]) index = np.where(index_param==i)[0]
if index.size > 1:
return newkern self.tie_param(index)
for i in prev_constr_pos:
self.constrain_positive(np.where(index_param==i)[0])
for i in prev_constr_neg:
self.constrain_neg(np.where(index_param==i)[0])
for j, i in enumerate(prev_constr_fix):
self.constrain_fixed(np.where(index_param==i)[0],prev_constr_fix_values[j])
for j, i in enumerate(prev_constr_bou):
self.constrain_bounded(np.where(index_param==i)[0],prev_constr_bou_low[j],prev_constr_bou_upp[j])
def _get_params(self): def _get_params(self):
return np.hstack([p._get_params() for p in self.parts]) return np.hstack([p._get_params() for p in self.parts])
@ -446,7 +487,7 @@ class kern(parameterised):
Xnew = np.vstack((xx.flatten(),yy.flatten())).T Xnew = np.vstack((xx.flatten(),yy.flatten())).T
Kx = self.K(Xnew,x,slices2=which_functions) Kx = self.K(Xnew,x,slices2=which_functions)
Kx = Kx.reshape(resolution,resolution).T Kx = Kx.reshape(resolution,resolution).T
pb.contour(xg,yg,Kx,vmin=Kx.min(),vmax=Kx.max(),cmap=pb.cm.jet) pb.contour(xg,yg,Kx,vmin=Kx.min(),vmax=Kx.max(),cmap=pb.cm.jet,*args,**kwargs)
pb.xlim(xmin[0],xmax[0]) pb.xlim(xmin[0],xmax[0])
pb.ylim(xmin[1],xmax[1]) pb.ylim(xmin[1],xmax[1])
pb.xlabel("x1") pb.xlabel("x1")

View file

@ -21,13 +21,14 @@ class periodic_Matern32(kernpart):
""" """
def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
assert D==1 assert D==1, "Periodic kernels are only defined for D=1"
self.name = 'periodic_Mat32' self.name = 'periodic_Mat32'
self.D = D self.D = D
if lengthscale is not None: if lengthscale is not None:
assert lengthscale.shape==(self.D,) lengthscale = np.asarray(lengthscale)
assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
else: else:
lengthscale = np.ones(self.D) lengthscale = np.ones(1)
self.lower,self.upper = lower, upper self.lower,self.upper = lower, upper
self.Nparam = 3 self.Nparam = 3
self.n_freq = n_freq self.n_freq = n_freq

View file

@ -21,13 +21,14 @@ class periodic_Matern52(kernpart):
""" """
def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
assert D==1 assert D==1, "Periodic kernels are only defined for D=1"
self.name = 'periodic_Mat52' self.name = 'periodic_Mat52'
self.D = D self.D = D
if lengthscale is not None: if lengthscale is not None:
assert lengthscale.shape==(self.D,) lengthscale = np.asarray(lengthscale)
assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
else: else:
lengthscale = np.ones(self.D) lengthscale = np.ones(1)
self.lower,self.upper = lower, upper self.lower,self.upper = lower, upper
self.Nparam = 3 self.Nparam = 3
self.n_freq = n_freq self.n_freq = n_freq

View file

@ -21,13 +21,14 @@ class periodic_exponential(kernpart):
""" """
def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
assert D==1 assert D==1, "Periodic kernels are only defined for D=1"
self.name = 'periodic_exp' self.name = 'periodic_exp'
self.D = D self.D = D
if lengthscale is not None: if lengthscale is not None:
assert lengthscale.shape==(self.D,) lengthscale = np.asarray(lengthscale)
assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
else: else:
lengthscale = np.ones(self.D) lengthscale = np.ones(1)
self.lower,self.upper = lower, upper self.lower,self.upper = lower, upper
self.Nparam = 3 self.Nparam = 3
self.n_freq = n_freq self.n_freq = n_freq
@ -45,7 +46,7 @@ class periodic_exponential(kernpart):
r = np.sqrt(r1**2 + r2**2) r = np.sqrt(r1**2 + r2**2)
psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2)) psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
return r,omega[:,0:1], psi return r,omega[:,0:1], psi
def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2): def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) ) Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + np.cos(phi1-phi2.T)*(self.upper-self.lower) Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + np.cos(phi1-phi2.T)*(self.upper-self.lower)

86
GPy/kern/product.py Normal file
View file

@ -0,0 +1,86 @@
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from kernpart import kernpart
import numpy as np
import hashlib
#from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
class product(kernpart):
"""
Computes the product of 2 kernels that are defined on the same space
:param k1, k2: the kernels to multiply
:type k1, k2: kernpart
:rtype: kernel object
"""
def __init__(self,k1,k2):
assert k1.D == k2.D, "Error: The input spaces of the kernels to multiply must have the same dimension"
self.D = k1.D
self.Nparam = k1.Nparam + k2.Nparam
self.name = k1.name + '<times>' + k2.name
self.k1 = k1
self.k2 = k2
self._set_params(np.hstack((k1._get_params(),k2._get_params())))
def _get_params(self):
"""return the value of the parameters."""
return self.params
def _set_params(self,x):
"""set the value of the parameters."""
self.k1._set_params(x[:self.k1.Nparam])
self.k2._set_params(x[self.k1.Nparam:])
self.params = x
def _get_param_names(self):
"""return parameter names."""
return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
def K(self,X,X2,target):
"""Compute the covariance matrix between X and X2."""
if X2 is None: X2 = X
target1 = np.zeros((X.shape[0],X2.shape[0]))
target2 = np.zeros((X.shape[0],X2.shape[0]))
self.k1.K(X,X2,target1)
self.k2.K(X,X2,target2)
target += target1 * target2
def Kdiag(self,X,target):
"""Compute the diagonal of the covariance matrix associated to X."""
target1 = np.zeros((X.shape[0],))
target2 = np.zeros((X.shape[0],))
self.k1.Kdiag(X,target1)
self.k2.Kdiag(X,target2)
target += target1 * target2
def dK_dtheta(self,partial,X,X2,target):
"""derivative of the covariance matrix with respect to the parameters."""
if X2 is None: X2 = X
K1 = np.zeros((X.shape[0],X2.shape[0]))
K2 = np.zeros((X.shape[0],X2.shape[0]))
self.k1.K(X,X2,K1)
self.k2.K(X,X2,K2)
k1_target = np.zeros(self.k1.Nparam)
k2_target = np.zeros(self.k2.Nparam)
self.k1.dK_dtheta(partial*K2, X, X2, k1_target)
self.k2.dK_dtheta(partial*K1, X, X2, k2_target)
target[:self.k1.Nparam] += k1_target
target[self.k1.Nparam:] += k2_target
def dK_dX(self,partial,X,X2,target):
"""derivative of the covariance matrix with respect to X."""
if X2 is None: X2 = X
K1 = np.zeros((X.shape[0],X2.shape[0]))
K2 = np.zeros((X.shape[0],X2.shape[0]))
self.k1.K(X,X2,K1)
self.k2.K(X,X2,K2)
self.k1.dK_dX(partial*K2, X, X2, target)
self.k2.dK_dX(partial*K1, X, X2, target)
def dKdiag_dX(self,X,target):
pass

View file

@ -4,7 +4,7 @@
from kernpart import kernpart from kernpart import kernpart
import numpy as np import numpy as np
import hashlib import hashlib
from scipy import integrate #from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
class product_orthogonal(kernpart): class product_orthogonal(kernpart):
""" """
@ -16,13 +16,12 @@ class product_orthogonal(kernpart):
""" """
def __init__(self,k1,k2): def __init__(self,k1,k2):
assert k1._get_param_names()[0] == 'variance' and k2._get_param_names()[0] == 'variance', "Error: The multipication of kernels is only defined when the first parameters of the kernels to multiply is the variance."
self.D = k1.D + k2.D self.D = k1.D + k2.D
self.Nparam = k1.Nparam + k2.Nparam - 1 self.Nparam = k1.Nparam + k2.Nparam
self.name = k1.name + '<times>' + k2.name self.name = k1.name + '<times>' + k2.name
self.k1 = k1 self.k1 = k1
self.k2 = k2 self.k2 = k2
self._set_params(np.hstack((k1._get_params()[0]*k2._get_params()[0], k1._get_params()[1:],k2._get_params()[1:]))) self._set_params(np.hstack((k1._get_params(),k2._get_params())))
def _get_params(self): def _get_params(self):
"""return the value of the parameters.""" """return the value of the parameters."""
@ -30,14 +29,14 @@ class product_orthogonal(kernpart):
def _set_params(self,x): def _set_params(self,x):
"""set the value of the parameters.""" """set the value of the parameters."""
self.k1._set_params(np.hstack((1.,x[1:self.k1.Nparam]))) self.k1._set_params(x[:self.k1.Nparam])
self.k2._set_params(np.hstack((1.,x[self.k1.Nparam:]))) self.k2._set_params(x[self.k1.Nparam:])
self.params = x self.params = x
def _get_param_names(self): def _get_param_names(self):
"""return parameter names.""" """return parameter names."""
return ['variance']+[self.k1.name + '_' + self.k1._get_param_names()[i+1] for i in range(self.k1.Nparam-1)] + [self.k2.name + '_' + self.k2._get_param_names()[i+1] for i in range(self.k2.Nparam-1)] return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
def K(self,X,X2,target): def K(self,X,X2,target):
"""Compute the covariance matrix between X and X2.""" """Compute the covariance matrix between X and X2."""
if X2 is None: X2 = X if X2 is None: X2 = X
@ -45,7 +44,7 @@ class product_orthogonal(kernpart):
target2 = np.zeros((X.shape[0],X2.shape[0])) target2 = np.zeros((X.shape[0],X2.shape[0]))
self.k1.K(X[:,0:self.k1.D],X2[:,0:self.k1.D],target1) self.k1.K(X[:,0:self.k1.D],X2[:,0:self.k1.D],target1)
self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2) self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2)
target += self.params[0]*target1 * target2 target += target1 * target2
def Kdiag(self,X,target): def Kdiag(self,X,target):
"""Compute the diagonal of the covariance matrix associated to X.""" """Compute the diagonal of the covariance matrix associated to X."""
@ -53,7 +52,7 @@ class product_orthogonal(kernpart):
target2 = np.zeros((X.shape[0],)) target2 = np.zeros((X.shape[0],))
self.k1.Kdiag(X[:,0:self.k1.D],target1) self.k1.Kdiag(X[:,0:self.k1.D],target1)
self.k2.Kdiag(X[:,self.k1.D:],target2) self.k2.Kdiag(X[:,self.k1.D:],target2)
target += self.params[0]*target1 * target2 target += target1 * target2
def dK_dtheta(self,partial,X,X2,target): def dK_dtheta(self,partial,X,X2,target):
"""derivative of the covariance matrix with respect to the parameters.""" """derivative of the covariance matrix with respect to the parameters."""
@ -68,13 +67,8 @@ class product_orthogonal(kernpart):
self.k1.dK_dtheta(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], k1_target) self.k1.dK_dtheta(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], k1_target)
self.k2.dK_dtheta(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], k2_target) self.k2.dK_dtheta(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], k2_target)
target[0] += np.sum(K1*K2*partial) target[:self.k1.Nparam] += k1_target
target[1:self.k1.Nparam] += self.params[0]* k1_target[1:] target[self.k1.Nparam:] += k2_target
target[self.k1.Nparam:] += self.params[0]* k2_target[1:]
def dKdiag_dtheta(self,partial,X,target):
"""derivative of the diagonal of the covariance matrix with respect to the parameters."""
target[0] += 1
def dK_dX(self,partial,X,X2,target): def dK_dX(self,partial,X,X2,target):
"""derivative of the covariance matrix with respect to X.""" """derivative of the covariance matrix with respect to X."""

Binary file not shown.

Before

Width:  |  Height:  |  Size: 38 KiB

View file

Before

Width:  |  Height:  |  Size: 38 KiB

After

Width:  |  Height:  |  Size: 38 KiB

Before After
Before After

Binary file not shown.

Before

Width:  |  Height:  |  Size: 55 KiB

After

Width:  |  Height:  |  Size: 37 KiB

Before After
Before After

Binary file not shown.

Before

Width:  |  Height:  |  Size: 84 KiB

After

Width:  |  Height:  |  Size: 78 KiB

Before After
Before After

Binary file not shown.

After

Width:  |  Height:  |  Size: 61 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 40 KiB

View file

@ -1,6 +1,14 @@
inference Package inference Package
================= =================
:mod:`SGD` Module
-----------------
.. automodule:: GPy.inference.SGD
:members:
:undoc-members:
:show-inheritance:
:mod:`optimization` Module :mod:`optimization` Module
-------------------------- --------------------------

View file

@ -113,6 +113,14 @@ kern Package
:undoc-members: :undoc-members:
:show-inheritance: :show-inheritance:
:mod:`product` Module
---------------------
.. automodule:: GPy.kern.product
:members:
:undoc-members:
:show-inheritance:
:mod:`product_orthogonal` Module :mod:`product_orthogonal` Module
-------------------------------- --------------------------------

View file

@ -13,20 +13,27 @@ First we import the libraries we will need ::
For most kernels, the dimension is the only mandatory parameter to define a kernel object. However, it is also possible to specify the values of the parameters. For example, the three following commands are valid for defining a squared exponential kernel (ie rbf or Gaussian) :: For most kernels, the dimension is the only mandatory parameter to define a kernel object. However, it is also possible to specify the values of the parameters. For example, the three following commands are valid for defining a squared exponential kernel (ie rbf or Gaussian) ::
ker1 = GPy.kern.rbf(1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.) ker1 = GPy.kern.rbf(1) # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.)
ker2 = GPy.kern.rbf(D=1, variance = 1.5, lengthscale=2.) ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=2.)
ker3 = GPy.kern.rbf(1, .5, .5) ker3 = GPy.kern.rbf(1, .5, .5)
A `plot` and a `print` functions are implemented to represent kernel objects :: A ``print`` and a ``plot`` functions are implemented to represent kernel objects. The commands ::
print ker1 print ker2
ker1.plot() ker1.plot()
ker2.plot() ker2.plot()
ker3.plot() ker3.plot()
.. figure:: Figures/tuto_kern_overview_basicdef.png return::
Name | Value | Constraints | Ties
-------------------------------------------------------
rbf_variance | 0.7500 | |
rbf_lengthscale | 2.0000 | |
.. figure:: Figures/tuto_kern_overview_basicplot.png
:align: center :align: center
:height: 350px :height: 300px
Implemented kernels Implemented kernels
=================== ===================
@ -37,33 +44,131 @@ Many kernels are already implemented in GPy. Here is a summary of most of them:
:align: center :align: center
:height: 800px :height: 800px
On the other hand, it is possible to use the `sympy` package to build new kernels. This will be the subject of another tutorial. On the other hand, it is possible to use the `sympy` package to build new kernels. This will be the subject of another tutorial.
Operations to combine kernel Operations to combine kernels
============================ =============================
In ``GPy``, kernel objects can be combined with the usual ``+`` and ``*`` operators. :: In ``GPy``, kernel objects can be added or multiplied. In both cases, two kinds of operations are possible since one can assume that the kernels to add/multiply are defined on the same space or on different subspaces. In other words, it is possible to use two kernels :math:`k_1,\ k_2` over :math:`\mathbb{R} \times \mathbb{R}` to create
k1 = GPy.kern.rbf(1,variance=1., lengthscale=2)
k2 = GPy.kern.Matern32(1,variance=1., lengthscale=2)
ker_add = k1 + k2 * a kernel over :math:`\mathbb{R} \times \mathbb{R}`: :math:`k(x,y) = k_1(x,y) \times k_2(x,y)`
print ker_add * a kernel over :math:`\mathbb{R}^2 \times \mathbb{R}^2`: :math:`k(\mathbf{x},\mathbf{y}) = k_1(x_1,y_1) \times k_2(x_2,y_2)`
ker_prod = k1 * k2 These two options are available in GPy under the name ``prod`` and ``prod_orthogonal`` (resp. ``add`` and ``add_orthogonal`` for the addition). Here is a quick example ::
print ker_prod
Note that by default, the operator ``+`` adds kernels defined on the same input space whereas ``*`` assumes that the kernels are defined on different input spaces. Here for example ``ker_add.D`` will return ``1`` whereas ``ker_prod.D`` will return ``2``. k1 = GPy.kern.rbf(1,1.,2.)
k2 = GPy.kern.Matern32(1, 0.5, 0.2)
In order to add kernels defined on the different input spaces, the required command is:: # Product of kernels
k_prod = k1.prod(k2)
k_prodorth = k1.prod_orthogonal(k2)
ker_add_orth = k1.add_orthogonal(k2) # Sum of kernels
k_add = k1.add(k2)
k_addorth = k1.add_orthogonal(k2)
.. figure:: Figures/tuto_kern_overview_add_orth.png .. # plots
pb.figure(figsize=(8,8))
pb.subplot(2,2,1)
k_prod.plot()
pb.title('prod')
pb.subplot(2,2,2)
k_prodorth.plot()
pb.title('prod_orthogonal')
pb.subplot(2,2,3)
k_add.plot()
pb.title('add')
pb.subplot(2,2,4)
k_addorth.plot()
pb.title('add_orthogonal')
pb.subplots_adjust(wspace=0.3, hspace=0.3)
.. figure:: Figures/tuto_kern_overview_multadd.png
:align: center :align: center
:height: 350px :height: 500px
Output of ``ker_add_orth.plot(plot_limits=[[-10,-10],[10,10]])``. A shortcut for ``add`` and ``prod`` is provided by the usual ``+`` and ``*`` operators. Here is another example where we create a periodic kernel with some decay ::
k1 = GPy.kern.rbf(1,1.,2)
k2 = GPy.kern.periodic_Matern52(1,variance=1e3, lengthscale=1, period = 1.5, lower=-5., upper = 5)
k = k1 * k2 # equivalent to k = k1.prod(k2)
print k
# Simulate sample paths
X = np.linspace(-5,5,501)[:,None]
Y = np.random.multivariate_normal(np.zeros(501),k.K(X),1)
.. # plot
pb.figure(figsize=(10,4))
pb.subplot(1,2,1)
k.plot()
pb.subplot(1,2,2)
pb.plot(X,Y.T)
pb.ylabel("Sample path")
pb.subplots_adjust(wspace=0.3)
.. figure:: Figures/tuto_kern_overview_multperdecay.png
:align: center
:height: 300px
In general, ``kern`` objects can be seen as a sum of ``kernparts`` objects, where the later are covariance functions denied on the same space. For example, the following code ::
k = (k1+k2)*(k1+k2)
print k.parts[0].name, '\n', k.parts[1].name, '\n', k.parts[2].name, '\n', k.parts[3].name
returns ::
rbf<times>rbf
rbf<times>periodic_Mat52
periodic_Mat52<times>rbf
periodic_Mat52<times>periodic_Mat52
Constraining the parameters
===========================
Various constrains can be applied to the parameters of a kernel
* ``constrain_fixed`` to fix the value of a parameter (the value will not be modified during optimisation)
* ``constrain_positive`` to make sure the parameter is greater than 0.
* ``constrain_bounded`` to impose the parameter to be in a given range.
* ``tie_param`` to impose the value of two (or more) parameters to be equal.
When calling one of these functions, the parameters to constrain can either by specified by a regular expression that matches its name or by a number that corresponds to the rank of the parameter. Here is an example ::
k1 = GPy.kern.rbf(1)
k2 = GPy.kern.Matern32(1)
k3 = GPy.kern.white(1)
k = k1 + k2 + k3
print k
k.constrain_positive('var')
k.constrain_fixed(np.array([1]),1.75)
k.tie_param('len')
k.unconstrain('white')
k.constrain_bounded('white',lower=1e-5,upper=.5)
print k
with output::
Name | Value | Constraints | Ties
---------------------------------------------------------
rbf_variance | 1.0000 | |
rbf_lengthscale | 1.0000 | |
Mat32_variance | 1.0000 | |
Mat32_lengthscale | 1.0000 | |
white_variance | 1.0000 | |
Name | Value | Constraints | Ties
----------------------------------------------------------
rbf_variance | 1.0000 | (+ve) |
rbf_lengthscale | 1.7500 | Fixed | (0)
Mat32_variance | 1.0000 | (+ve) |
Mat32_lengthscale | 1.7500 | | (0)
white_variance | 0.3655 | (1e-05, 0.5) |
Example : Building an ANOVA kernel Example : Building an ANOVA kernel
================================== ==================================
@ -78,23 +183,27 @@ Let us assume that we want to define an ANOVA kernel with a Matern 3/2 kernel fo
k_cst = GPy.kern.bias(1,variance=1.) k_cst = GPy.kern.bias(1,variance=1.)
k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3) k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3)
Kanova = (k_cst + k_mat) * (k_cst + k_mat) Kanova = (k_cst + k_mat).prod_orthogonal(k_cst + k_mat)
print Kanova print Kanova
Printing the resulting kernel outputs the following :: Printing the resulting kernel outputs the following ::
Name | Value | Constraints | Ties Name | Value | Constraints | Ties
--------------------------------------------------------------------------- ---------------------------------------------------------------------------
bias<times>bias_variance | 1.0000 | | bias<times>bias_bias_variance | 1.0000 | | (0)
bias<times>Mat52_variance | 1.0000 | | bias<times>bias_bias_variance | 1.0000 | | (3)
bias<times>Mat52_Mat52_lengthscale | 3.0000 | | (1) bias<times>Mat52_bias_variance | 1.0000 | | (0)
Mat52<times>bias_variance | 1.0000 | | bias<times>Mat52_Mat52_variance | 1.0000 | | (4)
Mat52<times>bias_Mat52_lengthscale | 3.0000 | | (0) bias<times>Mat52_Mat52_lengthscale | 3.0000 | | (5)
Mat52<times>Mat52_variance | 1.0000 | | Mat52<times>bias_Mat52_variance | 1.0000 | | (1)
Mat52<times>Mat52_Mat52_lengthscale | 3.0000 | | (0) Mat52<times>bias_Mat52_lengthscale | 3.0000 | | (2)
Mat52<times>Mat52_Mat52_lengthscale | 3.0000 | | (1) Mat52<times>bias_bias_variance | 1.0000 | | (3)
Mat52<times>Mat52_Mat52_variance | 1.0000 | | (1)
Mat52<times>Mat52_Mat52_lengthscale | 3.0000 | | (2)
Mat52<times>Mat52_Mat52_variance | 1.0000 | | (4)
Mat52<times>Mat52_Mat52_lengthscale | 3.0000 | | (5)
Note the ties between the lengthscales of ``Kanova`` to keep the number of lengthscales equal to 2. On the other hand, there are four variance terms in the new parameterization: one for each term of the right hand part of the above equation. We can illustrate the use of this kernel on a toy example:: Note the ties between the parameters of ``Kanova`` that reflect the links between the parameters of the kernparts objects. We can illustrate the use of this kernel on a toy example::
# sample inputs and outputs # sample inputs and outputs
X = np.random.uniform(-3.,3.,(40,2)) X = np.random.uniform(-3.,3.,(40,2))
@ -102,12 +211,12 @@ Note the ties between the lengthscales of ``Kanova`` to keep the number of lengt
# Create GP regression model # Create GP regression model
m = GPy.models.GP_regression(X,Y,Kanova) m = GPy.models.GP_regression(X,Y,Kanova)
pb.figure(figsize=(5,5))
m.plot() m.plot()
.. figure:: Figures/tuto_kern_overview_mANOVA.png .. figure:: Figures/tuto_kern_overview_mANOVA.png
:align: center :align: center
:height: 350px :height: 300px
As :math:`k_{ANOVA}` corresponds to the sum of 4 kernels, the best predictor can be splited in a sum of 4 functions As :math:`k_{ANOVA}` corresponds to the sum of 4 kernels, the best predictor can be splited in a sum of 4 functions
@ -119,7 +228,7 @@ As :math:`k_{ANOVA}` corresponds to the sum of 4 kernels, the best predictor can
The submodels can be represented with the option ``which_function`` of ``plot``: :: The submodels can be represented with the option ``which_function`` of ``plot``: ::
pb.figure(figsize=(20,5)) pb.figure(figsize=(20,3))
pb.subplots_adjust(wspace=0.5) pb.subplots_adjust(wspace=0.5)
pb.subplot(1,5,1) pb.subplot(1,5,1)
m.plot() m.plot()
@ -135,10 +244,11 @@ The submodels can be represented with the option ``which_function`` of ``plot``:
pb.ylabel("+ ",rotation='horizontal',fontsize='30') pb.ylabel("+ ",rotation='horizontal',fontsize='30')
m.plot(which_functions=[False,False,False,True]) m.plot(which_functions=[False,False,False,True])
.. pb.savefig('tuto_kern_overview_mANOVAdec.png',bbox_inches='tight')
.. figure:: Figures/tuto_kern_overview_mANOVAdec.png .. figure:: Figures/tuto_kern_overview_mANOVAdec.png
:align: center :align: center
:height: 200px :height: 250px
.. import pylab as pb .. import pylab as pb