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

Merge branch 'devel' of https://github.com/SheffieldML/GPy into devel

Browse source 
Conflicts:
	GPy/examples/dimensionality_reduction.py
This commit is contained in:
Neil Lawrence 2013-05-05 08:01:47 +01:00
commit 7ffcefc511
20 changed files with 981 additions and 409 deletions
Download patch file Download diff file Expand all files Collapse all files

37
GPy/core/model.py
View file

@ -6,7 +6,7 @@ from .. import likelihoods
from ..inference import optimization from ..inference import optimization
from ..util.linalg import jitchol from ..util.linalg import jitchol
from GPy.util.misc import opt_wrapper from GPy.util.misc import opt_wrapper
from parameterised import parameterised, truncate_pad from parameterised import parameterised
from scipy import optimize from scipy import optimize
import multiprocessing as mp import multiprocessing as mp
import numpy as np import numpy as np
@ -67,9 +67,14 @@ class model(parameterised):
# check constraints are okay # check constraints are okay
if isinstance(what, (priors.gamma, priors.log_Gaussian)): if isinstance(what, (priors.gamma, priors.log_Gaussian)):
assert not np.any(which[:, None] == self.constrained_negative_indices), "constraint and prior incompatible" constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain == 'positive']
assert not np.any(which[:, None] == self.constrained_bounded_indices), "constraint and prior incompatible" if len(constrained_positive_indices):
unconst = np.setdiff1d(which, self.constrained_positive_indices) constrained_positive_indices = np.hstack(constrained_positive_indices)
else:
constrained_positive_indices = np.zeros(shape=(0,))
bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices)
assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible"
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:"
print '\n'.join([n for i, n in enumerate(self._get_param_names()) if i in unconst]) print '\n'.join([n for i, n in enumerate(self._get_param_names()) if i in unconst])
@ -80,7 +85,6 @@ class model(parameterised):
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
@ -110,22 +114,16 @@ class model(parameterised):
return ret return ret
def _transform_gradients(self, g): def _transform_gradients(self, g):
"""
Takes a list of gradients and return an array of transformed gradients (positive/negative/tied/and so on)
"""
x = self._get_params() x = self._get_params()
g[self.constrained_positive_indices] = g[self.constrained_positive_indices] * x[self.constrained_positive_indices] for index, constraint in zip(self.constrained_indices, self.constraints):
g[self.constrained_negative_indices] = g[self.constrained_negative_indices] * x[self.constrained_negative_indices] g[index] = g[index] * constraint.gradfactor(x[index])
[np.put(g, i, g[i] * (x[i] - l) * (h - x[i]) / (h - l)) for i, l, h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]] [np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
if len(self.tied_indices) or len(self.constrained_fixed_indices): if len(self.tied_indices) or len(self.fixed_indices):
to_remove = np.hstack((self.constrained_fixed_indices + [t[1:] for t in self.tied_indices])) to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
return np.delete(g, to_remove) return np.delete(g, to_remove)
else: else:
return g return g
def randomize(self): def randomize(self):
""" """
Randomize the model. Randomize the model.
@ -209,7 +207,7 @@ class model(parameterised):
""" """
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.
""" """
positive_strings = ['variance', 'lengthscale', 'precision'] positive_strings = ['variance', 'lengthscale', 'precision', 'kappa']
param_names = self._get_param_names() param_names = self._get_param_names()
currently_constrained = self.all_constrained_indices() currently_constrained = self.all_constrained_indices()
to_make_positive = [] to_make_positive = []
@ -361,10 +359,7 @@ class model(parameterised):
numerical_gradient = (f1 - f2) / (2 * dx) numerical_gradient = (f1 - f2) / (2 * dx)
global_ratio = (f1 - f2) / (2 * np.dot(dx, gradient)) global_ratio = (f1 - f2) / (2 * np.dot(dx, gradient))
if (np.abs(1. - global_ratio) < tolerance) and not np.isnan(global_ratio): return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gradient - numerical_gradient).mean() - 1) < tolerance
return True
else:
return False
else: else:
# check the gradient of each parameter individually, and do some pretty printing # check the gradient of each parameter individually, and do some pretty printing
try: try:
@ -401,7 +396,7 @@ class model(parameterised):
ratio = (f1 - f2) / (2 * step * gradient) ratio = (f1 - f2) / (2 * step * gradient)
difference = np.abs((f1 - f2) / 2 / step - gradient) difference = np.abs((f1 - f2) / 2 / step - gradient)
if (np.abs(ratio - 1) < tolerance): if (np.abs(1. - ratio) < tolerance) or np.abs(difference) < tolerance:
formatted_name = "\033[92m {0} \033[0m".format(names[i]) formatted_name = "\033[92m {0} \033[0m".format(names[i])
else: else:
formatted_name = "\033[91m {0} \033[0m".format(names[i]) formatted_name = "\033[91m {0} \033[0m".format(names[i])

250
GPy/core/parameterised.py
View file

