estimating parameters of gaussian process

[Linas Mockus]

Hi,

Can somebody suggest how to speed up the estimation of GP parameters? In the code below I am using Squared Exponential kernel with one discrete variable arch (coded as 0 or 1) and two continuous variables (wavelength & bandwidth). I wonder if using Cholesky will help.  Or any other tricks? Sorry for two many spaces in the code...

Thanks,
Linas

data {
    int <lower = 0> N;
    vector [N] response;
    vector [N] arch;
    vector [N] wavelength;
    vector [N] bandwidth;
}
transformed data {
}
parameters {
    real <lower = 0> sigmaf;
    real <lower = 0> sigman;
    real <lower = 0> l;
    vector [8] c;
}
model {
    matrix [N, N] Sigma;
    vector [N] mu;

    // off-diagonal elements
    for (i in 1 : (N - 1)) {
        for (j in i : N) {
            Sigma [i, j] <- sigmaf * exp (-(pow (arch [i] - arch [j], 2) +
                pow (wavelength [i] - wavelength [j], 2) +
                pow (bandwidth [i] - bandwidth [j], 2)) / (2 * l));
            Sigma [j, i] <- Sigma [i, j];
        }
    }
    // diagonal elements
    for (i in 1 : N)
        Sigma [i, i] <- sigmaf + sigman; // + jitter
    // priors
    sigman ~ cauchy (0, 2.5);
    sigmaf ~ cauchy (0, 2.5);
    l ~ cauchy (0, 2.5);
    c ~ normal (0, 100);

    for (i in 1 : N) mu [i] <- c [1] +
        c [2] * arch [i] +
        c [3] * wavelength [i] +
        c [4] * bandwidth [i] +
        c [5] * arch [i] * wavelength [i] +
        c [6] * arch [i] * bandwidth [i] +
        c [7] * wavelength [i] * bandwidth [i] +
        c [8] * arch [i] * wavelength [i] * bandwidth [i];
    response ~ multi_normal (mu, Sigma);
}
generated quantities {
}

[Mike Lawrence]

You can pre-compute all the distances in transformed parameters.

It seems that student_t(4,0,2.5) is being recommended these days over cauchy(0,2.5) as the former is relatively heavy-tailed but not so much so that it causes computational issues like the former.

That normal(0,100) prior is probably slowing you down; is there a way to rescale the inputs so you can have a normal(0,1)? Besides, while I'm relatively new to GPs, this is the first time I've seen a GP with a non-zero prior on the mean; normally c would be set in transformed data as a vector of zeros.

The cholesky form shown in the manual should definitely improve ess/time.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Michael Betancourt]

Precomputing the distances can certainly help, but because the 

covariance matrix depends on the parameters you have to compute

it at every iteration and there won’t be much benefit to computing

the Cholesky decomposition yourself (the multi_normal does a

Cholesky behind the scenes).  

There are some weak nonidentifiabilites in these kernels that we’ve

been debating, but the difference between that student_t and the

Cauchy will still be relatively small.  Similarly, having a non-zero

mean function shouldn’t cause too much trouble.

The immediate concern, however, is the discrete “arch” variable.

Trying to model a discrete covariate like that with a GP is typically

a bad idea.  Ultimately the more you cover the range of covariate

values the better identified your GP will be (and the better it will

fit), and with only 0 and 1 you’re not really constraining the hyperparameters

all that much.

[Linas Mockus]

Thanks a lot.

Are there any tricks to constrain hyperparameters? Maybe use GPstuff for initial estimates and then switch to rstan?

Linas

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/QAjcwWT9noc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Michael Betancourt]

That won’t help.  Remember we’re talking about nonidentifiabilities here --

the problem isn’t getting a good set of hyperparmeters it’s finding _all_ of

them.

[Linas Mockus]

What about using analytical expressions for joint predictive? I could make mu=0 but I am really interested having covariates arch, wavelength, and bandwidth. I doubt that I could make arch more than three levels.

Linas

[Mike Lawrence]

For the discrete parameter issue, could you maybe model the data as the sum of two latent GPs, where one GP forms the intercept and the other forms the arch effect (- this Gap for arch level 0, - this GP arch level 1)? 

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Andrew Gelman]

Often you can put in reasonable values for hyperparameters, just from the context of your problem, i.e. prior info,

