Various doc string edits and consistency changes in comments.

This commit is contained in:
Neil Lawrence 2013-08-19 20:11:25 +02:00
parent 7bb7c0a08c
commit c140b2118a
5 changed files with 218 additions and 79 deletions

View file

@ -6,8 +6,8 @@ from GPy.core.model import Model
class GPBase(Model): class GPBase(Model):
""" """
Gaussian Process Model for holding shared behaviour between Gaussian process base model for holding shared behaviour between
sprase_GP and GP models sparse_GP and GP models.
""" """
def __init__(self, X, likelihood, kernel, normalize_X=False): def __init__(self, X, likelihood, kernel, normalize_X=False):
@ -35,8 +35,7 @@ class GPBase(Model):
def getstate(self): def getstate(self):
""" """
Get the current state of the class, Get the current state of the class, here we return everything that is needed to recompute the model.
here just all the indices, rest can get recomputed
""" """
return Model.getstate(self) + [self.X, return Model.getstate(self) + [self.X,
self.num_data, self.num_data,
@ -63,14 +62,6 @@ class GPBase(Model):
Plot the GP's view of the world, where the data is normalized and the Plot the GP's view of the world, where the data is normalized and the
likelihood is Gaussian. likelihood is Gaussian.
:param samples: the number of a posteriori samples to plot
:param which_data: which if the training data to plot (default all)
:type which_data: 'all' or a slice object to slice self.X, self.Y
:param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
:param which_parts: which of the kernel functions to plot (additively)
:type which_parts: 'all', or list of bools
:param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
Plot the posterior of the GP. Plot the posterior of the GP.
- In one dimension, the function is plotted with a shaded region identifying two standard deviations. - In one dimension, the function is plotted with a shaded region identifying two standard deviations.
- In two dimsensions, a contour-plot shows the mean predicted function - In two dimsensions, a contour-plot shows the mean predicted function
@ -78,6 +69,22 @@ class GPBase(Model):
Can plot only part of the data and part of the posterior functions Can plot only part of the data and part of the posterior functions
using which_data and which_functions using which_data and which_functions
:param samples: the number of a posteriori samples to plot
:param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
:param which_data: which if the training data to plot (default all)
:type which_data: 'all' or a slice object to slice self.X, self.Y
:param which_parts: which of the kernel functions to plot (additively)
:type which_parts: 'all', or list of bools
:param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
:type resolution: int
:param full_cov:
:type full_cov: bool
:param fignum: figure to plot on.
:type fignum: figure number
:param ax: axes to plot on.
:type ax: axes handle
""" """
if which_data == 'all': if which_data == 'all':
which_data = slice(None) which_data = slice(None)
@ -118,12 +125,41 @@ class GPBase(Model):
def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']): def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']):
""" """
TODO: Docstrings! Plot the GP with noise where the likelihood is Gaussian.
Plot the posterior of the GP.
- In one dimension, the function is plotted with a shaded region identifying two standard deviations.
- In two dimsensions, a contour-plot shows the mean predicted function
- In higher dimensions, we've no implemented this yet !TODO!
Can plot only part of the data and part of the posterior functions
using which_data and which_functions
:param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
:type plot_limits: np.array
:param which_data: which if the training data to plot (default all)
:type which_data: 'all' or a slice object to slice self.X, self.Y
:param which_parts: which of the kernel functions to plot (additively)
:type which_parts: 'all', or list of bools
:param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
:type resolution: int
:param levels: number of levels to plot in a contour plot.
:type levels: int
:param samples: the number of a posteriori samples to plot
:type samples: int
:param fignum: figure to plot on.
:type fignum: figure number
:param ax: axes to plot on.
:type ax: axes handle
:param fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input index i should be set to value v.
:type fixed_inputs: a list of tuples
:param linecol: color of line to plot.
:type linecol:
:param fillcol: color of fill
:type fillcol:
:param levels: for 2D plotting, the number of contour levels to use :param levels: for 2D plotting, the number of contour levels to use
is ax is None, create a new figure is ax is None, create a new figure
fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input index i should be set to value v.
""" """
# TODO include samples # TODO include samples
if which_data == 'all': if which_data == 'all':