@ -9,26 +9,7 @@ import cPickle
import os import os
from ..util.squashers import sigmoid from ..util.squashers import sigmoid
import warnings import warnings
import transformations
def truncate_pad(string, width, align='m'):
"""
A helper function to make aligned strings for parameterised.__str__
"""
width = max(width, 4)
if len(string) > width:
return string[:width - 3] + '...'
elif len(string) == width:
return string
elif len(string) < width:
diff = width - len(string)
if align == 'm':
return ' ' * np.floor(diff / 2.) + string + ' ' * np.ceil(diff / 2.)
elif align == 'l':
return string + ' ' * diff
elif align == 'r':
return ' ' * diff + string
else:
raise ValueError
class parameterised(object): class parameterised(object):
def __init__(self): def __init__(self):
@ -36,13 +17,10 @@ class parameterised(object):
This is the base class for model and kernel. Mostly just handles tieing and constraining of parameters This is the base class for model and kernel. Mostly just handles tieing and constraining of parameters
""" """
self.tied_indices = [] self.tied_indices = []
self.constrained_fixed_indices = [] self.fixed_indices = []
self.constrained_fixed_values = [] self.fixed_values = []
self.constrained_positive_indices = np.empty(shape=(0,), dtype=np.int64) self.constrained_indices = []
self.constrained_negative_indices = np.empty(shape=(0,), dtype=np.int64) self.constraints = []
self.constrained_bounded_indices = []
self.constrained_bounded_uppers = []
self.constrained_bounded_lowers = []
def pickle(self, filename, protocol= -1): def pickle(self, filename, protocol= -1):
f = file(filename, 'w') f = file(filename, 'w')
@ -50,10 +28,7 @@ class parameterised(object):
f.close() f.close()
def copy(self): def copy(self):
""" """Returns a (deep) copy of the current model """
Returns a (deep) copy of the current model
"""
return copy.deepcopy(self) return copy.deepcopy(self)
@property @property
@ -64,6 +39,7 @@ class parameterised(object):
:see_also: :py:func:`GPy.core.parameterised.params_transformed` :see_also: :py:func:`GPy.core.parameterised.params_transformed`
""" """
return self._get_params() return self._get_params()
@params.setter @params.setter
def params(self, params): def params(self, params):
self._set_params(params) self._set_params(params)
@ -76,6 +52,7 @@ class parameterised(object):
:see_also: :py:func:`GPy.core.parameterised.params` :see_also: :py:func:`GPy.core.parameterised.params`
""" """
return self._get_params_transformed() return self._get_params_transformed()
@params_transformed.setter @params_transformed.setter
def params_transformed(self, params): def params_transformed(self, params):
self._set_params_transformed(params) self._set_params_transformed(params)
@ -97,7 +74,9 @@ class parameterised(object):
def __getitem__(self, name, return_names=False): def __getitem__(self, name, return_names=False):
""" """
Get a model parameter by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. Get a model parameter by name. The name is applied as a regular
expression and all parameters that match that regular expression are
returned.
""" """
matches = self.grep_param_names(name) matches = self.grep_param_names(name)
if len(matches): if len(matches):
@ -110,7 +89,9 @@ class parameterised(object):
def __setitem__(self, name, val): def __setitem__(self, name, val):
""" """
Set model parameter(s) by name. The name is provided as a regular expression. All parameters matching that regular expression are set to ghe given value. Set model parameter(s) by name. The name is provided as a regular
expression. All parameters matching that regular expression are set to
the given value.
""" """
matches = self.grep_param_names(name) matches = self.grep_param_names(name)
if len(matches): if len(matches):
@ -119,8 +100,6 @@ class parameterised(object):
x = self.params x = self.params
x[matches] = val x[matches] = val
self.params = x self.params = x
# import ipdb;ipdb.set_trace()
# self.params[matches] = val
else: else:
raise AttributeError, "no parameter matches %s" % name raise AttributeError, "no parameter matches %s" % name
@ -140,13 +119,6 @@ class parameterised(object):
"""Unties all parameters by setting tied_indices to an empty list.""" """Unties all parameters by setting tied_indices to an empty list."""
self.tied_indices = [] self.tied_indices = []
def all_constrained_indices(self):
"""Return a np array of all the constrained indices"""
ret = [np.hstack(i) for i in [self.constrained_bounded_indices, self.constrained_positive_indices, self.constrained_negative_indices, self.constrained_fixed_indices] if len(i)]
if len(ret):
return np.hstack(ret)
else:
return []
def grep_param_names(self, expr): def grep_param_names(self, expr):
""" """
Arguments Arguments
@ -159,7 +131,7 @@ class parameterised(object):
Notes Notes
----- -----
Other objects are passed through - i.e. integers which were'nt meant for grepping Other objects are passed through - i.e. integers which weren't meant for grepping
""" """
if type(expr) in [str, np.string_, np.str]: if type(expr) in [str, np.string_, np.str]:
@ -171,99 +143,76 @@ class parameterised(object):
return expr return expr
def Nparam_transformed(self): def Nparam_transformed(self):
ties = 0 removed = 0
for ar in self.tied_indices: for tie in self.tied_indices:
ties += ar.size - 1 removed += tie.size - 1
return self.Nparam - len(self.constrained_fixed_indices) - ties
def constrain_positive(self, which): for fix in self.fixed_indices:
""" removed += fix.size
Set positive constraints.
Arguments
---------
which -- np.array(dtype=int), or regular expression object or string
"""
matches = self.grep_param_names(which)
assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained"
self.constrained_positive_indices = np.hstack((self.constrained_positive_indices, matches))
# check to ensure constraint is in place
x = self._get_params()
for i, xx in enumerate(x):
if (xx < 0) & (i in matches):
x[i] = -xx
self._set_params(x)
return len(self._get_params()) - removed
def unconstrain(self, which): def unconstrain(self, which):
"""Unconstrain matching parameters. does not untie parameters""" """Unconstrain matching parameters. does not untie parameters"""
matches = self.grep_param_names(which) matches = self.grep_param_names(which)
# positive/negative
self.constrained_positive_indices = np.delete(self.constrained_positive_indices, np.nonzero(np.sum(self.constrained_positive_indices[:, None] == matches[None, :], 1))[0]) #tranformed contraints:
self.constrained_negative_indices = np.delete(self.constrained_negative_indices, np.nonzero(np.sum(self.constrained_negative_indices[:, None] == matches[None, :], 1))[0]) for match in matches:
# bounded self.constrained_indices = [i[i<>match] for i in self.constrained_indices]
if len(self.constrained_bounded_indices):
self.constrained_bounded_indices = [np.delete(a, np.nonzero(np.sum(a[:, None] == matches[None, :], 1))[0]) for a in self.constrained_bounded_indices] #remove empty constraints
if np.hstack(self.constrained_bounded_indices).size: tmp = zip(*[(i,t) for i,t in zip(self.constrained_indices,self.constraints) if len(i)])
self.constrained_bounded_uppers, self.constrained_bounded_lowers, self.constrained_bounded_indices = zip(*[(u, l, i) for u, l, i in zip(self.constrained_bounded_uppers, self.constrained_bounded_lowers, self.constrained_bounded_indices) if i.size])
self.constrained_bounded_uppers, self.constrained_bounded_lowers, self.constrained_bounded_indices = list(self.constrained_bounded_uppers), list(self.constrained_bounded_lowers), list(self.constrained_bounded_indices)
else:
self.constrained_bounded_uppers, self.constrained_bounded_lowers, self.constrained_bounded_indices = [], [], []
# fixed:
for i, indices in enumerate(self.constrained_fixed_indices):
self.constrained_fixed_indices[i] = np.delete(indices, np.nonzero(np.sum(indices[:, None] == matches[None, :], 1))[0])
# remove empty elements
tmp = [(i, v) for i, v in zip(self.constrained_fixed_indices, self.constrained_fixed_values) if len(i)]
if tmp: if tmp:
self.constrained_fixed_indices, self.constrained_fixed_values = zip(*tmp) self.constrained_indices, self.constraints = zip(*[(i,t) for i,t in zip(self.constrained_indices,self.constraints) if len(i)])
self.constrained_fixed_indices, self.constrained_fixed_values = list(self.constrained_fixed_indices), list(self.constrained_fixed_values) self.constrained_indices, self.constraints = list(self.constrained_indices), list(self.constraints)
# fixed:
self.fixed_values = [np.delete(values, np.nonzero(np.sum(indices[:, None] == matches[None, :], 1))[0]) for indices,values in zip(self.fixed_indices,self.fixed_values)]
self.fixed_indices = [np.delete(indices, np.nonzero(np.sum(indices[:, None] == matches[None, :], 1))[0]) for indices in self.fixed_indices]
# remove empty elements
tmp = [(i, v) for i, v in zip(self.fixed_indices, self.fixed_values) if len(i)]
if tmp:
self.fixed_indices, self.fixed_values = zip(*tmp)
self.fixed_indices, self.fixed_values = list(self.fixed_indices), list(self.fixed_values)
else: else:
self.constrained_fixed_indices, self.constrained_fixed_values = [], [] self.fixed_indices, self.fixed_values = [], []
def constrain_negative(self, which): def constrain_negative(self, which):
""" """ Set negative constraints. """
Set negative constraints. self.constrain(which, transformations.negative_exponent())
:param which: which variables to constrain
:type which: regular expression string
"""
matches = self.grep_param_names(which)
assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained"
self.constrained_negative_indices = np.hstack((self.constrained_negative_indices, matches))
# check to ensure constraint is in place
x = self._get_params()
for i, xx in enumerate(x):
if (xx > 0.) and (i in matches):
x[i] = -xx
self._set_params(x)
def constrain_positive(self, which):
""" Set positive constraints. """
self.constrain(which, transformations.logexp())
def constrain_bounded(self, which,lower, upper): def constrain_bounded(self, which,lower, upper):
"""Set bounded constraints. """ Set bounded constraints. """
self.constrain(which, transformations.logistic(lower, upper))
def all_constrained_indices(self):
if len(self.constrained_indices) or len(self.fixed_indices):
return np.hstack(self.constrained_indices + self.fixed_indices)
else:
return np.empty(shape=(0,))
def constrain(self,which,transform):
assert isinstance(transform,transformations.transformation)
Arguments
---------
which -- np.array(dtype=int), or regular expression object or string
upper -- (float) the upper bound on the constraint
lower -- (float) the lower bound on the constraint
"""
matches = self.grep_param_names(which) matches = self.grep_param_names(which)
assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained" overlap = set(matches).intersection(set(self.all_constrained_indices()))
assert lower < upper, "lower bound must be smaller than upper bound!" if overlap:
self.constrained_bounded_indices.append(matches) self.unconstrain(np.asarray(list(overlap)))
self.constrained_bounded_uppers.append(upper) print 'Warning: re-constraining these parameters'
self.constrained_bounded_lowers.append(lower) pn = self._get_param_names()
# check to ensure constraint is in place for i in overlap:
x = self._get_params() print pn[i]
for i, xx in enumerate(x):
if ((xx <= lower) | (xx >= upper)) & (i in matches):
x[i] = sigmoid(xx) * (upper - lower) + lower
self._set_params(x)
self.constrained_indices.append(matches)
self.constraints.append(transform)
x = self._get_params()
x[matches] = transform.initialize(x[matches])
self._set_params(x)
def constrain_fixed(self, which, value=None): def constrain_fixed(self, which, value=None):
""" """
@ -280,42 +229,36 @@ class parameterised(object):
""" """
matches = self.grep_param_names(which) matches = self.grep_param_names(which)
assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained" assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained"
self.constrained_fixed_indices.append(matches) self.fixed_indices.append(matches)
if value != None: if value != None:
self.constrained_fixed_values.append(value) self.fixed_values.append(value)
else: else:
self.constrained_fixed_values.append(self._get_params()[self.constrained_fixed_indices[-1]]) self.fixed_values.append(self._get_params()[self.fixed_indices[-1]])
# self.constrained_fixed_values.append(value) # self.fixed_values.append(value)
self._set_params_transformed(self._get_params_transformed()) self._set_params_transformed(self._get_params_transformed())
def _get_params_transformed(self): def _get_params_transformed(self):
"""use self._get_params to get the 'true' parameters of the model, which are then tied, constrained and fixed""" """use self._get_params to get the 'true' parameters of the model, which are then tied, constrained and fixed"""
x = self._get_params() x = self._get_params()
x[self.constrained_positive_indices] = np.log(x[self.constrained_positive_indices]) [np.put(x,i,t.finv(x[i])) for i,t in zip(self.constrained_indices,self.constraints)]
x[self.constrained_negative_indices] = np.log(-x[self.constrained_negative_indices])
[np.put(x, i, np.log(np.clip(x[i] - l, 1e-10, np.inf) / np.clip(h - x[i], 1e-10, np.inf))) for i, l, h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
to_remove = self.constrained_fixed_indices + [t[1:] for t in self.tied_indices] to_remove = self.fixed_indices + [t[1:] for t in self.tied_indices]
if len(to_remove): if len(to_remove):
return np.delete(x, np.hstack(to_remove)) return np.delete(x, np.hstack(to_remove))
else: else:
return x return x
def _set_params_transformed(self, x): def _set_params_transformed(self, x):
""" takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self._set_params""" """ takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self._set_params"""
# work out how many places are fixed, and where they are. tricky logic! # work out how many places are fixed, and where they are. tricky logic!
Nfix_places = 0. fix_places = self.fixed_indices + [t[1:] for t in self.tied_indices]
if len(self.tied_indices): if len(fix_places):
Nfix_places += np.hstack(self.tied_indices).size - len(self.tied_indices) fix_places = np.hstack(fix_places)
if len(self.constrained_fixed_indices): Nfix_places = fix_places.size
Nfix_places += np.hstack(self.constrained_fixed_indices).size
if Nfix_places:
fix_places = np.hstack(self.constrained_fixed_indices + [t[1:] for t in self.tied_indices])
else: else:
fix_places = [] Nfix_places = 0
free_places = np.setdiff1d(np.arange(Nfix_places + x.size, dtype=np.int), fix_places) free_places = np.setdiff1d(np.arange(Nfix_places + x.size, dtype=np.int), fix_places)
@ -323,11 +266,12 @@ class parameterised(object):
xx = np.zeros(Nfix_places + free_places.size, dtype=np.float64) xx = np.zeros(Nfix_places + free_places.size, dtype=np.float64)
xx[free_places] = x xx[free_places] = x
[np.put(xx, i, v) for i, v in zip(self.constrained_fixed_indices, self.constrained_fixed_values)] [np.put(xx, i, v) for i, v in zip(self.fixed_indices, self.fixed_values)]
[np.put(xx, i, v) for i, v in [(t[1:], xx[t[0]]) for t in self.tied_indices] ] [np.put(xx, i, v) for i, v in [(t[1:], xx[t[0]]) for t in self.tied_indices] ]
xx[self.constrained_positive_indices] = np.exp(xx[self.constrained_positive_indices])
xx[self.constrained_negative_indices] = -np.exp(xx[self.constrained_negative_indices]) [np.put(xx,i,t.f(xx[i])) for i,t in zip(self.constrained_indices, self.constraints)]
[np.put(xx, i, low + sigmoid(xx[i]) * (high - low)) for i, low, high in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)] if hasattr(self,'debug'):
stop
self._set_params(xx) self._set_params(xx)
def _get_param_names_transformed(self): def _get_param_names_transformed(self):
@ -346,17 +290,13 @@ class parameterised(object):
remove = np.empty(shape=(0,), dtype=np.int) remove = np.empty(shape=(0,), dtype=np.int)
# also remove the fixed params # also remove the fixed params
if len(self.constrained_fixed_indices): if len(self.fixed_indices):
remove = np.hstack((remove, np.hstack(self.constrained_fixed_indices))) remove = np.hstack((remove, np.hstack(self.fixed_indices)))
# add markers to show that some variables are constrained # add markers to show that some variables are constrained
for i in self.constrained_positive_indices: for i,t in zip(self.constrained_indices,self.constraints):
n[i] = n[i] + '(+ve)'
for i in self.constrained_negative_indices:
n[i] = n[i] + '(-ve)'
for i, l, h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers):
for ii in i: for ii in i:
n[ii] = n[ii] + '(bounded)' n[ii] = n[ii] + t.__str__()
n = [nn for i, nn in enumerate(n) if not i in remove] n = [nn for i, nn in enumerate(n) if not i in remove]
return n return n
@ -374,16 +314,12 @@ class parameterised(object):
values = self._get_params() # map(str,self._get_params()) values = self._get_params() # map(str,self._get_params())
# sort out the constraints # sort out the constraints
constraints = [''] * len(names) constraints = [''] * len(names)
for i in self.constrained_positive_indices: for i,t in zip(self.constrained_indices,self.constraints):
constraints[i] = '(+ve)' for ii in i:
for i in self.constrained_negative_indices: constraints[ii] = t.__str__()
constraints[i] = '(-ve)' for i in self.fixed_indices:
for i in self.constrained_fixed_indices:
for ii in i: for ii in i:
constraints[ii] = 'Fixed' constraints[ii] = 'Fixed'
for i, u, l in zip(self.constrained_bounded_indices, self.constrained_bounded_uppers, self.constrained_bounded_lowers):
for ii in i:
constraints[ii] = '(' + str(l) + ', ' + str(u) + ')'
# sort out the ties # sort out the ties
ties = [''] * len(names) ties = [''] * len(names)
for i, tie in enumerate(self.tied_indices): for i, tie in enumerate(self.tied_indices):

