gpflow / gpflow Goto Github PK

The following code gives an error.

import numpy as np
import GPflow

rng = np.random.RandomState(1)
X = rng.randn( 10, 1)
Y = rng.randn( 10, 1 ) 
Z = rng.randn( 3,1 )
model = GPflow.svgp.SVGP( X=X, Y=Y, kern=GPflow.kernels.RBF(1) ,likelihood=GPflow.likelihoods.Gaussian() , Z=Z )
model.compute_log_likelihood()

FailedPreconditionError: Attempting to use uninitialized value tfY

It seems this may be down to recent changes?

Some GPflow users are interested in interdomain inducing points. At the moment it is a little fiddly to integrate their code with the main code base. We should think about making it easier without disrupting the main code base.

Here:

https://github.com/GPflow/GPflow/blob/master/testing/test_method_equivalence.py#L26

Some of the methods are given an additional white noise term in the kernel and some are not. Why?

There have been several improvements to TensorFlow which mean that we need fewer UserOps. Where possible we should switch over to using the official versions of Ops.

I get am import error when trying to import GPflow:
"ImportError: numpy.core.multiarray failed to import"
I have followed the pip installation advice (on a Ubuntu machine with Python 2.7.6) and downgraded numpy to version 1.9.3.

Hi team,

I am considering implementing in GPflow the following paper as a start to multi-output GPs:
https://papers.nips.cc/paper/3189-multi-task-gaussian-process-prediction.pdf

It's a kronecker product kernel where the between-output correlation is estimated using maximum likelihood.
I had implemented this a while back in Matlab and with autodifferentiation it should be a breeze to implement in GPflow.

Do you think this is worth doing or there are much better approaches out there we should implement instead? I like the kronecker approach because it's so simple and interpretable.

A couple of incorrect behaviours of code coverage that I have observed:

1. The code coverage shown on the front page of the repo is currently dispayed as the coverage of the last build even if that build is not master. So for example the coverage on the front page reads 95% when in fact master has 98% coverage because the last build was on the GPLVM pull request. In this case the code coverage is misleadingly low for new visitors to our front page.

2. A similar problem. The merge criterion that the coverage on a branch should not be less than master is currently that the coverage should not be less than the last autobuild. Thus if you build a branch with 50% test coverage and the next build has coverage of 75% then it will pass even if master has coverage of 98% which is very wrong behaviour that could damage the coverage we have worked so hard to obtain over time.

What would the best way be to combine mean functions? Ideally I would like to add addition and multiplication to the MeanFunction class, but I keep running into problems with the Parameter specifications (and I feel that your setup would allow a more elegant solution that implementing multiplicative/additive mean function classes).

It would be good to have the multi class likelihood sucessfully used in this paper

https://papers.nips.cc/paper/5875-mcmc-for-variationally-sparse-gaussian-processes.pdf

in GPflow . Ideally this would be done using fast ops written in TensorFlow .

Should the kernel evaluation functions also be decorated with Autoflow? It would allow for easy evaluation of kernels, which is great for illustration and debugging.

I did some early profiling with GPflow to try and understand the computational bottlenecks but I was somewhat held back by the availability of tools.

Profiling support has however improved. See this TensorFlow issue:

tensorflow/tensorflow#1824

Most of the documentation is currently in the discussion of the issue I raised. I haven't yet had time to try it.

We have done some work on getting minibatches working in GPflow. We should consolidate it into its simplest and cleanest form and then submit a pull request onto master.

We need to update the list of contributors in alphabetical order on the readme.

Although I prefer the Titsias method for a number of reasons, many people still want to compare to FITC which was first created by Snelson and Ghahramani 2006 Sparse Gaussian Processes using Pseudo-inputs. Therefore I am going to implement it in GPflow.

Sometimes it's necessary to reseat tensorflow's graph for memory reasons. This can be done with tf.reset_default_graph().

Problem is, this breaks all out AutoFlow functions! Here's a MWE:

tf.reset_default_graph()
k = GPflow.kernels.Matern32(1)
X, Y = np.random.randn(2, 10, 1)
m = GPflow.gpr.GPR(X, Y, kern=k)
m.predict_f(m.X)