[Mike Lawrence]

When you set mu=0 you will still be able to do inference on your variables' influence on the outcome. Because you've included these variables in the specification of the covariance matrix, the posterior samples on the mean function will be non-zero to the extent that the data call for it. 

On Tuesday, 12 July 2016, Linas Mockus <linasm...@gmail.com> wrote:

[Linas Mockus]

Maybe I should use mu=0 then? I like latent GPs - is it a mixture of two GPs or something else? Could you please point to some example?

It would be cool to set priors on hyperparameters. Any suggestions?

Linas

[Aki Vehtari]

Current Stan is slow for GPs due to big autodiff expression tree for covariance matrix operations, There is work in progress to implement internal covariance function which whould provide significant speed-up. The priors and mean functions doesn't have much effect on the speed, and Linas's code is otherwise fine except it would be better to use separate l for each covariate. Half-Cauchy prior can make HMC/NUTS slower as the covariance function parameters are weakly identified and the posterior is then likely to have a thivk tail, but changing half-Cauchy to something with less thick tail is not going to make this much faster.

On Tuesday, July 12, 2016 at 9:31:15 PM UTC+3, Michael Betancourt wrote:

The immediate concern, however, is the discrete “arch” variable.

Trying to model a discrete covariate like that with a GP is typically

a bad idea.

And I think one of the great properties is of GPs is that how easy it is to build hierarchical models by using discrete covariates! I've never had problems when combining discrete and continuous covariates.  The Linas's model above already does automatically the sum of two latent GPs, that Mike proposes :)

 

 Ultimately the more you cover the range of covariate

values the better identified your GP will be (and the better it will

fit), and with only 0 and 1 you’re not really constraining the hyperparameters

all that much.

For a single binary covariate the covariance matrix collapses to 2x2 covariance matrix. If exponentiated quadratic (aka squared exponential) is used it has two degrees of freedoms. I think this covariance matrix would be in many cases well identified and thus also the corresponding discerete GP ie normal distribution would be well identified, although there is dependency between sigmaf and l (and sigman) making them not so well identified. What I try to say is that you don't need to worry so much about the identifiability of hyperparameters and even less in discrete case :)

 

Aki

[Linas Mockus]

Thanks for all your input. I have modified the code to have separate ls for each covariate but it runs awfully slow. However, when I run two GP fits, one for arch=0, other for arch=1, the estimation time becames reasonable (maybe 10min or so). I wonder how to make sense from two GP models for different settings of discrete variables. And of course looking forward for new stan releases since as per David Duvenaud (http://www.cs.toronto.edu/~duvenaud/cookbook/index.html) GP should handle well the discrete variables (see How to use categorical variables in a Gaussian Process regression).

Linas

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/QAjcwWT9noc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Michael Betancourt]

In the degenerate 2x2 case the two parameters should be identifiable

(the variances constraining the marginal variance, the covariances

constraining the product of the marginal variance and length scale)

but they will still be nonlinearly correlated and without sufficient data

they can still be hard to fit!  At least in the world of _finite_ data.  ;-)

[Linas Mockus]

In this case there is bunch of data: 2X6X2 where first covariate is discrete and has two levels, while second and third covariates are continuous and have 6 and 2 levels correspondingly. Unfortunately design is unbalanced. In addition each set of covariates has 50 replicates - total 800 data points.

I am curious if encoding of discrete covariates matters: 0, 1 vs. -1, 1. Maybe scaling continuous covariates could help?

Thanks,

Linas

[Mike Lawrence]

For the variable that is continuous but you only measured it at 2 levels, I don't believe that there's any benefit in treating it as a second continuous dimension in the GP; indeed, it may even be computationally wasteful. My understanding is that with only 2 measurement points, you'll only ever find support for at max a linear effect anyway. So, really you have a single-dimension GP measured at 6 points crossed with two 2-level discrete variables. 

--

[Linas Mockus]

This is true, but I am estimating the model parameters and based on the estimates will try to figure out where to perform next experiment (sort of sequential DOE) to get another training point and hopefully improve the model. Ideally I would like to have GP model which includes discrete covariates. To circumvent computational difficulties I got two GP models for each level of discrete covariate but don't know how to make sense from those two models.

Linas