85
GPy/core/transformations.py Normal file
View file

@ -0,0 +1,85 @@
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
class transformation(object):
def __init__(self):
#set the domain. Suggest we use 'positive', 'bounded', etc
self.domain = 'undefined'
def f(self, x):
raise NotImplementedError
def finv(self,x):
raise NotImplementedError
def gradfactor(self,f):
""" df_dx evaluated at self.f(x)=f"""
raise NotImplementedError
def initialize(self,f):
""" produce a sensible initial values for f(x)"""
raise NotImplementedError
def __str__(self):
raise NotImplementedError
class logexp(transformation):
def __init__(self):
self.domain= 'positive'
def f(self,x):
return np.log(1. + np.exp(x))
def finv(self,f):
return np.log(np.exp(f) - 1.)
def gradfactor(self,f):
ef = np.exp(f)
return (ef - 1.)/ef
def initialize(self,f):
return np.abs(f)
def __str__(self):
return '(+ve)'
class exponent(transformation):
def __init__(self):
self.domain= 'positive'
def f(self,x):
return np.exp(x)
def finv(self,x):
return np.log(x)
def gradfactor(self,f):
return f
def initialize(self,f):
return np.abs(f)
def __str__(self):
return '(+ve)'
class negative_exponent(transformation):
def __init__(self):
self.domain= 'negative'
def f(self,x):
return -np.exp(x)
def finv(self,x):
return np.log(-x)
def gradfactor(self,f):
return f
def initialize(self,f):
return -np.abs(f)
def __str__(self):
return '(-ve)'
class logistic(transformation):
def __init__(self,lower,upper):
self.domain= 'bounded'
assert lower < upper
self.lower, self.upper = float(lower), float(upper)
self.difference = self.upper - self.lower
def f(self,x):
return self.lower + self.difference/(1.+np.exp(-x))
def finv(self,f):
return np.log(np.clip(f - self.lower, 1e-10, np.inf) / np.clip(self.upper - f, 1e-10, np.inf))
def gradfactor(self,f):
return (f-self.lower)*(self.upper-f)/self.difference
def initialize(self,f):
return self.f(f*0.)
def __str__(self):
return '({},{})'.format(self.lower,self.upper)

19
GPy/examples/dimensionality_reduction.py
View file