This is because the Autoflow constructor takes placeholders as arguments. After reset_default_graph, those placeholders are not on the same graph as everything else. I'm going to propose a fix where we just give AutoFlow enough information to make the placeholders, not the placeholders themselves.

It would be good to have the experiment with FITC documented in an IPython notebook.

To some extent this is a continuation of #37

Opening this here to get thoughts on the best design for integrating Blitzkriging into GPflow - I have the tensorflow ops in place to get it all working but would appreciate the maintainers thoughts on how to structure the code.

In the short term I think I will implement the functionality I need based after the functions in svgp.py - for example, I'd make a set of new KL divergence functions with the calling convention of kullback_leiblers.gauss_kl(...).

Do you think it would be better to integrate the kronecker approximation into SVGP directly, or as a separate class? I'm tending towards a separate class to avoid multiple levels of if statements but I'm receptive to other ideas.

Am I correct in thinking that if I make a variable a Param it will get optimised automatically?

At the moment in this line:

https://github.com/GPflow/GPflow/blob/master/GPflow/model.py#L291

looks like the tolerance was set tentatively. What would be the best way to pass through the tolerance to the various scipy optimizers?

I need to be able to manipulate the stopping tolerance for some experiments I have been doing in #38 and #37 .

Do we have in mind a particular way to serialize and de-serialize models so that they can be stored? Pickle doesn't currently work.

Parameterized objects within a ParamList do not get the tf_mode set correctly because the parameter class incorrectly uses this list comprehension:

[child._begin_tf_mode() for key, child in self.__dict__.items() if isinstance(child, Parameterized) and key is not '_parent']

Which duplicates the definition of self.sorted_params:

    @property
    def sorted_params(self):
        #snip
        params = [child for key, child in self.__dict__.items() if isinstance(child, (Param, Parameterized)) and key is not '_parent']
        return sorted(params, key=id)

Do we want to maybe include the ability to get the diagonal terms of the predictive covariance? Clearly this can be intractable for large vectors but there may be other cases where it is useful.

It would be useful to return predictive covariance as well as predictive variance. Clearly this is intractable for large predictive sets but still worth thinking about in my opinion for smaller ones.

Your branch of tensorflow support GPU？

Greetings,

First of all, nice work on this library! =)
I tried to replace some GPy models with GPflow ones in some stuff I am doing but I noticed I was getting worse results (in terms of MSE and other metrics). So I did a small test to check if GPy and GPflow were outputting the same values:

https://gist.github.com/beckdaniel/ac44ae59d223fde3e5c7

I expected the code above to give the same predict values but this is not happening. The model is supposed to be a standard GP regression with RBF kernel and zero mean function. No optimization is performed: all hyperparameters have their values kept at 1.0.

Did I miss something when I instantiated my GPflow model above? Is this behaviour expected?

In the Gaussian likelihood case I think there are instances where it would be good to work with the collapsed evidence lower bound of Titisias ( equation 16 of http://www.aueb.gr/users/mtitsias/papers/sparseGPv2.pdf ) rather than the uncollapsed "Big Data" bound used by Hensman et at ( http://auai.org/uai2013/prints/papers/244.pdf equation 4) . Perhaps we should consider implementing it?

Error is: ValueError: Cannot feed value of shape (4,) for Tensor u'Variable:0', which has shape '(5,)'
Error happens in predict_f - full error trace below.

Minimal working example below.
Issue on master branch as of April 24, 2016.

Code MWE - based on simple_regression.py

import GPflow
import tensorflow as tf
import os
import numpy as np
def getData():
    rng = np.random.RandomState( 1 )
    N = 30
    X = rng.rand(N,1)
    Y = np.sin(12*X) + 0.66*np.cos(25*X) + rng.randn(N,1)*0.1 + 3
    return X,Y
if __name__ == '__main__':
    X,Y = getData()
    k = GPflow.kernels.Matern52(1)
    meanf = GPflow.mean_functions.Linear(1,0)
    m = GPflow.gpr.GPR(X, Y, k, meanf)
    m.likelihood.variance = 0.01
    m._compile()
    print "Here are the parameters before optimization"
    print m
    m.kern.variance.fixed = True
    #m._compile() # If we compile again the code below works
    [mu,var] = m.predict_f(X)
    print mu
    print 'done'

Trace:

Traceback (most recent call last):
  File "/Users/mqbssaby/PrivateProjects/BranchedGP/runfile.py", line 29, in <module>
    [mu,var] = m.predict_f(X)
  File "/Users/mqbssaby/pythonlibs/GPflow/GPflow/model.py", line 82, in runnable
    return instance._session.run(graph, feed_dict=feed_dict)
  File "/Users/mqbssaby/anaconda/lib/python2.7/site-packages/tensorflow/python/client/session.py", line 340, in run
    run_metadata_ptr)
  File "/Users/mqbssaby/anaconda/lib/python2.7/site-packages/tensorflow/python/client/session.py", line 553, in _run
    % (np_val.shape, subfeed_t.name, str(subfeed_t.get_shape())))
