apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-11 13:02:38 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add eq_ode1 covariance.

Browse source
This commit is contained in:
Neil Lawrence 2013-09-20 14:20:58 +01:00
parent f2750ff8ed
commit d9f895603b
3 changed files with 67 additions and 39 deletions
Download patch file Download diff file Expand all files Collapse all files

45
GPy/kern/constructors.py
View file

@ -140,6 +140,33 @@ def white(input_dim,variance=1.):
part = parts.white.White(input_dim,variance) part = parts.white.White(input_dim,variance)
return kern(input_dim, [part]) return kern(input_dim, [part])
def eq_ode1(output_dim, W=None, rank=1, kappa=None, length_scale=1., decay=None, delay=None):
"""Covariance function for first order differential equation driven by an exponentiated quadratic covariance.
This outputs of this kernel have the form
.. math::
\frac{\text{d}y_j}{\text{d}t} = \sum_{i=1}^R w_{j,i} f_i(t-\delta_j) +\sqrt{\kappa_j}g_j(t) - d_jy_j(t)
where :math:`R` is the rank of the system, :math:`w_{j,i}` is the sensitivity of the :math:`j`th output to the :math:`i`th latent function, :math:`d_j` is the decay rate of the :math:`j`th output and :math:`f_i(t)` and :math:`g_i(t)` are independent latent Gaussian processes goverened by an exponentiated quadratic covariance.
:param output_dim: number of outputs driven by latent function.
:type output_dim: int
:param W: sensitivities of each output to the latent driving function.
:type W: ndarray (output_dim x rank).
:param rank: If rank is greater than 1 then there are assumed to be a total of rank latent forces independently driving the system, each with identical covariance.
:type rank: int
:param decay: decay rates for the first order system.
:type decay: array of length output_dim.
:param delay: delay between latent force and output response.
:type delay: array of length output_dim.
:param kappa: diagonal term that allows each latent output to have an independent component to the response.
:type kappa: array of length output_dim.
.. Note: see first order differential equation examples in GPy.examples.regression for some usage.
"""
part = parts.eq_ode1.Eq_ode1(output_dim, W, rank, kappa, length_scale, decay, delay)
return kern(2, [part])
def exponential(input_dim,variance=1., lengthscale=None, ARD=False): def exponential(input_dim,variance=1., lengthscale=None, ARD=False):
""" """
@ -346,7 +373,7 @@ def symmetric(k):
k_.parts = [symmetric.Symmetric(p) for p in k.parts] k_.parts = [symmetric.Symmetric(p) for p in k.parts]
return k_ return k_
def coregionalize(num_outputs,W_columns=1, W=None, kappa=None): def coregionalize(output_dim,W_columns=1, W=None, kappa=None):
""" """
Coregionlization matrix B, of the form: Coregionlization matrix B, of the form:
.. math:: .. math::
@ -358,18 +385,18 @@ def coregionalize(num_outputs,W_columns=1, W=None, kappa=None):
it is obtainded as the tensor product between a kernel k(x,y) and B. it is obtainded as the tensor product between a kernel k(x,y) and B.
:param num_outputs: the number of outputs to corregionalize :param output_dim: the number of outputs to corregionalize
:type num_outputs: int :type output_dim: int
:param W_columns: number of columns of the W matrix (this parameter is ignored if parameter W is not None) :param W_columns: number of columns of the W matrix (this parameter is ignored if parameter W is not None)
:type W_colunns: int :type W_colunns: int
:param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalisation matrix B :param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalisation matrix B
:type W: numpy array of dimensionality (num_outpus, W_columns) :type W: numpy array of dimensionality (num_outpus, W_columns)
:param kappa: a vector which allows the outputs to behave independently :param kappa: a vector which allows the outputs to behave independently
:type kappa: numpy array of dimensionality (num_outputs,) :type kappa: numpy array of dimensionality (output_dim,)
:rtype: kernel object :rtype: kernel object
""" """
p = parts.coregionalize.Coregionalize(num_outputs,W_columns,W,kappa) p = parts.coregionalize.Coregionalize(output_dim,W_columns,W,kappa)
return kern(1,[p]) return kern(1,[p])
@ -429,12 +456,12 @@ def hierarchical(k):
_parts = [parts.hierarchical.Hierarchical(k.parts)] _parts = [parts.hierarchical.Hierarchical(k.parts)]
return kern(k.input_dim+len(k.parts),_parts) return kern(k.input_dim+len(k.parts),_parts)
def build_lcm(input_dim, num_outputs, kernel_list = [], W_columns=1,W=None,kappa=None): def build_lcm(input_dim, output_dim, kernel_list = [], W_columns=1,W=None,kappa=None):
""" """
Builds a kernel of a linear coregionalization model Builds a kernel of a linear coregionalization model
:input_dim: Input dimensionality :input_dim: Input dimensionality
:num_outputs: Number of outputs :output_dim: Number of outputs
:kernel_list: List of coregionalized kernels, each element in the list will be multiplied by a different corregionalization matrix :kernel_list: List of coregionalized kernels, each element in the list will be multiplied by a different corregionalization matrix
:type kernel_list: list of GPy kernels :type kernel_list: list of GPy kernels
:param W_columns: number tuples of the corregionalization parameters 'coregion_W' :param W_columns: number tuples of the corregionalization parameters 'coregion_W'
@ -448,11 +475,11 @@ def build_lcm(input_dim, num_outputs, kernel_list = [], W_columns=1,W=None,kappa
k.input_dim = input_dim k.input_dim = input_dim
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.") warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
k_coreg = coregionalize(num_outputs,W_columns,W,kappa) k_coreg = coregionalize(output_dim,W_columns,W,kappa)
kernel = kernel_list[0]**k_coreg.copy() kernel = kernel_list[0]**k_coreg.copy()
for k in kernel_list[1:]: for k in kernel_list[1:]:
k_coreg = coregionalize(num_outputs,W_columns,W,kappa) k_coreg = coregionalize(output_dim,W_columns,W,kappa)
kernel += k**k_coreg.copy() kernel += k**k_coreg.copy()
return kernel return kernel

3
GPy/kern/parts/__init__.py
View file

@ -2,10 +2,11 @@ import bias
import Brownian import Brownian
import coregionalize import coregionalize
import exponential import exponential
import eq_ode1
import finite_dimensional import finite_dimensional
import fixed import fixed
import gibbs import gibbs
#import hetero #hetero.py is not commited: omitting for now. JH. import hetero #hetero.py is not commited: omitting for now. JH.
import hierarchical import hierarchical
import independent_outputs import independent_outputs
import linear import linear

58
GPy/kern/parts/coregionalize.py
View file

@ -13,7 +13,7 @@ class Coregionalize(Kernpart):
This covariance has the form This covariance has the form
.. math:: .. math::
\mathbf{B} = \mathbf{W}\mathbf{W}^\top + kappa \mathbf{I} \mathbf{B} = \mathbf{W}\mathbf{W}^\top + \text{diag}(kappa)
An intrinsic/linear coregionalization covariance function of the form An intrinsic/linear coregionalization covariance function of the form
.. math:: .. math::
@ -22,33 +22,33 @@ class Coregionalize(Kernpart):
it is obtained as the tensor product between a covariance function it is obtained as the tensor product between a covariance function
k(x,y) and B. k(x,y) and B.
:param num_outputs: number of outputs to coregionalize :param output_dim: number of outputs to coregionalize
:type num_outputs: int :type output_dim: int
:param W_columns: number of columns of the W matrix (this parameter is ignored if parameter W is not None) :param W_columns: number of columns of the W matrix (this parameter is ignored if parameter W is not None)
:type W_colunns: int :type W_colunns: int
:param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalization matrix B :param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalization matrix B
:type W: numpy array of dimensionality (num_outpus, W_columns) :type W: numpy array of dimensionality (num_outpus, W_columns)
:param kappa: a vector which allows the outputs to behave independently :param kappa: a vector which allows the outputs to behave independently
:type kappa: numpy array of dimensionality (num_outputs,) :type kappa: numpy array of dimensionality (output_dim,)
.. Note: see coregionalization examples in GPy.examples.regression for some usage. .. Note: see coregionalization examples in GPy.examples.regression for some usage.
""" """
def __init__(self,num_outputs,W_columns=1, W=None, kappa=None): def __init__(self, output_dim, rank=1, W=None, kappa=None):
self.input_dim = 1 self.input_dim = 1
self.name = 'coregion' self.name = 'coregion'
self.num_outputs = num_outputs self.output_dim = output_dim
self.W_columns = W_columns self.rank = rank
if W is None: if W is None:
self.W = 0.5*np.random.randn(self.num_outputs,self.W_columns)/np.sqrt(self.W_columns) self.W = 0.5*np.random.randn(self.output_dim,self.rank)/np.sqrt(self.rank)
else: else:
assert W.shape==(self.num_outputs,self.W_columns) assert W.shape==(self.output_dim,self.rank)
self.W = W self.W = W
if kappa is None: if kappa is None:
kappa = 0.5*np.ones(self.num_outputs) kappa = 0.5*np.ones(self.output_dim)
else: else:
assert kappa.shape==(self.num_outputs,) assert kappa.shape==(self.output_dim,)
self.kappa = kappa self.kappa = kappa
self.num_params = self.num_outputs*(self.W_columns + 1) self.num_params = self.output_dim*(self.rank + 1)
self._set_params(np.hstack([self.W.flatten(),self.kappa])) self._set_params(np.hstack([self.W.flatten(),self.kappa]))
def _get_params(self): def _get_params(self):
@ -56,12 +56,12 @@ class Coregionalize(Kernpart):
def _set_params(self,x): def _set_params(self,x):
assert x.size == self.num_params assert x.size == self.num_params
self.kappa = x[-self.num_outputs:] self.kappa = x[-self.output_dim:]
self.W = x[:-self.num_outputs].reshape(self.num_outputs,self.W_columns) self.W = x[:-self.output_dim].reshape(self.output_dim,self.rank)
self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa) self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa)
def _get_param_names(self): def _get_param_names(self):
return sum([['W%i_%i'%(i,j) for j in range(self.W_columns)] for i in range(self.num_outputs)],[]) + ['kappa_%i'%i for i in range(self.num_outputs)] return sum([['W%i_%i'%(i,j) for j in range(self.rank)] for i in range(self.output_dim)],[]) + ['kappa_%i'%i for i in range(self.output_dim)]
def K(self,index,index2,target): def K(self,index,index2,target):
index = np.asarray(index,dtype=np.int) index = np.asarray(index,dtype=np.int)
@ -79,26 +79,26 @@ class Coregionalize(Kernpart):
if index2 is None: if index2 is None:
code=""" code="""
for(int i=0;i<N; i++){ for(int i=0;i<N; i++){
target[i+i*N] += B[index[i]+num_outputs*index[i]]; target[i+i*N] += B[index[i]+output_dim*index[i]];
for(int j=0; j<i; j++){ for(int j=0; j<i; j++){
target[j+i*N] += B[index[i]+num_outputs*index[j]]; target[j+i*N] += B[index[i]+output_dim*index[j]];
target[i+j*N] += target[j+i*N]; target[i+j*N] += target[j+i*N];
} }
} }
""" """
N,B,num_outputs = index.size, self.B, self.num_outputs N,B,output_dim = index.size, self.B, self.output_dim
weave.inline(code,['target','index','N','B','num_outputs']) weave.inline(code,['target','index','N','B','output_dim'])
else: else:
index2 = np.asarray(index2,dtype=np.int) index2 = np.asarray(index2,dtype=np.int)
code=""" code="""
for(int i=0;i<num_inducing; i++){ for(int i=0;i<num_inducing; i++){
for(int j=0; j<N; j++){ for(int j=0; j<N; j++){
target[i+j*num_inducing] += B[num_outputs*index[j]+index2[i]]; target[i+j*num_inducing] += B[output_dim*index[j]+index2[i]];
} }
} }
""" """
N,num_inducing,B,num_outputs = index.size,index2.size, self.B, self.num_outputs N,num_inducing,B,output_dim = index.size,index2.size, self.B, self.output_dim
weave.inline(code,['target','index','index2','N','num_inducing','B','num_outputs']) weave.inline(code,['target','index','index2','N','num_inducing','B','output_dim'])
def Kdiag(self,index,target): def Kdiag(self,index,target):
@ -115,12 +115,12 @@ class Coregionalize(Kernpart):
code=""" code="""
for(int i=0; i<num_inducing; i++){ for(int i=0; i<num_inducing; i++){
for(int j=0; j<N; j++){ for(int j=0; j<N; j++){
dL_dK_small[index[j] + num_outputs*index2[i]] += dL_dK[i+j*num_inducing]; dL_dK_small[index[j] + output_dim*index2[i]] += dL_dK[i+j*num_inducing];
} }
} }
""" """
N, num_inducing, num_outputs = index.size, index2.size, self.num_outputs N, num_inducing, output_dim = index.size, index2.size, self.output_dim
weave.inline(code, ['N','num_inducing','num_outputs','dL_dK','dL_dK_small','index','index2']) weave.inline(code, ['N','num_inducing','output_dim','dL_dK','dL_dK_small','index','index2'])
dkappa = np.diag(dL_dK_small) dkappa = np.diag(dL_dK_small)
dL_dK_small += dL_dK_small.T dL_dK_small += dL_dK_small.T
@ -137,8 +137,8 @@ class Coregionalize(Kernpart):
ii,jj = ii.T, jj.T ii,jj = ii.T, jj.T
dL_dK_small = np.zeros_like(self.B) dL_dK_small = np.zeros_like(self.B)
for i in range(self.num_outputs): for i in range(self.output_dim):
for j in range(self.num_outputs): for j in range(self.output_dim):
tmp = np.sum(dL_dK[(ii==i)*(jj==j)]) tmp = np.sum(dL_dK[(ii==i)*(jj==j)])
dL_dK_small[i,j] = tmp dL_dK_small[i,j] = tmp
@ -150,8 +150,8 @@ class Coregionalize(Kernpart):
def dKdiag_dtheta(self,dL_dKdiag,index,target): def dKdiag_dtheta(self,dL_dKdiag,index,target):
index = np.asarray(index,dtype=np.int).flatten() index = np.asarray(index,dtype=np.int).flatten()
dL_dKdiag_small = np.zeros(self.num_outputs) dL_dKdiag_small = np.zeros(self.output_dim)
for i in range(self.num_outputs): for i in range(self.output_dim):
dL_dKdiag_small[i] += np.sum(dL_dKdiag[index==i]) dL_dKdiag_small[i] += np.sum(dL_dKdiag[index==i])
dW = 2.*self.W*dL_dKdiag_small[:,None] dW = 2.*self.W*dL_dKdiag_small[:,None]
dkappa = dL_dKdiag_small dkappa = dL_dKdiag_small