@ -82,7 +82,7 @@ def BGPLVM_oil(optimize=True, N=100, Q=10, M=15, max_f_eval=300):
m.ensure_default_constraints() m.ensure_default_constraints()
y = m.likelihood.Y[0, :] y = m.likelihood.Y[0, :]
fig,(latent_axes,sense_axes) = plt.subplots(1,2) fig, (latent_axes, hist_axes) = plt.subplots(1, 2)
plt.sca(latent_axes) plt.sca(latent_axes)
m.plot_latent() m.plot_latent()
data_show = GPy.util.visualize.vector_show(y) data_show = GPy.util.visualize.vector_show(y)
@ -181,15 +181,26 @@ def bgplvm_simulation_matlab_compare():
from GPy.models import mrd from GPy.models import mrd
from GPy import kern from GPy import kern
reload(mrd); reload(kern) reload(mrd); reload(kern)
k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k,
# X=mu, # X=mu,
# X_variance=S, # X_variance=S,
_debug=True) _debug=True)
m.ensure_default_constraints() m.ensure_default_constraints()
m.auto_scale_factor = True m.auto_scale_factor = True
m['noise'] = .01 # Y.var() / 100. m['noise'] = Y.var() / 100.
m['{}_variance'.format(k.parts[0].name)] = .01 m['linear_variance'] = .01
# lscstr = 'X_variance'
# m[lscstr] = .01
# m.unconstrain(lscstr); m.constrain_fixed(lscstr, .1)
# cstr = 'white'
# m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.)
# cstr = 'noise'
# m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.)
return m return m
def bgplvm_simulation(burnin='scg', plot_sim=False, def bgplvm_simulation(burnin='scg', plot_sim=False,

275
GPy/inference/conjugate_gradient_descent.py Normal file
View file

@ -0,0 +1,275 @@
'''
Created on 24 Apr 2013
@author: maxz
'''
from GPy.inference.gradient_descent_update_rules import FletcherReeves
from Queue import Empty
from multiprocessing import Value
from multiprocessing.queues import Queue
from multiprocessing.synchronize import Event
from scipy.optimize.linesearch import line_search_wolfe1, line_search_wolfe2
from threading import Thread
import numpy
import sys
RUNNING = "running"
CONVERGED = "converged"
MAXITER = "maximum number of iterations reached"
MAX_F_EVAL = "maximum number of function calls reached"
LINE_SEARCH = "line search failed"
KBINTERRUPT = "interrupted"
class _Async_Optimization(Thread):
def __init__(self, f, df, x0, update_rule, runsignal, SENTINEL,
report_every=10, messages=0, maxiter=5e3, max_f_eval=15e3,
gtol=1e-6, outqueue=None, *args, **kw):
"""
Helper Process class for async optimization
f_call and df_call are Multiprocessing Values, for synchronized assignment
"""
self.f_call = Value('i', 0)
self.df_call = Value('i', 0)
self.f = self.f_wrapper(f, self.f_call)
self.df = self.f_wrapper(df, self.df_call)
self.x0 = x0
self.update_rule = update_rule
self.report_every = report_every
self.messages = messages
self.maxiter = maxiter
self.max_f_eval = max_f_eval
self.gtol = gtol
self.SENTINEL = SENTINEL
self.runsignal = runsignal
# self.parent = parent
# self.result = None
self.outq = outqueue
super(_Async_Optimization, self).__init__(target=self.run,
name="CG Optimization",
*args, **kw)
# def __enter__(self):
# return self
#
# def __exit__(self, type, value, traceback):
# return isinstance(value, TypeError)
def f_wrapper(self, f, counter):
def f_w(*a, **kw):
counter.value += 1
return f(*a, **kw)
return f_w
def callback(self, *a):
if self.outq is not None:
self.outq.put(a)
# self.parent and self.parent.callback(*a, **kw)
pass
# print "callback done"
def callback_return(self, *a):
self.callback(*a)
self.callback(self.SENTINEL)
self.runsignal.clear()
def run(self, *args, **kwargs):
raise NotImplementedError("Overwrite this with optimization (for async use)")
pass
class _CGDAsync(_Async_Optimization):
def reset(self, xi, *a, **kw):
gi = -self.df(xi, *a, **kw)
si = gi
ur = self.update_rule(gi)
return gi, ur, si
def run(self, *a, **kw):
status = RUNNING
fi = self.f(self.x0)
fi_old = fi + 5000
gi, ur, si = self.reset(self.x0, *a, **kw)
xi = self.x0
xi_old = numpy.nan
it = 0
while it < self.maxiter:
if not self.runsignal.is_set():
break
if self.f_call.value > self.max_f_eval:
status = MAX_F_EVAL
gi = -self.df(xi, *a, **kw)
if numpy.dot(gi.T, gi) < self.gtol:
status = CONVERGED
break
if numpy.isnan(numpy.dot(gi.T, gi)):
if numpy.any(numpy.isnan(xi_old)):
status = CONVERGED
break
self.reset(xi_old)
gammai = ur(gi)
if gammai < 1e-6 or it % xi.shape[0] == 0:
gi, ur, si = self.reset(xi, *a, **kw)
si = gi + gammai * si
alphai, _, _, fi2, fi_old2, gfi = line_search_wolfe1(self.f,
self.df,
xi,
si, gi,
fi, fi_old)
if alphai is not None and fi2 < fi:
fi, fi_old = fi2, fi_old2
else:
alphai, _, _, fi, fi_old, gfi = \
line_search_wolfe2(self.f, self.df,
xi, si, gi,
fi, fi_old)
if alphai is None:
# This line search also failed to find a better solution.
status = LINE_SEARCH
break
if gfi is not None:
gi = gfi
if numpy.isnan(fi) or fi_old < fi:
gi, ur, si = self.reset(xi, *a, **kw)
else:
xi += numpy.dot(alphai, si)
if self.messages:
sys.stdout.write("\r")
sys.stdout.flush()
sys.stdout.write("iteration: {0:> 6g} f:{1:> 12e} |g|:{2:> 12e}".format(it, fi, numpy.dot(gi.T, gi)))
if it % self.report_every == 0:
self.callback(xi, fi, gi, it, self.f_call.value, self.df_call.value, status)
it += 1
else:
status = MAXITER
self.callback_return(xi, fi, gi, it, self.f_call.value, self.df_call.value, status)
self.result = [xi, fi, gi, it, self.f_call.value, self.df_call.value, status]
class Async_Optimize(object):
callback = lambda *x: None
runsignal = Event()
SENTINEL = "SENTINEL"
def async_callback_collect(self, q):
while self.runsignal.is_set():
try:
for ret in iter(lambda: q.get(timeout=1), self.SENTINEL):
self.callback(*ret)
except Empty:
pass
def opt_async(self, f, df, x0, callback, update_rule=FletcherReeves,
messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6,
report_every=10, *args, **kwargs):
self.runsignal.set()
c = None
outqueue = None
if callback:
outqueue = Queue()
self.callback = callback
c = Thread(target=self.async_callback_collect, args=(outqueue,))
c.start()
p = _CGDAsync(f, df, x0, update_rule, self.runsignal, self.SENTINEL,
report_every=report_every, messages=messages, maxiter=maxiter,
max_f_eval=max_f_eval, gtol=gtol, outqueue=outqueue, *args, **kwargs)
p.start()
return p, c
def opt(self, f, df, x0, callback=None, update_rule=FletcherReeves,
messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6,
report_every=10, *args, **kwargs):
p, c = self.opt_async(f, df, x0, callback, update_rule, messages,
maxiter, max_f_eval, gtol,
report_every, *args, **kwargs)
while self.runsignal.is_set():
try:
p.join(1)
# c.join(1)
except KeyboardInterrupt:
# print "^C"
self.runsignal.clear()
p.join()
c.join()
if c and c.is_alive():
print "WARNING: callback still running, optimisation done!"
return p.result
class CGD(Async_Optimize):
'''
Conjugate gradient descent algorithm to minimize
function f with gradients df, starting at x0
with update rule update_rule
if df returns tuple (grad, natgrad) it will optimize according
to natural gradient rules
'''
opt_name = "Conjugate Gradient Descent"
def opt_async(self, *a, **kw):
"""
opt_async(self, f, df, x0, callback, update_rule=FletcherReeves,
messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6,
report_every=10, *args, **kwargs)
callback gets called every `report_every` iterations
callback(xi, fi, gi, iteration, function_calls, gradient_calls, status_message)
if df returns tuple (grad, natgrad) it will optimize according
to natural gradient rules
f, and df will be called with
f(xi, *args, **kwargs)
df(xi, *args, **kwargs)
**returns**
-----------
Started `Process` object, optimizing asynchronously
**calls**
---------
callback(x_opt, f_opt, g_opt, iteration, function_calls, gradient_calls, status_message)
at end of optimization!
"""
return super(CGD, self).opt_async(*a, **kw)
def opt(self, *a, **kw):
"""
opt(self, f, df, x0, callback=None, update_rule=FletcherReeves,
messages=0, maxiter=5e3, max_f_eval=15e3, gtol=1e-6,
report_every=10, *args, **kwargs)
Minimize f, calling callback every `report_every` iterations with following syntax:
callback(xi, fi, gi, iteration, function_calls, gradient_calls, status_message)
if df returns tuple (grad, natgrad) it will optimize according
to natural gradient rules
f, and df will be called with
f(xi, *args, **kwargs)
df(xi, *args, **kwargs)
**returns**
---------
x_opt, f_opt, g_opt, iteration, function_calls, gradient_calls, status_message
at end of optimization
"""
return super(CGD, self).opt(*a, **kw)

43
GPy/inference/gradient_descent_update_rules.py Normal file
View file

@ -0,0 +1,43 @@
'''
Created on 24 Apr 2013
@author: maxz
'''
import numpy
class GDUpdateRule():
_gradnat = None
_gradnatold = None
def __init__(self, initgrad, initgradnat=None):
self.grad = initgrad
if initgradnat:
self.gradnat = initgradnat
else:
self.gradnat = initgrad
# self.grad, self.gradnat
def _gamma(self):
raise NotImplemented("""Implement gamma update rule here,
you can use self.grad and self.gradold for parameters, as well as
self.gradnat and self.gradnatold for natural gradients.""")
def __call__(self, grad, gradnat=None, si=None, *args, **kw):
"""
Return gamma for given gradients and optional natural gradients
"""
if not gradnat:
gradnat = grad
self.gradold = self.grad
self.gradnatold = self.gradnat
self.grad = grad
self.gradnat = gradnat
self.si = si
return self._gamma(*args, **kw)
class FletcherReeves(GDUpdateRule):
'''
Fletcher Reeves update rule for gamma
'''
def _gamma(self, *a, **kw):
tmp = numpy.dot(self.grad.T, self.gradnat)
if tmp:
return tmp / numpy.dot(self.gradold.T, self.gradnatold)
return tmp

58
GPy/kern/kern.py
View file

@ -48,11 +48,7 @@ class kern(parameterised):
def plot_ARD(self, ax=None): def plot_ARD(self, ax=None):
""" """If an ARD kernel is present, it bar-plots the ARD parameters"""
If an ARD kernel is present, it bar-plots the ARD parameters
"""
if ax is None: if ax is None:
ax = pb.gca() ax = pb.gca()
for p in self.parts: for p in self.parts:
@ -71,12 +67,10 @@ class kern(parameterised):
def _transform_gradients(self, g): def _transform_gradients(self, g):
x = self._get_params() x = self._get_params()
g[self.constrained_positive_indices] = g[self.constrained_positive_indices] * x[self.constrained_positive_indices] [np.put(x,i,x*t.gradfactor(x[i])) for i,t in zip(self.constrained_indices, self.constraints)]
g[self.constrained_negative_indices] = g[self.constrained_negative_indices] * x[self.constrained_negative_indices]
[np.put(g, i, g[i] * (x[i] - l) * (h - x[i]) / (h - l)) for i, l, h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]] [np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
if len(self.tied_indices) or len(self.constrained_fixed_indices): if len(self.tied_indices) or len(self.fixed_indices):
to_remove = np.hstack((self.constrained_fixed_indices + [t[1:] for t in self.tied_indices])) to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
return np.delete(g, to_remove) return np.delete(g, to_remove)
else: else:
return g return g
@ -93,13 +87,10 @@ class kern(parameterised):
assert self.D == other.D assert self.D == other.D
newkern = kern(self.D, self.parts + other.parts, self.input_slices + other.input_slices) newkern = kern(self.D, self.parts + other.parts, self.input_slices + other.input_slices)
# transfer constraints: # transfer constraints:
newkern.constrained_positive_indices = np.hstack((self.constrained_positive_indices, self.Nparam + other.constrained_positive_indices)) newkern.constrained_indices = self.constrained_indices + [i+self.Nparam for i in other.constrained_indices]
newkern.constrained_negative_indices = np.hstack((self.constrained_negative_indices, self.Nparam + other.constrained_negative_indices)) newkern.constraints = self.constraints + other.constraints
newkern.constrained_bounded_indices = self.constrained_bounded_indices + [self.Nparam + x for x in other.constrained_bounded_indices] newkern.fixed_indices = self.fixed_indices + [self.Nparam + x for x in other.fixed_indices]
newkern.constrained_bounded_lowers = self.constrained_bounded_lowers + other.constrained_bounded_lowers newkern.fixed_values = self.fixed_values + other.fixed_values
newkern.constrained_bounded_uppers = self.constrained_bounded_uppers + other.constrained_bounded_uppers
newkern.constrained_fixed_indices = self.constrained_fixed_indices + [self.Nparam + x for x in other.constrained_fixed_indices]
newkern.constrained_fixed_values = self.constrained_fixed_values + other.constrained_fixed_values
newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices] newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices]
return newkern return newkern
@ -126,13 +117,12 @@ class kern(parameterised):
newkern = kern(D, self.parts + other.parts, self_input_slices + other_input_slices) newkern = kern(D, self.parts + other.parts, self_input_slices + other_input_slices)
# transfer constraints: # transfer constraints:
newkern.constrained_positive_indices = np.hstack((self.constrained_positive_indices, self.Nparam + other.constrained_positive_indices)) newkern.constrained_indices = self.constrained_indices + [x+self.Nparam for x in other.constrained_indices]
newkern.constrained_negative_indices = np.hstack((self.constrained_negative_indices, self.Nparam + other.constrained_negative_indices)) newkern.constraints = self.constraints + other.constraints
newkern.constrained_bounded_indices = self.constrained_bounded_indices + [self.Nparam + x for x in other.constrained_bounded_indices] newkern.fixed_indices = self.fixed_indices + [self.Nparam + x for x in other.fixed_indices]
newkern.constrained_bounded_lowers = self.constrained_bounded_lowers + other.constrained_bounded_lowers newkern.fixed_values = self.fixed_values + other.fixed_values
newkern.constraints = self.constraints + other.constraints
newkern.constrained_bounded_uppers = self.constrained_bounded_uppers + other.constrained_bounded_uppers newkern.constrained_bounded_uppers = self.constrained_bounded_uppers + other.constrained_bounded_uppers
newkern.constrained_fixed_indices = self.constrained_fixed_indices + [self.Nparam + x for x in other.constrained_fixed_indices]
newkern.constrained_fixed_values = self.constrained_fixed_values + other.constrained_fixed_values
newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices] newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices]
return newkern return newkern
@ -208,15 +198,11 @@ class kern(parameterised):
# Get the ties and constrains of the kernels before the multiplication # 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_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_ind = [K1.constrained_indices] + [K1.Nparam + i for i in K2.constrained_indices]
prev_constr_neg = np.append(K1.constrained_negative_indices, K1.Nparam + K2.constrained_negative_indices) prev_constr = K1.constraints + K2.constraints
prev_constr_fix = K1.constrained_fixed_indices + [arr + K1.Nparam for arr in K2.constrained_fixed_indices] prev_constr_fix = K1.fixed_indices + [arr + K1.Nparam for arr in K2.fixed_indices]
prev_constr_fix_values = K1.constrained_fixed_values + K2.constrained_fixed_values prev_constr_fix_values = K1.fixed_values + K2.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:
@ -228,14 +214,8 @@ class kern(parameterised):
index = np.where(index_param == i)[0] index = np.where(index_param == i)[0]
if index.size > 1: if index.size > 1:
self.tie_params(index) self.tie_params(index)
for i in prev_constr_pos: for i,t in zip(prev_constr_ind,prev_constr):
self.constrain_positive(np.where(index_param == i)[0]) self.constrain(np.where(index_param == i)[0],t)
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])

70
GPy/kern/linear.py
View file