ValueError: Cannot feed value of shape (4,) for Tensor u'Variable:0', which has shape '(5,)'

As per suggestion of @jameshensman it might be a good idea to start a page explaining how to contribute to GPflow .

Is there a way of getting the likelihood of the data given the model? The standard GPy ways seem not to work. I need to extract the likelihood of my regression model as a numeric value.

When fixing q_sqrt in SVGP i cannot compile the model as it breaks with the following error message

Traceback (most recent call last):
  File "/path/Debug.py", line 24, in <module>
    m._compile()
  File "/path/GPflow/GPflow/svgp.py", line 81, in _compile
    opt_step = GPModel._compile(self, optimizer)
  File "/path/GPflow/GPflow/model.py", line 154, in _compile
    f = self.build_likelihood() + self.build_prior()
  File "/path/GPflow/GPflow/svgp.py", line 112, in build_likelihood
    fmean, fvar = conditionals.gaussian_gp_predict_whitened(self._tfX, self.Z, self.kern, self.q_mu, self.q_sqrt, self.num_latent)
  File "/path/GPflow/GPflow/conditionals.py", line 100, in gaussian_gp_predict_whitened
    return conditional(Xnew, X, kern, q_mu, num_columns=num_columns, full_cov=full_cov, q_sqrt=q_sqrt, whiten=True)
  File "/path/GPflow/GPflow/conditionals.py", line 72, in conditional
    if q_sqrt.get_shape().ndims==2:
AttributeError: 'numpy.ndarray' object has no attribute 'get_shape'

Here is some code to reproduce the error:

import GPflow
import numpy as np
X = np.reshape(np.linspace(0, 1, 100), (100,1))
Y = np.reshape(np.linspace(0, 1, 100), (100,1))
Z = Y[:10]
m = SVGP(X, Y, GPflow.kernels.RBF(1), GPflow.likelihoods.Gaussian(), Z)
m.q_sqrt.fixed = True
m._compile()

I have had two cases in the last two weeks where are Codeship build has failed and then I have re-run it and it has worked fine.

In one case the error message was a server error which looked to have nothing to do with our code.

In the second case a command that is normally fine timed out suggesting that Codeship hung.

I would be willing to put up with this happening occasionally but I suggest that if others encounter similar problems we log it here and decide whether we think this acceptable given that there are other providers of build services on GitHub. In making such a decision we should also consider whether we have resolved #101 .

What would be the preferred way to sample from a GP with mean function/compositional kernel? I tired the various sample/predict and also the evalK() from #63 but I can't make it work.

When I attempt to perform GP regression on a dataset where dim(X)=17 and dim(y)=1, I get an exception as follows after about 40 iterations of the L-BFGS-B solver:

Traceback (most recent call last):
  File "/usr/local/lib/python3.5/dist-packages/tensorflow/python/client/session.py", line 655, in _do_call
    return fn(*args)
  File "/usr/local/lib/python3.5/dist-packages/tensorflow/python/client/session.py", line 637, in _run_fn
    status, run_metadata)
  File "/usr/lib/python3.5/contextlib.py", line 66, in __exit__
    next(self.gen)
  File "/usr/local/lib/python3.5/dist-packages/tensorflow/python/framework/errors.py", line 450, in raise_exception_on_not_ok_status
    pywrap_tensorflow.TF_GetCode(status))