View file

@ -24,15 +24,17 @@ class Model(Parameterized):
self.preferred_optimizer = 'scg' self.preferred_optimizer = 'scg'
# self._set_params(self._get_params()) has been taken out as it should only be called on leaf nodes # self._set_params(self._get_params()) has been taken out as it should only be called on leaf nodes
def log_likelihood(self): def log_likelihood(self):
raise NotImplementedError, "this needs to be implemented to use the Model class" raise NotImplementedError, "this needs to be implemented to use the model class"
def _log_likelihood_gradients(self): def _log_likelihood_gradients(self):
raise NotImplementedError, "this needs to be implemented to use the Model class" raise NotImplementedError, "this needs to be implemented to use the model class"
def getstate(self): def getstate(self):
""" """
Get the current state of the class. Get the current state of the class.
Inherited from Parameterized, so add those parameters to the state Inherited from Parameterized, so add those parameters to the state
:return: list of states from the model.
""" """
return Parameterized.getstate(self) + \ return Parameterized.getstate(self) + \
[self.priors, self.optimization_runs, [self.priors, self.optimization_runs,
@ -41,8 +43,10 @@ class Model(Parameterized):
def setstate(self, state): def setstate(self, state):
""" """
set state from previous call to getstate set state from previous call to getstate
call Parameterized with the rest of the state call Parameterized with the rest of the state
:param state: the state of the model.
:type state: list as returned from getstate.
""" """
self.preferred_optimizer = state.pop() self.preferred_optimizer = state.pop()
self.sampling_runs = state.pop() self.sampling_runs = state.pop()
@ -52,21 +56,22 @@ class Model(Parameterized):
def set_prior(self, regexp, what): def set_prior(self, regexp, what):
""" """
Sets priors on the Model parameters. Sets priors on the model parameters.
Arguments
---------
regexp -- string, regexp, or integer array
what -- instance of a Prior class
Notes Notes
----- -----
Asserts that the Prior is suitable for the constraint. If the Asserts that the prior is suitable for the constraint. If the
wrong constraint is in place, an error is raised. If no wrong constraint is in place, an error is raised. If no
constraint is in place, one is added (warning printed). constraint is in place, one is added (warning printed).
For tied parameters, the Prior will only be "counted" once, thus For tied parameters, the prior will only be "counted" once, thus
a Prior object is only inserted on the first tied index a prior object is only inserted on the first tied index
:param regexp: regular expression of parameters on which priors need to be set.
:type param: string, regexp, or integer array
:param what: prior to set on parameter.
:type what: GPy.core.Prior type
""" """
if self.priors is None: if self.priors is None:
self.priors = [None for i in range(self._get_params().size)] self.priors = [None for i in range(self._get_params().size)]
@ -76,12 +81,12 @@ class Model(Parameterized):
# check tied situation # check tied situation
tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))] tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))]
if len(tie_partial_matches): if len(tie_partial_matches):
raise ValueError, "cannot place Prior across partial ties" raise ValueError, "cannot place prior across partial ties"
tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ] tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ]
if len(tie_matches) > 1: if len(tie_matches) > 1:
raise ValueError, "cannot place Prior across multiple ties" raise ValueError, "cannot place prior across multiple ties"
elif len(tie_matches) == 1: elif len(tie_matches) == 1:
which = which[:1] # just place a Prior object on the first parameter which = which[:1] # just place a prior object on the first parameter
# check constraints are okay # check constraints are okay
@ -93,7 +98,7 @@ class Model(Parameterized):
else: else:
constrained_positive_indices = np.zeros(shape=(0,)) constrained_positive_indices = np.zeros(shape=(0,))
bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices) bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices)
assert not np.any(which[:, None] == bad_constraints), "constraint and Prior incompatible" assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible"
unconst = np.setdiff1d(which, constrained_positive_indices) unconst = np.setdiff1d(which, constrained_positive_indices)
if len(unconst): if len(unconst):
print "Warning: constraining parameters to be positive:" print "Warning: constraining parameters to be positive:"
@ -101,17 +106,22 @@ class Model(Parameterized):
print '\n' print '\n'
self.constrain_positive(unconst) self.constrain_positive(unconst)
elif what.domain is REAL: elif what.domain is REAL:
assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and Prior incompatible" assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and prior incompatible"
else: else:
raise ValueError, "Prior not recognised" raise ValueError, "prior not recognised"
# store the Prior in a local list # store the prior in a local list
for w in which: for w in which:
self.priors[w] = what self.priors[w] = what
def get_gradient(self, name, return_names=False): def get_gradient(self, name, return_names=False):
""" """
Get Model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. Get model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned.
:param name: the name of parameters required (as a regular expression).
:type name: regular expression
:param return_names: whether or not to return the names matched (default False)
:type return_names: bool
""" """
matches = self.grep_param_names(name) matches = self.grep_param_names(name)
if len(matches): if len(matches):
@ -151,14 +161,14 @@ class Model(Parameterized):
def randomize(self): def randomize(self):
""" """
Randomize the Model. Randomize the model.
Make this draw from the Prior if one exists, else draw from N(0,1) Make this draw from the prior if one exists, else draw from N(0,1)
""" """
# first take care of all parameters (from N(0,1)) # first take care of all parameters (from N(0,1))
x = self._get_params_transformed() x = self._get_params_transformed()
x = np.random.randn(x.size) x = np.random.randn(x.size)
self._set_params_transformed(x) self._set_params_transformed(x)
# now draw from Prior where possible # now draw from prior where possible
x = self._get_params() x = self._get_params()
if self.priors is not None: if self.priors is not None:
[np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) if not p is None] [np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) if not p is None]
@ -168,21 +178,30 @@ class Model(Parameterized):
def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs): def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs):
""" """
Perform random restarts of the Model, and set the Model to the best Perform random restarts of the model, and set the model to the best
seen solution. seen solution.
If the robust flag is set, exceptions raised during optimizations will If the robust flag is set, exceptions raised during optimizations will
be handled silently. If _all_ runs fail, the Model is reset to the be handled silently. If _all_ runs fail, the model is reset to the
existing parameter values. existing parameter values.
Notes Notes
----- -----
:param num_restarts: number of restarts to use (default 10)
:type num_restarts: int
:param robust: whether to handle exceptions silently or not (default False)
:type robust: bool
:param parallel: whether to run each restart as a separate process. It relies on the multiprocessing module.
:type parallel: bool
:param num_processes: number of workers in the multiprocessing pool
:type numprocesses: int
**kwargs are passed to the optimizer. They can be: **kwargs are passed to the optimizer. They can be:
:max_f_eval: maximum number of function evaluations :param max_f_eval: maximum number of function evaluations
:messages: whether to display during optimisation :type max_f_eval: int
:verbose: whether to show informations about the current restart :param max_iters: maximum number of iterations
:parallel: whether to run each restart as a separate process. It relies on the multiprocessing module. :type max_iters: int
:num_processes: number of workers in the multiprocessing pool :param messages: whether to display during optimisation
:type messages: bool
..Note: If num_processes is None, the number of workes in the multiprocessing pool is automatically ..Note: If num_processes is None, the number of workes in the multiprocessing pool is automatically
set to the number of processors on the current machine. set to the number of processors on the current machine.
@ -230,8 +249,13 @@ class Model(Parameterized):
self._set_params_transformed(initial_parameters) self._set_params_transformed(initial_parameters)
def ensure_default_constraints(self): def ensure_default_constraints(self):
""" """
Ensure that any variables which should clearly be positive have been constrained somehow. Ensure that any variables which should clearly be positive
have been constrained somehow. The method performs a regular
expression search on parameter names looking for the terms
'variance', 'lengthscale', 'precision' and 'kappa'. If any of
these terms are present in the name the parameter is
constrained positive.
""" """
positive_strings = ['variance', 'lengthscale', 'precision', 'kappa'] positive_strings = ['variance', 'lengthscale', 'precision', 'kappa']
# param_names = self._get_param_names() # param_names = self._get_param_names()
@ -246,11 +270,15 @@ class Model(Parameterized):
def objective_function(self, x): def objective_function(self, x):
""" """
The objective function passed to the optimizer. It combines the likelihood and the priors. The objective function passed to the optimizer. It combines
the likelihood and the priors.
Failures are handled robustly. The algorithm will try several times to Failures are handled robustly. The algorithm will try several times to
return the objective, and will raise the original exception if it return the objective, and will raise the original exception if it
the objective cannot be computed. the objective cannot be computed.
:param x: the parameters of the model.
:parameter type: np.array
""" """
try: try:
self._set_params_transformed(x) self._set_params_transformed(x)
@ -267,8 +295,11 @@ class Model(Parameterized):
Gets the gradients from the likelihood and the priors. Gets the gradients from the likelihood and the priors.
Failures are handled robustly. The algorithm will try several times to Failures are handled robustly. The algorithm will try several times to
return the objective, and will raise the original exception if it return the gradients, and will raise the original exception if it
the objective cannot be computed. the objective cannot be computed.
:param x: the parameters of the model.
:parameter type: np.array
""" """
try: try:
self._set_params_transformed(x) self._set_params_transformed(x)
@ -282,6 +313,13 @@ class Model(Parameterized):
return obj_grads return obj_grads
def objective_and_gradients(self, x): def objective_and_gradients(self, x):
"""
Compute the objective function of the model and the gradient of the model at the point given by x.
:param x: the point at which gradients are to be computed.
:type np.array:
"""
try: try:
self._set_params_transformed(x) self._set_params_transformed(x)
obj_f = -self.log_likelihood() - self.log_prior() obj_f = -self.log_likelihood() - self.log_prior()
@ -297,11 +335,13 @@ class Model(Parameterized):
def optimize(self, optimizer=None, start=None, **kwargs): def optimize(self, optimizer=None, start=None, **kwargs):
""" """
Optimize the Model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors. Optimize the model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors.
kwargs are passed to the optimizer. They can be: kwargs are passed to the optimizer. They can be:
:max_f_eval: maximum number of function evaluations :param max_f_eval: maximum number of function evaluations
:type max_f_eval: int
:messages: whether to display during optimisation :messages: whether to display during optimisation
:type messages: bool
:param optimzer: which optimizer to use (defaults to self.preferred optimizer) :param optimzer: which optimizer to use (defaults to self.preferred optimizer)
:type optimzer: string TODO: valid strings? :type optimzer: string TODO: valid strings?
""" """
@ -327,14 +367,14 @@ class Model(Parameterized):
self.optimization_runs.append(sgd) self.optimization_runs.append(sgd)
def Laplace_covariance(self): def Laplace_covariance(self):
"""return the covariance matric of a Laplace approximatino at the current (stationary) point""" """return the covariance matrix of a Laplace approximation at the current (stationary) point."""
# TODO add in the Prior contributions for MAP estimation # TODO add in the prior contributions for MAP estimation
# TODO fix the hessian for tied, constrained and fixed components # TODO fix the hessian for tied, constrained and fixed components
if hasattr(self, 'log_likelihood_hessian'): if hasattr(self, 'log_likelihood_hessian'):
A = -self.log_likelihood_hessian() A = -self.log_likelihood_hessian()
else: else:
print "numerically calculating hessian. please be patient!" print "numerically calculating Hessian. please be patient!"
x = self._get_params() x = self._get_params()
def f(x): def f(x):
self._set_params(x) self._set_params(x)
@ -348,8 +388,8 @@ class Model(Parameterized):
return A return A
def Laplace_evidence(self): def Laplace_evidence(self):
"""Returns an estiamte of the Model evidence based on the Laplace approximation. """Returns an estiamte of the model evidence based on the Laplace approximation.
Uses a numerical estimate of the hessian if none is available analytically""" Uses a numerical estimate of the Hessian if none is available analytically."""
A = self.Laplace_covariance() A = self.Laplace_covariance()
try: try:
hld = np.sum(np.log(np.diag(jitchol(A)[0]))) hld = np.sum(np.log(np.diag(jitchol(A)[0])))
@ -370,27 +410,28 @@ class Model(Parameterized):
log_prior = self.log_prior() log_prior = self.log_prior()
obj_funct = '\nLog-likelihood: {0:.3e}'.format(log_like) obj_funct = '\nLog-likelihood: {0:.3e}'.format(log_like)
if len(''.join(strs)) != 0: if len(''.join(strs)) != 0:
obj_funct += ', Log Prior: {0:.3e}, LL+Prior = {0:.3e}'.format(log_prior, log_like + log_prior) obj_funct += ', Log prior: {0:.3e}, LL+prior = {0:.3e}'.format(log_prior, log_like + log_prior)
obj_funct += '\n\n' obj_funct += '\n\n'
s[0] = obj_funct + s[0] s[0] = obj_funct + s[0]
s[0] += "|{h:^{col}}".format(h='Prior', col=width) s[0] += "|{h:^{col}}".format(h='prior', col=width)
s[1] += '-' * (width + 1) s[1] += '-' * (width + 1)
for p in range(2, len(strs) + 2): for p in range(2, len(strs) + 2):
s[p] += '|{Prior:^{width}}'.format(Prior=strs[p - 2], width=width) s[p] += '|{prior:^{width}}'.format(prior=strs[p - 2], width=width)
return '\n'.join(s) return '\n'.join(s)
def checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3): def checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3):
""" """
Check the gradient of the Model by comparing to a numerical estimate. Check the gradient of the ,odel by comparing to a numerical
If the verbose flag is passed, invividual components are tested (and printed) estimate. If the verbose flag is passed, invividual
components are tested (and printed)
:param verbose: If True, print a "full" checking of each parameter :param verbose: If True, print a "full" checking of each parameter
:type verbose: bool :type verbose: bool
:param step: The size of the step around which to linearise the objective :param step: The size of the step around which to linearise the objective
:type step: float (defaul 1e-6) :type step: float (default 1e-6)
:param tolerance: the tolerance allowed (see note) :param tolerance: the tolerance allowed (see note)
:type tolerance: float (default 1e-3) :type tolerance: float (default 1e-3)
@ -467,15 +508,15 @@ class Model(Parameterized):
def input_sensitivity(self): def input_sensitivity(self):
""" """
return an array describing the sesitivity of the Model to each input return an array describing the sesitivity of the model to each input
NB. Right now, we're basing this on the lengthscales (or NB. Right now, we're basing this on the lengthscales (or
variances) of the kernel. TODO: proper sensitivity analysis variances) of the kernel. TODO: proper sensitivity analysis
where we integrate across the Model inputs and evaluate the where we integrate across the model inputs and evaluate the
effect on the variance of the Model output. """ effect on the variance of the model output. """
if not hasattr(self, 'kern'): if not hasattr(self, 'kern'):
raise ValueError, "this Model has no kernel" raise ValueError, "this model has no kernel"
k = [p for p in self.kern.parts if p.name in ['rbf', 'linear', 'rbf_inv']] k = [p for p in self.kern.parts if p.name in ['rbf', 'linear', 'rbf_inv']]
if (not len(k) == 1) or (not k[0].ARD): if (not len(k) == 1) or (not k[0].ARD):