@ -140,13 +140,27 @@ class linear(kernpart):
returns N,M,M matrix returns N,M,M matrix
""" """
self._psi_computations(Z, mu, S) self._psi_computations(Z, mu, S)
#psi2 = self.ZZ*np.square(self.variances)*self.mu2_S[:, None, None, :] # psi2_old = self.ZZ * np.square(self.variances) * self.mu2_S[:, None, None, :]
# target += psi2.sum(-1) # target += psi2.sum(-1)
target += np.tensordot(self.ZZ[None,:,:,:]*np.square(self.variances),self.mu2_S[:, None, None, :],((3),(3))).squeeze().T # slow way of doing it, but right
# psi2_real = rm np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
# for n in range(mu.shape[0]):
# for m_prime in range(Z.shape[0]):
# for m in range(Z.shape[0]):
# tmp = self._Z[m:m + 1] * self.variances
# tmp = np.dot(tmp, (tdot(self._mu[n:n + 1].T) + np.diag(S[n])))
# psi2_real[n, m, m_prime] = np.dot(tmp, (
# self._Z[m_prime:m_prime + 1] * self.variances).T)
# mu2_S = (self._mu[:, None, :] * self._mu[:, :, None])
# mu2_S[:, np.arange(self.D), np.arange(self.D)] += self._S
# psi2 = (self.ZA[None, :, None, :] * mu2_S[:, None]).sum(-1)
# psi2 = (psi2[:, :, None] * self.ZA[None, None]).sum(-1)
# psi2_tensor = np.tensordot(self.ZZ[None, :, :, :] * np.square(self.variances), self.mu2_S[:, None, None, :], ((3), (3))).squeeze().T
target += self._psi2
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target): def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
self._psi_computations(Z, mu, S) self._psi_computations(Z, mu, S)
tmp = (dL_dpsi2[:,:,:,None]*(2.*self.ZZ*self.mu2_S[:,None,None,:]*self.variances)) tmp = dL_dpsi2[:, :, :, None] * (self.ZAinner[:, :, None, :] * (2 * Z)[None, None, :, :])
if self.ARD: if self.ARD:
target += tmp.sum(0).sum(0).sum(0) target += tmp.sum(0).sum(0).sum(0)
else: else:
@ -155,15 +169,35 @@ class linear(kernpart):
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S): def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
"""Think N,M,M,Q """ """Think N,M,M,Q """
self._psi_computations(Z, mu, S) self._psi_computations(Z, mu, S)
tmp = self.ZZ*np.square(self.variances) # M,M,Q AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]
target_mu += (dL_dpsi2[:,:,:,None]*tmp*2.*mu[:,None,None,:]).sum(1).sum(1) AZZA = AZZA + AZZA.swapaxes(1, 2)
target_S += (dL_dpsi2[:,:,:,None]*tmp).sum(1).sum(1) target_S += (dL_dpsi2[:, :, :, None] * self.ZA[None, :, None, :] * self.ZA[None, None, :, :]).sum(1).sum(1)
dpsi2_dmu = (dL_dpsi2[:, :, :, None] * np.tensordot(mu, AZZA, (-1, 0))).sum(1).sum(1)
target_mu += dpsi2_dmu
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target): def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target):
self._psi_computations(Z, mu, S) self._psi_computations(Z, mu, S)
mu2_S = np.sum(self.mu2_S,0)# Q, # mu2_S = np.sum(self.mu2_S, 0) # Q,
target += (dL_dpsi2[:,:,:,None] * (self.mu2_S[:,None,None,:]*(Z*np.square(self.variances)[None,:])[None,None,:,:])).sum(0).sum(1) # import ipdb;ipdb.set_trace()
#TODO: tensordot would gain some time here # psi2_dZ_real = np.zeros((mu.shape[0], Z.shape[0], Z.shape[1]))
# for n in range(mu.shape[0]):
# for m in range(Z.shape[0]):
# tmp = self.variances * (tdot(self._mu[n:n + 1].T) + np.diag(S[n]))
# psi2_dZ_real[n, m, :] = np.dot(tmp, (
# self._Z[m:m + 1] * self.variances).T).T
# tmp = self._Z[m:m + 1] * self.variances
# tmp = np.dot(tmp, (tdot(self._mu[n:n + 1].T) + np.diag(S[n])))
# psi2_dZ_real[n, m, :] = tmp * self.variances
# for m_prime in range(Z.shape[0]):
# if m == m_prime:
# psi2_dZ_real[n, m, :] *= 2
# prod = (dL_dpsi2[:, :, :, None] * np.eye(Z.shape[0])[None, :, :, None] * (self.ZAinner * self.variances).swapaxes(0, 1)[:, :, None, :])
# psi2_dZ = prod.swapaxes(1, 2) + prod
psi2_dZ = dL_dpsi2[:, :, :, None] * self.variances * self.ZAinner[:, :, None, :]
target += psi2_dZ.sum(0).sum(0)
# import ipdb;ipdb.set_trace()
# psi2_dZ_old = (dL_dpsi2[:, :, :, None] * (self.mu2_S[:, None, None, :] * (Z * np.square(self.variances)[None, :])[None, None, :, :])).sum(0).sum(1)
# target += (dL_dpsi2[:, :, :, None] * psi2_dZ_real[:, :, None, :]).sum(0).sum(0) * 2 # (self.variances * np.dot(self.inner, self.ZA.T)).sum(1)
#---------------------------------------# #---------------------------------------#
# Precomputations # # Precomputations #
@ -181,12 +215,22 @@ class linear(kernpart):
def _psi_computations(self, Z, mu, S): def _psi_computations(self, Z, mu, S):
# here are the "statistics" for psi1 and psi2 # here are the "statistics" for psi1 and psi2
if not np.all(Z==self._Z): Zv_changed = not (np.array_equal(Z, self._Z) and np.array_equal(self.variances, self._variances))
muS_changed = not (np.array_equal(mu, self._mu) and np.array_equal(S, self._S))
if Zv_changed:
# Z has changed, compute Z specific stuff # Z has changed, compute Z specific stuff
# self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,Q # self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,Q
self.ZZ = np.empty((Z.shape[0],Z.shape[0],Z.shape[1]),order='F') # self.ZZ = np.empty((Z.shape[0], Z.shape[0], Z.shape[1]), order='F')
[tdot(Z[:,i:i+1],self.ZZ[:,:,i].T) for i in xrange(Z.shape[1])] # [tdot(Z[:, i:i + 1], self.ZZ[:, :, i].T) for i in xrange(Z.shape[1])]
self.ZA = Z * self.variances
self._Z = Z.copy() self._Z = Z.copy()
if not (np.all(mu==self._mu) and np.all(S==self._S)): self._variances = self.variances.copy()
if muS_changed:
self.mu2_S = np.square(mu) + S self.mu2_S = np.square(mu) + S
self.inner = (mu[:, None, :] * mu[:, :, None])
diag_indices = np.diag_indices(mu.shape[1], 2)
self.inner[:, diag_indices[0], diag_indices[1]] += S
self._mu, self._S = mu.copy(), S.copy() self._mu, self._S = mu.copy(), S.copy()
if Zv_changed or muS_changed:
self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [M x N x Q]!
self._psi2 = np.dot(self.ZAinner, self.ZA.T)

4
GPy/kern/rbf.py
View file

@ -188,12 +188,12 @@ class rbf(kernpart):
self._X2 = None self._X2 = None
X = X/self.lengthscale X = X/self.lengthscale
Xsquare = np.sum(np.square(X),1) Xsquare = np.sum(np.square(X),1)
self._K_dist2 = (-2.*tdot(X) + Xsquare[:,None] + Xsquare[None,:]) self._K_dist2 = -2.*tdot(X) + (Xsquare[:,None] + Xsquare[None,:])
else: else:
self._X2 = X2.copy() self._X2 = X2.copy()
X = X/self.lengthscale X = X/self.lengthscale
X2 = X2/self.lengthscale X2 = X2/self.lengthscale
self._K_dist2 = (-2.*np.dot(X, X2.T) + np.sum(np.square(X),1)[:,None] + np.sum(np.square(X2),1)[None,:]) self._K_dist2 = -2.*np.dot(X, X2.T) + (np.sum(np.square(X),1)[:,None] + np.sum(np.square(X2),1)[None,:])
self._K_dvar = np.exp(-0.5*self._K_dist2) self._K_dvar = np.exp(-0.5*self._K_dist2)
def _psi_computations(self,Z,mu,S): def _psi_computations(self,Z,mu,S):

5
GPy/likelihoods/EP.py
View file

@ -196,8 +196,9 @@ class EP(likelihood):
self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
self.v_tilde[i] = self.v_tilde[i] + Delta_v self.v_tilde[i] = self.v_tilde[i] + Delta_v
#Posterior distribution parameters update #Posterior distribution parameters update
LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
L = jitchol(LLT) #L = jitchol(LLT)
cholupdate(L,Kmn[:,i]*np.sqrt(Delta_tau))
V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1) V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1)
Sigma_diag = np.sum(V*V,-2) Sigma_diag = np.sum(V*V,-2)
si = np.sum(V.T*V[:,i],-1) si = np.sum(V.T*V[:,i],-1)

8
GPy/likelihoods/likelihood_functions.py
View file

@ -53,9 +53,11 @@ class probit(likelihood_function):
mu = mu.flatten() mu = mu.flatten()
var = var.flatten() var = var.flatten()
mean = stats.norm.cdf(mu/np.sqrt(1+var)) mean = stats.norm.cdf(mu/np.sqrt(1+var))
p_025 = np.zeros(mu.shape) norm_025 = [stats.norm.ppf(.025,m,v) for m,v in zip(mu,var)]
p_975 = np.ones(mu.shape) norm_975 = [stats.norm.ppf(.975,m,v) for m,v in zip(mu,var)]
return mean, np.nan*var, p_025, p_975 # TODO: better values here (mean is okay) p_025 = stats.norm.cdf(norm_025/np.sqrt(1+var))
p_975 = stats.norm.cdf(norm_975/np.sqrt(1+var))
return mean, np.nan*var, p_025, p_975 # TODO: var
class Poisson(likelihood_function): class Poisson(likelihood_function):
""" """

28
GPy/models/Bayesian_GPLVM.py
View file