tensorflow.python.framework.errors.InvalidArgumentError: Input matrix is not invertible.
     [[Node: MatrixTriangularSolve = MatrixTriangularSolve[T=DT_DOUBLE, adjoint=false, lower=true, _device="/job:localhost/replica:0/task:0/cpu:0"](Cholesky, sub_1)]]

However, GPy can compute this no problem (identical data, identical kernel setup), resulting in this:

Name : GP regression
Objective : -369.9790741292725
Number of Parameters : 11
Number of Optimization Parameters : 11
Updates : True
Parameters:
  GP_regression.                |              value  |  constraints  |  priors
  sum.std_periodic.variance     |    0.0202908485016  |      +ve      |        
  sum.std_periodic.period       |      961.078688147  |      +ve      |        
  sum.std_periodic.lengthscale  |      624.779412263  |      +ve      |        
  sum.rbf.variance              |  0.000746482801573  |      +ve      |        
  sum.rbf.lengthscale           |      2.94584474279  |      +ve      |        
  sum.mul.linear.variances      |     0.157037643245  |      +ve      |        
  sum.mul.Mat52.variance        |     0.157037643064  |      +ve      |        
  sum.mul.Mat52.lengthscale     |               (3,)  |      +ve      |        
  Gaussian_noise.variance       |    0.0137070123953  |      +ve      |  

It's not clear what is causing this to happen - the problem setup is essentially identical. It's frustrating because the regression itself is extremely good - when it computes. Am I missing something obvious here? Cheers.

Unless I've missed a thread we don't currently have the ability to fix a set of parameters all at once.

For instance it would be nice to say

model.kern.fixed=True

which would fix all the hyper parameters in the kernel, rather than going through fixing each one in turn which is what we currently have to do.

It would be a good idea to think about how we whiten matrices to stop the cholesky failing.

We often use something like:

A) chol( K + 1e-6 eye )

Is this actually what we want? We might for instance instead do

B) chol( K + 1e-6 eye * lambda)

where lambda takes into account the scale of K e.g lambda = mean(K) or somesuch like.

In any case even if we do want A) then we should consider pulling it out into a utlity function in my opinion. Thoughts?

Hi all,

I'm worried Autoflow still isn't working as expected. I'm on the latest commit. Below is an experiment where I add an extra inducing point. I think the bound calculation is good, but the predictions are different if I copy the parameters over to a new model.

Haven't had the chance to investigate thoroughly yet, but I thought I'd share.

[edit] More concise piece of code plus a workaround. The problem is that after adding the new inducing point, the prediction py2 is the same as the prediction py1. Only after deleting the prediction graph, do I get a correct prediction.

import numpy as np
import matplotlib.pyplot as plt
import GPflow
import datasets

Y, X, _, _ = datasets.load_data('snelson1d')
m = GPflow.sgpr.SGPR(X, Y, GPflow.kernels.RBF(1), Z=np.random.rand(2, 1)*np.max(X))

py1, pyv1 = m.predict_y([[3]])

Z = m.Z._array.copy()
m.Z._array = np.vstack((Z, [3]))  # Add an extra inducing point
m._compile()                               # Compile to ensure correct parameter mapping

py2, pyv2 = m.predict_y([[3]])

delattr(m, '_predict_y_graph')
py3, pyv3 = m.predict_y([[3]])

print py1
print py2
print py3

Hello there,

We are trying a very simple example and don't understand why the variance of a linear kernel is learnt to be 0. Here is what we do:

import GPflow
import numpy as np
import scipy as sp
import tensorflow as tf

# Simple training set
m = 2
n = 4
X = np.array([
[ 1.,  3.],
[ 2.,  4.],
[ 1.,  1.],
[ 2.,  1.]])
Y = np.array([1, 1, 0, 0])
Y = np.expand_dims(Y, axis=1)