View file

@ -14,7 +14,10 @@ class kern(Parameterized):
""" """
This is the main kernel class for GPy. It handles multiple (additive) kernel functions, and keeps track of variaous things like which parameters live where. This is the main kernel class for GPy. It handles multiple (additive) kernel functions, and keeps track of variaous things like which parameters live where.
The technical code for kernels is divided into _parts_ (see e.g. rbf.py). This obnject contains a list of parts, which are computed additively. For multiplication, special _prod_ parts are used. The technical code for kernels is divided into _parts_ (see
e.g. rbf.py). This object contains a list of parts, which are
computed additively. For multiplication, special _prod_ parts
are used.
:param input_dim: The dimensionality of the kernel's input space :param input_dim: The dimensionality of the kernel's input space
:type input_dim: int :type input_dim: int
@ -149,7 +152,7 @@ class kern(Parameterized):
return g return g
def compute_param_slices(self): def compute_param_slices(self):
"""create a set of slices that can index the parameters of each part""" """create a set of slices that can index the parameters of each part."""
self.param_slices = [] self.param_slices = []
count = 0 count = 0
for p in self.parts: for p in self.parts:
@ -200,11 +203,19 @@ class kern(Parameterized):
""" """
return self.prod(other) return self.prod(other)
def __pow__(self, other, tensor=False):
"""
Shortcut for tensor `prod`.
"""
return self.prod(other, tensor=True)
def prod(self, other, tensor=False): def prod(self, other, tensor=False):
""" """
multiply two kernels (either on the same space, or on the tensor product of the input space) multiply two kernels (either on the same space, or on the tensor product of the input space).
:param other: the other kernel to be added :param other: the other kernel to be added
:type other: GPy.kern :type other: GPy.kern
:param tensor: whether or not to use the tensor space (default is false).
:type tensor: bool
""" """
K1 = self.copy() K1 = self.copy()
K2 = other.copy() K2 = other.copy()
@ -273,7 +284,7 @@ class kern(Parameterized):
[p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)] [p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)]
def _get_param_names(self): def _get_param_names(self):
# this is a bit nasty: we wat to distinguish between parts with the same name by appending a count # this is a bit nasty: we want to distinguish between parts with the same name by appending a count
part_names = np.array([k.name for k in self.parts], dtype=np.str) part_names = np.array([k.name for k in self.parts], dtype=np.str)
counts = [np.sum(part_names == ni) for i, ni in enumerate(part_names)] counts = [np.sum(part_names == ni) for i, ni in enumerate(part_names)]
cum_counts = [np.sum(part_names[i:] == ni) for i, ni in enumerate(part_names)] cum_counts = [np.sum(part_names[i:] == ni) for i, ni in enumerate(part_names)]
@ -295,11 +306,13 @@ class kern(Parameterized):
def dK_dtheta(self, dL_dK, X, X2=None): def dK_dtheta(self, dL_dK, X, X2=None):
""" """
:param dL_dK: An array of dL_dK derivaties, dL_dK Compute the gradient of the covariance function with respect to the parameters.
:type dL_dK: Np.ndarray (N x num_inducing)
:param dL_dK: An array of gradients of the objective function with respect to the covariance function.
:type dL_dK: Np.ndarray (num_samples x num_inducing)
:param X: Observed data inputs :param X: Observed data inputs
:type X: np.ndarray (N x input_dim) :type X: np.ndarray (num_samples x input_dim)
:param X2: Observed dara inputs (optional, defaults to X) :param X2: Observed data inputs (optional, defaults to X)
:type X2: np.ndarray (num_inducing x input_dim) :type X2: np.ndarray (num_inducing x input_dim)
""" """
assert X.shape[1] == self.input_dim assert X.shape[1] == self.input_dim
@ -312,6 +325,14 @@ class kern(Parameterized):
return self._transform_gradients(target) return self._transform_gradients(target)
def dK_dX(self, dL_dK, X, X2=None): def dK_dX(self, dL_dK, X, X2=None):
"""Compute the gradient of the covariance function with respect to X.
:param dL_dK: An array of gradients of the objective function with respect to the covariance function.
:type dL_dK: np.ndarray (num_samples x num_inducing)
:param X: Observed data inputs
:type X: np.ndarray (num_samples x input_dim)
:param X2: Observed data inputs (optional, defaults to X)
:type X2: np.ndarray (num_inducing x input_dim)"""
if X2 is None: if X2 is None:
X2 = X X2 = X
target = np.zeros_like(X) target = np.zeros_like(X)
@ -322,6 +343,7 @@ class kern(Parameterized):
return target return target
def Kdiag(self, X, which_parts='all'): def Kdiag(self, X, which_parts='all'):
"""Compute the diagonal of the covariance function for inputs X."""
if which_parts == 'all': if which_parts == 'all':
which_parts = [True] * self.Nparts which_parts = [True] * self.Nparts
assert X.shape[1] == self.input_dim assert X.shape[1] == self.input_dim
@ -330,6 +352,7 @@ class kern(Parameterized):
return target return target
def dKdiag_dtheta(self, dL_dKdiag, X): def dKdiag_dtheta(self, dL_dKdiag, X):
"""Compute the gradient of the diagonal of the covariance function with respect to the parameters."""
assert X.shape[1] == self.input_dim assert X.shape[1] == self.input_dim
assert dL_dKdiag.size == X.shape[0] assert dL_dKdiag.size == X.shape[0]
target = np.zeros(self.num_params) target = np.zeros(self.num_params)
@ -373,16 +396,18 @@ class kern(Parameterized):
return target return target
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S): def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
"""return shapes are N,num_inducing,input_dim""" """return shapes are num_samples,num_inducing,input_dim"""
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1])) target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
[p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] [p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
return target_mu, target_S return target_mu, target_S
def psi2(self, Z, mu, S): def psi2(self, Z, mu, S):
""" """
Computer the psi2 statistics for the covariance function.
:param Z: np.ndarray of inducing inputs (num_inducing x input_dim) :param Z: np.ndarray of inducing inputs (num_inducing x input_dim)
:param mu, S: np.ndarrays of means and variances (each N x input_dim) :param mu, S: np.ndarrays of means and variances (each num_samples x input_dim)
:returns psi2: np.ndarray (N,num_inducing,num_inducing) :returns psi2: np.ndarray (num_samples,num_inducing,num_inducing)
""" """
target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0])) target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
[p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)] [p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
@ -406,6 +431,7 @@ class kern(Parameterized):
return target return target
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S): def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S):
"""Gradient of the psi2 statistics with respect to the parameters."""
target = np.zeros(self.num_params) target = np.zeros(self.num_params)
[p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)] [p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)]
@ -509,3 +535,36 @@ class kern(Parameterized):
pb.title("k(x1,x2 ; %0.1f,%0.1f)" % (x[0, 0], x[0, 1])) pb.title("k(x1,x2 ; %0.1f,%0.1f)" % (x[0, 0], x[0, 1]))
else: else:
raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions" raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"
def objective_and_gradients_dK_dtheta(self, param, X, X2=None):
self._set_param(param)
K = self.K(X, X2)
f = K.sum()
dL_dK = np.ones_like(K)
g = self.dK_dtheta(param, dL_dK, X, X2)
return f, g
def objective_and_gradients_dK_dX(self, param, X, X2=None):
self._set_param(param)
K = self.K(X, X2)
f = K.sum()
dL_dK = np.ones_like(K)
g = self.dK_dX(param, dL_dK, X, X2)
return f, g
def objective_and_gradients_dKdiag_dtheta(self, param, X, X2=None):
self._set_param(param)
Kdiag = self.Kdiag(X)
f = Kdiag.sum()
dL_dK = np.ones_like(Kdiag)
g = self.dKdiag_dtheta(param, dL_dK, X)
return f, g
def objective_and_gradients_dKdiag_dX(self, param, X, X2=None):
self._set_param(param)
Kdiag = self.Kdiag(X)
f = Kdiag.sum()
dL_dK = np.ones_like(Kdiag)
g = self.dK_dX(param, dL_dK, X)
return f, g

View file

@ -23,3 +23,7 @@ Can we comment up some of the list incomprehensions in hierarchical.py??
Need to tidy up classification.py, Need to tidy up classification.py,
many examples include help that doesn't apply many examples include help that doesn't apply
(it is suggested that you can try different approximation types) (it is suggested that you can try different approximation types)
Shall we overload the ** operator to have tensor products?
People aren't filling the doc strings in as they go *everyone* needs to get in the habit of this (and modifying them as they edit, or correcting them when there is a problem).

View file

@ -42,8 +42,7 @@ class KernelTests(unittest.TestCase):
self.assertTrue(m.checkgrad()) self.assertTrue(m.checkgrad())
if __name__ == "__main__": if __name__ == "__main__":
print "Running unit tests, please be (very) patient..." print "Running unit tests, please be (very) patient..."
unittest.main() unittest.main()