@ -54,6 +54,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
self._savedgradients = [] self._savedgradients = []
self._savederrors = [] self._savederrors = []
self._savedpsiKmm = [] self._savedpsiKmm = []
sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_variance=X_variance, **kwargs) sparse_GP.__init__(self, X, Gaussian(Y), kernel, Z=Z, X_variance=X_variance, **kwargs)
@property @property
@ -95,9 +96,9 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
print "\rWARNING: Caught LinAlgError, continueing without setting " print "\rWARNING: Caught LinAlgError, continueing without setting "
if self._debug: if self._debug:
self._savederrors.append(self.f_call) self._savederrors.append(self.f_call)
# if save_count > 10: if save_count > 10:
# raise raise
# self._set_params(self.oldps[-1], save_old=False, save_count=save_count + 1) self._set_params(self.oldps[-1], save_old=False, save_count=save_count + 1)
def dKL_dmuS(self): def dKL_dmuS(self):
dKL_dS = (1. - (1. / (self.X_variance))) * 0.5 dKL_dS = (1. - (1. / (self.X_variance))) * 0.5
@ -258,28 +259,28 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
ax2.text(.5, .5, r"$\mathbf{X}$", alpha=.5, transform=ax2.transAxes, ax2.text(.5, .5, r"$\mathbf{X}$", alpha=.5, transform=ax2.transAxes,
ha='center', va='center') ha='center', va='center')
figs[-1].canvas.draw() figs[-1].canvas.draw()
figs[-1].tight_layout(rect=(0, 0, 1, .9)) figs[-1].tight_layout(rect=(0, 0, 1, .86))
# ax3 = pylab.subplot2grid(splotshape, (3, 0), 2, 4, sharex=ax2) # ax3 = pylab.subplot2grid(splotshape, (3, 0), 2, 4, sharex=ax2)
figs.append(pylab.figure("BGPLVM DEBUG S", figsize=(12, 4))) figs.append(pylab.figure("BGPLVM DEBUG S", figsize=(12, 4)))
ax3 = self._debug_get_axis(figs) ax3 = self._debug_get_axis(figs)
ax3.text(.5, .5, r"$\mathbf{S}$", alpha=.5, transform=ax3.transAxes, ax3.text(.5, .5, r"$\mathbf{S}$", alpha=.5, transform=ax3.transAxes,
ha='center', va='center') ha='center', va='center')
figs[-1].canvas.draw() figs[-1].canvas.draw()
figs[-1].tight_layout(rect=(0, 0, 1, .9)) figs[-1].tight_layout(rect=(0, 0, 1, .86))
# ax4 = pylab.subplot2grid(splotshape, (5, 0), 2, 2) # ax4 = pylab.subplot2grid(splotshape, (5, 0), 2, 2)
figs.append(pylab.figure("BGPLVM DEBUG Z", figsize=(6, 4))) figs.append(pylab.figure("BGPLVM DEBUG Z", figsize=(6, 4)))
ax4 = self._debug_get_axis(figs) ax4 = self._debug_get_axis(figs)
ax4.text(.5, .5, r"$\mathbf{Z}$", alpha=.5, transform=ax4.transAxes, ax4.text(.5, .5, r"$\mathbf{Z}$", alpha=.5, transform=ax4.transAxes,
ha='center', va='center') ha='center', va='center')
figs[-1].canvas.draw() figs[-1].canvas.draw()
figs[-1].tight_layout(rect=(0, 0, 1, .9)) figs[-1].tight_layout(rect=(0, 0, 1, .86))
# ax5 = pylab.subplot2grid(splotshape, (5, 2), 2, 2) # ax5 = pylab.subplot2grid(splotshape, (5, 2), 2, 2)
figs.append(pylab.figure("BGPLVM DEBUG theta", figsize=(6, 4))) figs.append(pylab.figure("BGPLVM DEBUG theta", figsize=(6, 4)))
ax5 = self._debug_get_axis(figs) ax5 = self._debug_get_axis(figs)
ax5.text(.5, .5, r"${\theta}$", alpha=.5, transform=ax5.transAxes, ax5.text(.5, .5, r"${\theta}$", alpha=.5, transform=ax5.transAxes,
ha='center', va='center') ha='center', va='center')
figs[-1].canvas.draw() figs[-1].canvas.draw()
figs[-1].tight_layout(rect=(0, 0, 1, .9)) figs[-1].tight_layout(rect=(.15, 0, 1, .86))
figs.append(pylab.figure("BGPLVM DEBUG Kmm", figsize=(12, 6))) figs.append(pylab.figure("BGPLVM DEBUG Kmm", figsize=(12, 6)))
fig = figs[-1] fig = figs[-1]
ax6 = fig.add_subplot(121) ax6 = fig.add_subplot(121)
@ -308,6 +309,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
Slatentgrads = ax3.quiver(xlatent, S, Ulatent, Sg, color=colors, Slatentgrads = ax3.quiver(xlatent, S, Ulatent, Sg, color=colors,
units=quiver_units, scale_units=quiver_scale_units, units=quiver_units, scale_units=quiver_scale_units,
scale=quiver_scale) scale=quiver_scale)
ax3.set_ylim(0, 1.)
xZ = np.tile(np.arange(0, Z.shape[0])[:, None], Z.shape[1]) xZ = np.tile(np.arange(0, Z.shape[0])[:, None], Z.shape[1])
UZ = np.zeros_like(Z) UZ = np.zeros_like(Z)
@ -343,16 +345,16 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
# loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.15, 1, 1.15), # loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.15, 1, 1.15),
# borderaxespad=0, mode="expand") # borderaxespad=0, mode="expand")
ax2.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], ax2.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
borderaxespad=0, mode="expand") borderaxespad=0, mode="expand")
ax3.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], ax3.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
borderaxespad=0, mode="expand") borderaxespad=0, mode="expand")
ax4.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], ax4.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
borderaxespad=0, mode="expand") borderaxespad=0, mode="expand")
ax5.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)], ax5.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.01, 1, 1.01), loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
borderaxespad=0, mode="expand") borderaxespad=0, mode="expand")
Lleg = ax1.legend() Lleg = ax1.legend()
Lleg.draggable() Lleg.draggable()
@ -427,11 +429,11 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
cbarkmmdl.update_normal(imkmmdl) cbarkmmdl.update_normal(imkmmdl)
ax2.relim() ax2.relim()
ax3.relim() # ax3.relim()
ax4.relim() ax4.relim()
ax5.relim() ax5.relim()
ax2.autoscale() ax2.autoscale()
ax3.autoscale() # ax3.autoscale()
ax4.autoscale() ax4.autoscale()
ax5.autoscale() ax5.autoscale()

13
GPy/models/GP.py
View file

@ -35,6 +35,9 @@ class GP(model):
self.N, self.Q = self.X.shape self.N, self.Q = self.X.shape
assert isinstance(kernel, kern.kern) assert isinstance(kernel, kern.kern)
self.kern = kernel self.kern = kernel
self.likelihood = likelihood
assert self.X.shape[0] == self.likelihood.data.shape[0]
self.N, self.D = self.likelihood.data.shape
# here's some simple normalization for the inputs # here's some simple normalization for the inputs
if normalize_X: if normalize_X:
@ -47,12 +50,8 @@ class GP(model):
self._Xmean = np.zeros((1, self.X.shape[1])) self._Xmean = np.zeros((1, self.X.shape[1]))
self._Xstd = np.ones((1, self.X.shape[1])) self._Xstd = np.ones((1, self.X.shape[1]))
self.likelihood = likelihood if not hasattr(self,'has_uncertain_inputs'):
# assert self.X.shape[0] == self.likelihood.Y.shape[0] self.has_uncertain_inputs = False
# self.N, self.D = self.likelihood.Y.shape
assert self.X.shape[0] == self.likelihood.data.shape[0]
self.N, self.D = self.likelihood.data.shape
model.__init__(self) model.__init__(self)
def dL_dZ(self): def dL_dZ(self):
@ -232,7 +231,7 @@ class GP(model):
else: else:
raise NotImplementedError, "Cannot define a frame with more than two input dimensions" raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
def plot(self, samples=0, plot_limits=None, which_data='all', which_functions='all', resolution=None, levels=20): def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20):
""" """
TODO: Docstrings! TODO: Docstrings!
:param levels: for 2D plotting, the number of contour levels to use :param levels: for 2D plotting, the number of contour levels to use

93
GPy/models/sparse_GP.py
View file