# Build and train the model
kernel = GPflow.kernels.Linear(m)
probit = GPflow.likelihoods.Bernoulli()
model = GPflow.vgp.VGP(X=X, Y=Y, likelihood=probit, kern=kernel)
model.optimize()
# print model

# Predict a new data point
x_new = np.array([1.5, 5])
model.predict_y([x_new])

If we display the model, we see that model.kern.linear.variance is 0. Which leads to a predictive mean of 0.5 all the time.

Is that the expected behavior ?

Installing GPflow should be as easy as

pip install gpflow

via the python package index. It might also be nice to have an anaconda package.

1. Implement extensive scoping in the code. This will make the job of looking at the computational graph in tensorboard much easier. For example each kernel could define its own scope. We would need to agree on a systematic way of how to add scoping to the code.
2. Logging statements. To debug issues I have extensively been using the print commmand, for example to print the first 10 entries of a tensor and its shape:
   lik = tf.Print(lik,[pointToUpdate, tf.shape(lik), lik],message='point - lik=',name='likeddebug',summarize=10) In my code, I add Print nodes under a debug flag. Perhaps something similar would be useful for GPflow.
3. Use the summary statistics to log information. I have been playing with histogram_summary but tensorflow offers storing scalars and images as well. To keep a log of the value of a tensor simply do:
   tf.histogram_summary('probvector', tf.squeeze(probVector))
   Then collect all the summaries and write to file

tf.merge_all_summaries()
tf.train.SummaryWriter("logs",  self._session.graph)

This unit tests are taking too long on Travis as it stands leading to a build failure.

Another failure case of Autoflow:

import numpy as np
import GPflow
import tensorflow as tf

X = rnd.rand(10, 1)
Y = rnd.rand(10, 1)
mv = GPflow.sgpr.SGPR(X.copy(), Y.copy(), GPflow.kernels.RBF(1), Z=np.random.rand(4, 1)*np.max(X))
mf = GPflow.sgpr.GPRFITC(X.copy(), Y.copy(), GPflow.kernels.RBF(1), Z=np.random.rand(4, 1)*np.max(X))

mv.predict_y([[3]])
mf.predict_y([[3]])
mv.predict_y([[3]])

The cause: When running several models in parallel, AutoFlow shares its self.tf_args variable between models. So after mf.predict_y has been called, the AutoFlow object will now contain the tensorflow arguments for mf. These are the wrong placeholders when mv.predict_y is called again. The root of problems here, is that the decorator is called when the GPModel class is created, so the decorator is shared between all objects that inherit from it.

A quick fix would be to remove the caching mechanism. A different fix would be to store tf_args in instance instead of self.

I think the latter solution is preferable, so I'll submit a pull request for it.

I'm quite keen to get rid of the q_sqrt overparameterisation in SVGP. Currently a full square matrix is stored, even though only the bottom half is used. This has caused me some annoyances when optimising large amounts of inducing points. It slows down BFGS unnecessarily and makes the sizes of my parameter history files twice as large as they need to be.

I'm willing to put in some work to fix this, any tips on how to proceed? I could contribute to your current Tensorflow branch, or alternatively, the main Tensorflow branch. I think this is tangentially mentioned in #10.

At the moment we are able to merge a PR if it has no merge conflicts with the master branch.

The problem with this is that we run the tests on code that isn't actually what ends up master. Hence the PR tests could pass but master could fail.

We had a recent example of this when the unit tests became too long for Travis and this wasn't spotted before it took down master.

I wonder: Is it possible to get Travis to merge from master then run the tests? This would mean that what you see when you look at the results of tests is what you get.

It would be good to have continuous integration rather than checking the tests manually. It would also be good to run coverage at the same time. Actually at the time of writing the code coverage of the unit tests is very high.

Do the authors have any advice on integrating a GP from GPflow as part of a large graph? My two applications are: optimising a loss function with respect to predictive locations, and warping the inputs to the GP.

For the former, this seems to work:

    k = GPflow.kernels.Matern32(1)
    m = GPflow.gpr.GPR(X, Y, k)
    m.likelihood.variance = 0.01
    m.optimize(AdamOptimizer(0.1))
    init  = np.random.rand(3,)
    with m.tf_mode():
        xp = tf.Variable(init.reshape(-1,1))
        nu = np.min(Y)
        mean, cov = m.build_predict(xp, full_cov=True)
        m._session.run(tf.initialize_variables([p]))