@ -9,6 +9,11 @@ from .. import kern
from GP import GP from GP import GP
from scipy import linalg from scipy import linalg
def backsub_both_sides(L,X):
""" Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky"""
tmp,_ = linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(X),lower=1,trans=1)
return linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(tmp.T),lower=1,trans=1)[0].T
class sparse_GP(GP): class sparse_GP(GP):
""" """
Variational sparse GP model Variational sparse GP model
@ -68,51 +73,64 @@ class sparse_GP(GP):
sf = self.scale_factor sf = self.scale_factor
sf2 = sf**2 sf2 = sf**2
#The rather complex computations of psi2_beta_scaled #invert Kmm
self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
#The rather complex computations of self.A
if self.likelihood.is_heteroscedastic: if self.likelihood.is_heteroscedastic:
assert self.likelihood.D == 1 #TODO: what if the likelihood is heterscedatic and there are multiple independent outputs? assert self.likelihood.D == 1 #TODO: what if the likelihood is heterscedatic and there are multiple independent outputs?
if self.has_uncertain_inputs: if self.has_uncertain_inputs:
self.psi2_beta_scaled = (self.psi2*(self.likelihood.precision.flatten().reshape(self.N,1,1)/sf2)).sum(0) psi2_beta_scaled = (self.psi2*(self.likelihood.precision.flatten().reshape(self.N,1,1)/sf2)).sum(0)
evals, evecs = linalg.eigh(psi2_beta_scaled)
clipped_evals = np.clip(evals,0.,1e6) # TODO: make clipping configurable
if not np.allclose(evals, clipped_evals):
print "Warning: clipping posterior eigenvalues"
tmp = evecs*np.sqrt(clipped_evals)
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1)
self.A = tdot(tmp)
else: else:
tmp = self.psi1*(np.sqrt(self.likelihood.precision.flatten().reshape(1,self.N))/sf) tmp = self.psi1*(np.sqrt(self.likelihood.precision.flatten().reshape(1,self.N))/sf)
#self.psi2_beta_scaled = np.dot(tmp,tmp.T) #self.psi2_beta_scaled = tdot(tmp)
self.psi2_beta_scaled = tdot(tmp) tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1)
self.A = tdot(tmp)
else: else:
if self.has_uncertain_inputs: if self.has_uncertain_inputs:
self.psi2_beta_scaled = (self.psi2*(self.likelihood.precision/sf2)).sum(0) psi2_beta_scaled = (self.psi2*(self.likelihood.precision/sf2)).sum(0)
evals, evecs = linalg.eigh(psi2_beta_scaled)
clipped_evals = np.clip(evals,0.,1e6) # TODO: make clipping configurable
if not np.allclose(evals, clipped_evals):
print "Warning: clipping posterior eigenvalues"
tmp = evecs*np.sqrt(clipped_evals)
#self.psi2_beta_scaled = tdot(tmp)
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1)
self.A = tdot(tmp)
else: else:
tmp = self.psi1*(np.sqrt(self.likelihood.precision)/sf) tmp = self.psi1*(np.sqrt(self.likelihood.precision)/sf)
#self.psi2_beta_scaled = np.dot(tmp,tmp.T) #self.psi2_beta_scaled = tdot(tmp)
self.psi2_beta_scaled = tdot(tmp) tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp),lower=1)
self.A = tdot(tmp)
self.Kmmi, self.Lm, self.Lmi, self.Kmm_logdet = pdinv(self.Kmm)
self.V = (self.likelihood.precision/self.scale_factor)*self.likelihood.Y
#Compute A = L^-1 psi2 beta L^-T
#self. A = mdot(self.Lmi,self.psi2_beta_scaled,self.Lmi.T)
tmp = linalg.lapack.flapack.dtrtrs(self.Lm,self.psi2_beta_scaled.T,lower=1)[0]
self.A = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1)[0]
#invert B and compute C. C is the posterior covariance of u
self.B = np.eye(self.M)/sf2 + self.A self.B = np.eye(self.M)/sf2 + self.A
self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B) self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B)
self.psi1V = np.dot(self.psi1, self.V)
tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.Bi),lower=1,trans=1)[0] tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.Bi),lower=1,trans=1)[0]
self.C = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0] self.C = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0]
self.V = (self.likelihood.precision/self.scale_factor)*self.likelihood.Y
self.psi1V = np.dot(self.psi1, self.V)
#back substutue C into psi1V #back substutue C into psi1V
tmp,info1 = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.psi1V),lower=1,trans=0) tmp,info1 = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.psi1V),lower=1,trans=0)
self._LBi_Lmi_psi1V,_ = linalg.lapack.flapack.dtrtrs(self.LB,np.asfortranarray(tmp),lower=1,trans=0)
self._P = tdot(tmp)
tmp,info2 = linalg.lapack.flapack.dpotrs(self.LB,tmp,lower=1) tmp,info2 = linalg.lapack.flapack.dpotrs(self.LB,tmp,lower=1)
self.Cpsi1V,info3 = linalg.lapack.flapack.dtrtrs(self.Lm,tmp,lower=1,trans=1) self.Cpsi1V,info3 = linalg.lapack.flapack.dtrtrs(self.Lm,tmp,lower=1,trans=1)
#self.Cpsi1V = np.dot(self.C,self.psi1V) #self.Cpsi1V = np.dot(self.C,self.psi1V)
self.Cpsi1VVpsi1 = np.dot(self.Cpsi1V,self.psi1V.T)
self.E = tdot(self.Cpsi1V/sf) self.E = tdot(self.Cpsi1V/sf)
# Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertin inputs case # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertin inputs case
self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten() self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten()
self.dL_dpsi1 = np.dot(self.Cpsi1V,self.V.T) self.dL_dpsi1 = np.dot(self.Cpsi1V,self.V.T)
@ -130,27 +148,23 @@ class sparse_GP(GP):
self.dL_dpsi2 = None self.dL_dpsi2 = None
else: else:
#self.dL_dpsi2 = 0.5 * self.likelihood.precision * self.D * self.Kmmi # dB
#self.dL_dpsi2 += - 0.5 * self.likelihood.precision/sf2 * self.D * self.C # dC
#self.dL_dpsi2 += - 0.5 * self.likelihood.precision * self.E # dD
self.dL_dpsi2 = 0.5*self.likelihood.precision*(self.D*(self.Kmmi - self.C/sf2) -self.E) self.dL_dpsi2 = 0.5*self.likelihood.precision*(self.D*(self.Kmmi - self.C/sf2) -self.E)
if self.has_uncertain_inputs: if self.has_uncertain_inputs:
#repeat for each of the N psi_2 matrices #repeat for each of the N psi_2 matrices
self.dL_dpsi2 = np.repeat(self.dL_dpsi2[None,:,:],self.N,axis=0) self.dL_dpsi2 = np.repeat(self.dL_dpsi2[None,:,:],self.N,axis=0)
else: else:
#subsume back into psi1 (==Kmn)
self.dL_dpsi1 += 2.*np.dot(self.dL_dpsi2,self.psi1) self.dL_dpsi1 += 2.*np.dot(self.dL_dpsi2,self.psi1)
self.dL_dpsi2 = None self.dL_dpsi2 = None
# Compute dL_dKmm # Compute dL_dKmm
#self.dL_dKmm_old = -0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB tmp = tdot(self._LBi_Lmi_psi1V)
#self.dL_dKmm += -0.5 * self.D * (- self.C/sf2 - 2.*mdot(self.C, self.psi2_beta_scaled, self.Kmmi) + self.Kmmi) # dC self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.D*np.eye(self.M) + tmp)
#self.dL_dKmm += np.dot(np.dot(self.E*sf2, self.psi2_beta_scaled) - self.Cpsi1VVpsi1, self.Kmmi) + 0.5*self.E # dD tmp = -0.5*self.DBi_plus_BiPBi/sf2
tmp = linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(self.B),lower=1,trans=1)[0] tmp += -0.5*self.B*sf2*self.D
self.dL_dKmm = -0.5*self.D*sf2*linalg.lapack.flapack.dtrtrs(self.Lm,np.asfortranarray(tmp.T),lower=1,trans=1)[0] #dA tmp += self.D*np.eye(self.M)
tmp = np.dot(self.D*self.C + self.E*sf2,self.psi2_beta_scaled) - self.Cpsi1VVpsi1 self.dL_dKmm = backsub_both_sides(self.Lm,tmp)
tmp = linalg.lapack.flapack.dpotrs(self.Lm,np.asfortranarray(tmp.T),lower=1)[0].T
self.dL_dKmm += 0.5*(self.D*self.C/sf2 + self.E) +tmp # d(C+D)
#the partial derivative vector for the likelihood #the partial derivative vector for the likelihood
if self.likelihood.Nparams ==0: if self.likelihood.Nparams ==0:
@ -166,8 +180,9 @@ class sparse_GP(GP):
#likelihood is not heterscedatic #likelihood is not heterscedatic
self.partial_for_likelihood = - 0.5 * self.N*self.D*self.likelihood.precision + 0.5 * self.likelihood.trYYT*self.likelihood.precision**2 self.partial_for_likelihood = - 0.5 * self.N*self.D*self.likelihood.precision + 0.5 * self.likelihood.trYYT*self.likelihood.precision**2
self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum()*self.likelihood.precision**2 - np.trace(self.A)*self.likelihood.precision*sf2) self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum()*self.likelihood.precision**2 - np.trace(self.A)*self.likelihood.precision*sf2)
self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self.likelihood.precision #self.partial_for_likelihood += 0.5 * self.D * trace_dot(self.Bi,self.A)*self.likelihood.precision
self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1)) #self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.sum(np.square(self._LBi_Lmi_psi1V)))
self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.A,self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V)))
@ -181,7 +196,7 @@ class sparse_GP(GP):
A = -0.5*self.N*self.D*(np.log(2.*np.pi) + np.log(self.likelihood._variance)) -0.5*self.likelihood.precision*self.likelihood.trYYT A = -0.5*self.N*self.D*(np.log(2.*np.pi) + np.log(self.likelihood._variance)) -0.5*self.likelihood.precision*self.likelihood.trYYT
B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2) B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2)
C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2)) C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2))
D = 0.5*np.trace(self.Cpsi1VVpsi1) D = 0.5*np.sum(np.square(self._LBi_Lmi_psi1V))
return A+B+C+D return A+B+C+D
def _set_params(self, p): def _set_params(self, p):
@ -189,14 +204,14 @@ class sparse_GP(GP):
self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam]) self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam])
self.likelihood._set_params(p[self.Z.size+self.kern.Nparam:]) self.likelihood._set_params(p[self.Z.size+self.kern.Nparam:])
self._compute_kernel_matrices() self._compute_kernel_matrices()
if self.auto_scale_factor: #if self.auto_scale_factor:
self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision) # self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision)
#if self.auto_scale_factor: #if self.auto_scale_factor:
#if self.likelihood.is_heteroscedastic: #if self.likelihood.is_heteroscedastic:
# self.scale_factor = max(1,np.sqrt(self.psi2_beta_scaled.sum(0).mean())) #self.scale_factor = max(100,np.sqrt(self.psi2_beta_scaled.sum(0).mean()))
#else: #else:
#self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision) #self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision)
#self.scale_factor = 1. self.scale_factor = 1.
self._computations() self._computations()
def _get_params(self): def _get_params(self):

12
GPy/testing/__init__.py
View file

@ -0,0 +1,12 @@
"""
MaxZ
"""
import unittest
import sys
def deepTest(reason):
if 'deep' in sys.argv:
return lambda x:x
return unittest.skip("Not deep scanning, enable deepscan by adding 'deep' argument")

112
GPy/testing/cgd_tests.py Normal file
View file

@ -0,0 +1,112 @@
'''
Created on 26 Apr 2013
@author: maxz
'''
import unittest
import numpy
from GPy.inference.conjugate_gradient_descent import CGD, RUNNING
import pylab
import time
from scipy.optimize.optimize import rosen, rosen_der
class Test(unittest.TestCase):
def testMinimizeSquare(self):
N = 100
A = numpy.random.rand(N) * numpy.eye(N)
b = numpy.random.rand(N) * 0
f = lambda x: numpy.dot(x.T.dot(A), x) - numpy.dot(x.T, b)
df = lambda x: numpy.dot(A, x) - b
opt = CGD()
restarts = 10
for _ in range(restarts):
try:
x0 = numpy.random.randn(N) * 300
res = opt.opt(f, df, x0, messages=0,
maxiter=1000, gtol=1e-10)
assert numpy.allclose(res[0], 0, atol=1e-3)
break
except:
# RESTART
pass
else:
raise AssertionError("Test failed for {} restarts".format(restarts))
def testRosen(self):
N = 20
f = rosen
df = rosen_der
opt = CGD()
restarts = 10
for _ in range(restarts):
try:
x0 = numpy.random.randn(N) * .5
res = opt.opt(f, df, x0, messages=0,
maxiter=5e2, gtol=1e-2)
assert numpy.allclose(res[0], 1, atol=.1)
break
except:
# RESTART
pass
else:
raise AssertionError("Test failed for {} restarts".format(restarts))
if __name__ == "__main__":
# import sys;sys.argv = ['',
# 'Test.testMinimizeSquare',
# 'Test.testRosen',
# ]
# unittest.main()
N = 2
A = numpy.random.rand(N) * numpy.eye(N)
b = numpy.random.rand(N) * 0
# f = lambda x: numpy.dot(x.T.dot(A), x) - numpy.dot(x.T, b)
# df = lambda x: numpy.dot(A, x) - b
f = rosen
df = rosen_der
x0 = numpy.random.randn(N) * .5
opt = CGD()
fig = pylab.figure("cgd optimize")
if fig.axes:
ax = fig.axes[0]
ax.cla()
else:
ax = fig.add_subplot(111, projection='3d')
interpolation = 40
x, y = numpy.linspace(-1, 1, interpolation)[:, None], numpy.linspace(-1, 1, interpolation)[:, None]
X, Y = numpy.meshgrid(x, y)
fXY = numpy.array([f(numpy.array([x, y])) for x, y in zip(X.flatten(), Y.flatten())]).reshape(interpolation, interpolation)
ax.plot_wireframe(X, Y, fXY)
xopts = [x0.copy()]
optplts, = ax.plot3D([x0[0]], [x0[1]], zs=f(x0), marker='o', color='r')
raw_input("enter to start optimize")
res = [0]
def callback(*r):
xopts.append(r[0].copy())
# time.sleep(.3)
optplts._verts3d = [numpy.array(xopts)[:, 0], numpy.array(xopts)[:, 1], [f(xs) for xs in xopts]]
fig.canvas.draw()
if r[-1] != RUNNING:
res[0] = r
p, c = opt.opt_async(f, df, x0.copy(), callback, messages=True, maxiter=1000,
report_every=20, gtol=1e-12)
pylab.ion()
pylab.show()
pass

57
GPy/testing/kern_psi_stat_tests.py
View file

@ -6,19 +6,31 @@ Created on 26 Apr 2013
import unittest import unittest
import GPy import GPy
import numpy as np import numpy as np
import pylab import sys
from .. import testing
__test__ = False __test__ = True
np.random.seed(0)
def ard(p):
try:
if p.ARD:
return "ARD"
except:
pass
return ""
@testing.deepTest
class Test(unittest.TestCase): class Test(unittest.TestCase):
D = 9 D = 9
M = 5 M = 4
Nsamples = 3e6 N = 3
Nsamples = 6e6
def setUp(self): def setUp(self):
self.kerns = ( self.kerns = (
GPy.kern.rbf(self.D), GPy.kern.rbf(self.D, ARD=True), GPy.kern.rbf(self.D), GPy.kern.rbf(self.D, ARD=True),
GPy.kern.linear(self.D), GPy.kern.linear(self.D, ARD=True), GPy.kern.linear(self.D, ARD=False), GPy.kern.linear(self.D, ARD=True),
GPy.kern.linear(self.D) + GPy.kern.bias(self.D), GPy.kern.linear(self.D) + GPy.kern.bias(self.D),
GPy.kern.rbf(self.D) + GPy.kern.bias(self.D), GPy.kern.rbf(self.D) + GPy.kern.bias(self.D),
GPy.kern.linear(self.D) + GPy.kern.bias(self.D) + GPy.kern.white(self.D), GPy.kern.linear(self.D) + GPy.kern.bias(self.D) + GPy.kern.white(self.D),
@ -43,16 +55,26 @@ class Test(unittest.TestCase):
for kern in self.kerns: for kern in self.kerns:
Nsamples = 100 Nsamples = 100
psi1 = kern.psi1(self.Z, self.q_x_mean, self.q_x_variance) psi1 = kern.psi1(self.Z, self.q_x_mean, self.q_x_variance)
K_ = np.zeros((self.N, self.M)) K_ = np.zeros((Nsamples, self.M))
diffs = [] diffs = []
for i, q_x_sample_stripe in enumerate(np.array_split(self.q_x_samples, self.Nsamples / Nsamples)): for i, q_x_sample_stripe in enumerate(np.array_split(self.q_x_samples, self.Nsamples / Nsamples)):
K = kern.K(q_x_sample_stripe, self.Z) K = kern.K(q_x_sample_stripe, self.Z)
K_ += K K_ += K
diffs.append(((psi1 - (K_ / (i + 1))) ** 2).mean()) diffs.append(((psi1 - (K_ / (i + 1)))).mean())
K_ /= self.Nsamples / Nsamples K_ /= self.Nsamples / Nsamples
# pylab.figure("+".join([p.name for p in kern.parts]) + "psi1") msg = "psi1: " + "+".join([p.name + ard(p) for p in kern.parts])
try:
# pylab.figure(msg)
# pylab.plot(diffs) # pylab.plot(diffs)
self.assertTrue(np.allclose(psi1.flatten() , K.mean(0), rtol=1e-1)) self.assertTrue(np.allclose(psi1.squeeze(), K_,
rtol=1e-1, atol=.1),
msg=msg + ": not matching")
# sys.stdout.write(".")
except:
# import ipdb;ipdb.set_trace()
# kern.psi2(self.Z, self.q_x_mean, self.q_x_variance)
# sys.stdout.write("E") # msg + ": not matching"
pass
def test_psi2(self): def test_psi2(self):
for kern in self.kerns: for kern in self.kerns:
@ -64,20 +86,27 @@ class Test(unittest.TestCase):
K = kern.K(q_x_sample_stripe, self.Z) K = kern.K(q_x_sample_stripe, self.Z)
K = (K[:, :, None] * K[:, None, :]).mean(0) K = (K[:, :, None] * K[:, None, :]).mean(0)
K_ += K K_ += K
diffs.append(((psi2 - (K_ / (i + 1))) ** 2).mean()) diffs.append(((psi2 - (K_ / (i + 1)))).mean())
K_ /= self.Nsamples / Nsamples K_ /= self.Nsamples / Nsamples
msg = "psi2: {}".format("+".join([p.name + ard(p) for p in kern.parts]))
try: try:
# pylab.figure("+".join([p.name for p in kern.parts]) + "psi2") # pylab.figure(msg)
# pylab.plot(diffs) # pylab.plot(diffs)
self.assertTrue(np.allclose(psi2.squeeze(), K_, self.assertTrue(np.allclose(psi2.squeeze(), K_,
rtol=1e-1, atol=.1), rtol=1e-1, atol=.1),
msg="{}: not matching".format("+".join([p.name for p in kern.parts]))) msg=msg + ": not matching")
# sys.stdout.write(".")
except: except:
print "{}: not matching".format(kern.parts[0].name) # import ipdb;ipdb.set_trace()
# kern.psi2(self.Z, self.q_x_mean, self.q_x_variance)
# sys.stdout.write("E")
print msg + ": not matching"
pass
if __name__ == "__main__": if __name__ == "__main__":
import sys;sys.argv = ['', import sys;sys.argv = ['',
'Test.test_psi0', 'Test.test_psi0',
'Test.test_psi1', 'Test.test_psi1',
'Test.test_psi2'] 'Test.test_psi2',
]
unittest.main() unittest.main()

60
GPy/testing/psi_stat_tests.py
View file

@ -6,7 +6,6 @@ Created on 22 Apr 2013
import unittest import unittest
import numpy import numpy
from GPy.models.Bayesian_GPLVM import Bayesian_GPLVM
import GPy import GPy
import itertools import itertools
from GPy.core import model from GPy.core import model
@ -48,16 +47,16 @@ class PsiStatModel(model):
thetagrad = self.kern.__getattribute__("d" + self.which + "_dtheta")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance).flatten() thetagrad = self.kern.__getattribute__("d" + self.which + "_dtheta")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance).flatten()
return numpy.hstack((psimu.flatten(), psiS.flatten(), psiZ.flatten(), thetagrad)) return numpy.hstack((psimu.flatten(), psiS.flatten(), psiZ.flatten(), thetagrad))
class Test(unittest.TestCase): class DPsiStatTest(unittest.TestCase):
Q = 5 Q = 5
N = 50 N = 50
M = 10 M = 10
D = 10 D = 20
X = numpy.random.randn(N, Q) X = numpy.random.randn(N, Q)
X_var = .5 * numpy.ones_like(X) + .4 * numpy.clip(numpy.random.randn(*X.shape), 0, 1) X_var = .5 * numpy.ones_like(X) + .4 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
Z = numpy.random.permutation(X)[:M] Z = numpy.random.permutation(X)[:M]
Y = X.dot(numpy.random.randn(Q, D)) Y = X.dot(numpy.random.randn(Q, D))
kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q)] # kernels = [GPy.kern.linear(Q, ARD=True, variances=numpy.random.rand(Q)), GPy.kern.rbf(Q, ARD=True), GPy.kern.bias(Q)]
kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q), kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
GPy.kern.linear(Q) + GPy.kern.bias(Q), GPy.kern.linear(Q) + GPy.kern.bias(Q),
@ -67,7 +66,10 @@ class Test(unittest.TestCase):
for k in self.kernels: for k in self.kernels:
m = PsiStatModel('psi0', X=self.X, X_variance=self.X_var, Z=self.Z, m = PsiStatModel('psi0', X=self.X, X_variance=self.X_var, Z=self.Z,
M=self.M, kernel=k) M=self.M, kernel=k)
try:
assert m.checkgrad(), "{} x psi0".format("+".join(map(lambda x: x.name, k.parts))) assert m.checkgrad(), "{} x psi0".format("+".join(map(lambda x: x.name, k.parts)))
except:
import ipdb;ipdb.set_trace()
# def testPsi1(self): # def testPsi1(self):
# for k in self.kernels: # for k in self.kernels:
@ -106,31 +108,31 @@ if __name__ == "__main__":
import sys import sys
interactive = 'i' in sys.argv interactive = 'i' in sys.argv
if interactive: if interactive:
N, M, Q, D = 30, 5, 4, 30 # N, M, Q, D = 30, 5, 4, 30
X = numpy.random.rand(N, Q) # X = numpy.random.rand(N, Q)
k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) # k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
K = k.K(X) # K = k.K(X)
Y = numpy.random.multivariate_normal(numpy.zeros(N), K, D).T # Y = numpy.random.multivariate_normal(numpy.zeros(N), K, D).T
Y -= Y.mean(axis=0) # Y -= Y.mean(axis=0)
k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001) # k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M) # m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
m.ensure_default_constraints() # m.ensure_default_constraints()
m.randomize() # m.randomize()
# self.assertTrue(m.checkgrad()) # # self.assertTrue(m.checkgrad())
numpy.random.seed(0)
Q = 5 Q = 5
N = 50 N = 50
M = 10 M = 10
D = 10 D = 15
X = numpy.random.randn(N, Q) X = numpy.random.randn(N, Q)
X_var = .5 * numpy.ones_like(X) + .4 * numpy.clip(numpy.random.randn(*X.shape), 0, 1) X_var = .5 * numpy.ones_like(X) + .1 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
Z = numpy.random.permutation(X)[:M] Z = numpy.random.permutation(X)[:M]
Y = X.dot(numpy.random.randn(Q, D)) Y = X.dot(numpy.random.randn(Q, D))
kernel = GPy.kern.bias(Q) # kernel = GPy.kern.bias(Q)
#
kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q), # kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
GPy.kern.linear(Q) + GPy.kern.bias(Q), # GPy.kern.linear(Q) + GPy.kern.bias(Q),
GPy.kern.rbf(Q) + GPy.kern.bias(Q)] # GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
# for k in kernels: # for k in kernels:
# m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z, # m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
@ -143,11 +145,13 @@ if __name__ == "__main__":
# M=M, kernel=kernel) # M=M, kernel=kernel)
# m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z, # m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
# M=M, kernel=kernel) # M=M, kernel=kernel)
m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z, # m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
M=M, kernel=GPy.kern.rbf(Q)) # M=M, kernel=GPy.kern.rbf(Q))
m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z, m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
M=M, kernel=GPy.kern.linear(Q) + GPy.kern.bias(Q)) M=M, kernel=GPy.kern.linear(Q, ARD=True, variances=numpy.random.rand(Q)))
m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z, m3.ensure_default_constraints()
M=M, kernel=GPy.kern.rbf(Q) + GPy.kern.bias(Q)) # + GPy.kern.bias(Q))
# m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
# M=M, kernel=GPy.kern.rbf(Q) + GPy.kern.bias(Q))
else: else:
unittest.main() unittest.main()

BIN
GPy/util/datasets/BGPLVMSimulation.mat Normal file
View file

Binary file not shown.

29
GPy/util/linalg.py
View file

@ -72,7 +72,7 @@ def jitchol(A,maxtries=5):
raise linalg.LinAlgError, "not pd: negative diagonal elements" raise linalg.LinAlgError, "not pd: negative diagonal elements"
jitter= diagA.mean()*1e-6 jitter= diagA.mean()*1e-6
for i in range(1,maxtries+1): for i in range(1,maxtries+1):
print '\rWarning: adding jitter of {:.10e} '.format(jitter), print 'Warning: adding jitter of {:.10e}'.format(jitter)
try: try:
return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True) return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True)
except: except:
@ -276,3 +276,30 @@ def symmetrify_murray(A):
nn = A.shape[0] nn = A.shape[0]
A[[range(nn),range(nn)]] /= 2.0 A[[range(nn),range(nn)]] /= 2.0
def cholupdate(L,x):
"""
update the LOWER cholesky factor of a pd matrix IN PLACE
if L is the lower chol. of K, then this function computes L_
where L_ is the lower chol of K + x*x^T
"""
support_code = """
#include <math.h>
"""
code="""
double r,c,s;
int j,i;
for(j=0; j<N; j++){
r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
c = r / L(j,j);
s = x(j) / L(j,j);
L(j,j) = r;
for (i=j+1; i<N; i++){
L(i,j) = (L(i,j) + s*x(i))/c;
x(i) = c*x(i) - s*L(i,j);
}
}
"""
x = x.copy()
N = x.size
weave.inline(code, support_code=support_code, arg_names=['N','L','x'], type_converters=weave.converters.blitz)