but it requires me to call m.optimize() before constructing my further loss functions with build_predict, and I have to be careful about when and which variables I initialise since it seems to decouple the tf.Variables from the model, and give me garbage output after training.

For me, the ideal system would be a way to keep everything symbolic and inside whatever graph/session I am already using.

Xwarp = complex_function(X)
Xpwarp = complex_funciton(Xp)
k = GPflow.kernels.Matern32(2)
m = GPflow.gpr(Xwarp, Y, k)
mu, cov = m.predict(Xpwarp)
lik = m.likelihood()
loss = loss_function(mu, cov)
step = tf.train.AdamOptimizer(0.1).minimize(loss, [Xp])

Is there a way to do this in current GPflow?

On my desktop running the tests sequentially takes:

540s

Of this: 324s are spent in the mean functions test.

@jameshensman I don't have much familiarity with this test as it is new. Do you have any ideas about how we could keep coverage and reduce the time?

Obviously this would make a major difference to our build times.

@jameshensman mentioned that logging parameters during sampling/optimisation could be a useful feature to have. One suggestion was for the code to keep track of a list of dicts and return it after sampling.

What I've been using in my own code is an external logging class with a callback that is to be called by the optimiser. One class does more than just logging. I meant it to be a class that helps logging, storing and quitting long parallel optimisation runs, that are necessary for running lots of experiments for papers. I think it's nice, since it handles:

• Printing progress at intervals that can be adjusted by the user
• Logging parameters at intervals that can be adjusted by the user
• Storing the logged parameters to disk at intervals that can be adjusted by the user
• A timeout that kills optimisation if it takes too long

You can take a look at the current code here: https://gist.github.com/markvdw/216f6d4cfd040251e3b7e8c7551f559c. I'm happy to adapt it to fit into GPflow, if there's interest in more control over optimisation. It shouldn't be hard to get this to work with logging sampling as well.

This question is directed a @jameshensman specifically.

In the explanation of Titsias' method given in notebooks/SGPR notes.ipynb it says the following

Now to obtain a better conditioned matrix for inversion we rotate by L where LL^T = Kuu .

Could you explain a bit more about why the rotated matrix will be better conditioned please? I'm happy to add the answer to the notebook if you like.

I am creating lots of Models in a loop and would like to suppress the "compiling tensorflow function..." output it generates. I've modified _compile in model.py to contain a verbosity parameter that defaults to false, however in most instances the output is quite useful and I would like to have both options.

I am a bit lost in the autoflow decorator/param specifications...Where exactly is _compile called so that I can add an optional verbosity parameter?

I would like to reuse a same model for (a lot of) different data sets, but I believe in the current implementation GPflow needs compilation every time the data has changed.
( is it correct?)

I think TensorFlow itself has a framework to prevent such an unnecessary compilation, placeholder.
Do you guys have a plan to update GPflow to handle data as placeholder?

This is something that has been bugging me with GPflow, namely how fixed parameters are treated.

if you have a param.fixed = True then during
_compile-> make_tf_array
it is added in the graph as a tf.constant (see Param:make_tf_array)
However if we want to do somethings like
m.kern.param.fixed=True
m.kern.param = 2
m.optimize() -> this will trigger a compile that adds a constant
if we now do
m.kern.param = 3 this will also trigger a compile during the next optimize (see Param.setattr)
m.optimize()

However tensorflow offers another mechanism that will save some computation, namely placeholders.
As we have the free state, we could have a fixed state.
When _compile gets called it creates tf.Variables for the former and tf.PlaceHolders for the latter.

Then in our above scenario, recompilation would not be needed.
Note that this would apply to transforms and priors as well (for example we want to return the optimization with different limits on the logistic transform).

This will be only important when cost of recompilation is high, that is in complex models with large graphs.
Also note that by avoiding recompilation we do not pollute the default graph space in tensorflow with redundant graphs.

Perhaps this is not that important at this point but wanted to bring it to your attention for comments.