Logistic GP

[Aki Vehtari]

As a comparison to our Laplace approximated logistic GP <http://ba.stat.cmu.edu/abstracts/Vehtari.php>, we tested MCMC, too. Last autumun I asked from Betancourt why HMC is very slow for that model, and he suggested testing simplex transformation. As Stan has simplex implemented, I decided to test this first with Stan. Here's a code and example data. It's slow even with just 21 data points. For me it was interesting just to test if I can write this model with Stan, and I don't need for anything, but I'd be interested to learn why it is slow.

data {

  int<lower=1> N;     // number of bins

  vector[N] x;        // bin centers

  int<lower=0> y[N];  // number of observations per bin

}

parameters {

  simplex[N] p;       // densities

  real<lower=0> ms2;  // magnitude of GP

  real<lower=0> l;    // lengthscale of GP

  real b1;            // mean function parameter

  real<lower=0> b2;   // mean function parameter

}

transformed parameters {

  vector[N] mu;       // mean function

  vector[N] f;        // latent function

  mu <- b1*x - b2*x .* x;

  f <- log(p);        

}

model {

  // covariance matrix

  matrix[N,N] Sigma;

  for (i in 1:N) 

    for (j in 1:N)

      Sigma[i,j] <- ms2*exp(-pow(x[i] - x[j],2)/l)

                    + if_else(i==j, 0.0001, 0.0);

  for (i in 1:N)

    for (j in (i+1):N)

      Sigma[j,i] <- Sigma[i,j];

  // GP prior

  f ~ multi_normal(mu,Sigma);

  // priors for covariance function parameters

  ms2 ~ student_t(4,0,10);

  l ~ student_t(4,0,4);

  // Jacobian

  increment_log_prob(- f);

  // log-prob

  for (n in 1:N)

    increment_log_prob(y[n]*f[n]);

}

I use MatlabStan, but I guess you can understand following commands and translate to R

nn=21;

xx=[-10:10];

yy=[0 0 0 0 1 2 3 3 2 4 4 5 3 4 2 0 1 1 0 0 0];

dat=struct('y',yy,'N',nn,'x',xx);

fit = stan('model_code',code,'data',dat,'sample_file','logitgp','file_overwrite',true,'verbose',true);

samples = fit.extract('permuted',true);

plot(xx,mean(samples.p),xx,yy./sum(yy))

Aki

[Herra Huu]

Simplex transformation doesn't seem to do that much here, in fact I think the sampling become even a bit slower. But using reparameterizations "tricks" for multivariate normal distribution (p. 182-184 in Stan manual) did help a lot (about 4s / 2000 samples compared to 100s with the original model). I'm still learning about Gaussian processes, so I'm not 100% sure if the model is correct though (results seems to be ok). But anyway, here is the code:

data {
  int<lower=1> N;     // number of bins
  vector[N] x;        // bin centers
  int<lower=0> y[N];  // number of observations per bin
}
parameters {

  real<lower=0> ms2;  // magnitude of GP
  real<lower=0> l;    // lengthscale of GP
  real b1;            // mean function parameter
  real<lower=0> b2;   // mean function parameter

  vector[N] alpha;

}
transformed parameters {
  vector[N] mu;    

  vector[N] p;
  matrix[N,N] L;
  matrix[N,N] Sigma;

  mu <- b1*x - b2*x .* x;

  for (i in 1:N)
    for (j in 1:N)
    {
      Sigma[i,j] <- ms2*exp(-pow(x[i] - x[j],2)/l) + if_else(i==j, 0.0001, 0.0);

      Sigma[j,i] <- Sigma[i,j];
    }

  L <- cholesky_decompose(Sigma);
  p <- softmax(mu + L * alpha);       
}
model {
  alpha ~ normal(0,1);

  // priors for covariance function parameters
  ms2 ~ student_t(4,0,10);
  l ~ student_t(4,0,4);

  // log-prob
  for (n in 1:N)

    increment_log_prob(y[n]*log(p[n]));
}

[Bob Carpenter]

Thanks for reporting back.

The quantity to measure is effective sample size per second, not just
iteration speed.

Various reparameterizations affect both iterations per second or
effective samples per iteration, and it's the net effect that matters.

- Bob


On Mar 12, 2014, at 9:23 PM, Herra Huu <heikki.uu...@gmail.com> wrote:

> Simplex transformation doesn't seem to do that much here, in fact I think the sampling become even a bit slower. But using reparameterizations "tricks" for multivariate normal distribution (p. 182-183 in Stan manual) did help a lot (about 4s / 2000 samples compared to 100s with the original model). I'm still learning about Gaussian processes, so I'm not 100% sure if the model is correct though (results seems to be ok). But anyway, here is the code:
>
> stan_code <- 'data {

> int<lower=1> N; // number of bins
> vector[N] x; // bin centers
> int<lower=0> y[N]; // number of observations per bin
> }
> parameters {

> real<lower=0> ms2; // magnitude of GP
> real<lower=0> l; // lengthscale of GP
> real b1; // mean function parameter
> real<lower=0> b2; // mean function parameter

> vector[N] alpha;

> }
> transformed parameters {
> vector[N] mu;

> vector[N] log_p;

> matrix[N,N] L;
> matrix[N,N] Sigma;
>

> mu <- b1*x - b2*x .* x;

> for (i in 1:N)
> for (j in 1:N)
> {
> Sigma[i,j] <- ms2*exp(-pow(x[i] - x[j],2)/l) + if_else(i==j, 0.0001, 0.0);

> Sigma[j,i] <- Sigma[i,j];
> }

> L <- cholesky_decompose(Sigma);
> log_p <- log_softmax(mu + L * alpha);
> }
> model {
> alpha ~ normal(0,1);
>

> // priors for covariance function parameters
> ms2 ~ student_t(4,0,10);
> l ~ student_t(4,0,4);
>

> // log-prob
> for (n in 1:N)

> increment_log_prob(y[n]*log_p[n]);
> }'
>
> --
> 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.
> For more options, visit https://groups.google.com/d/optout.

[Herra Huu]

Yes, this is definitely true. Should have stated in the original post, that n_eff estimates for the parameters seems to a bit higher as well.

[Aki Vehtari]

Huu: Thanks for the Cholesky version. I knew the trick, but assumed that with n=20, it would not make so big difference.

Bob: Why in Huu's version there is softmax but in your posting the quotation of Huu's code has log_softmax? Did you edit that?

Huu's new code is still impractically slow with n=40.

The n_eff depends very much on much "jitter" is added to the diagonal

 > if_else(i==j, 0.0001, 0.0)

as larger jitter makes the latent function less constrained by softmax.  The Markov chain is run in n-dimensional space (ignoring now the covariance function parameters), but the latent function is constrained in n-1 dimensional manifold. Jitter allows the Markov chain to move a little bit away from the manifold. Iideally the amount of jitter should be zero, but for example, if the jitter is set to 0.0000001, n_eff is approximately 20 with n=4000.

In autumn Michael said that had obtained huge speed-up by using simplex in some of his models and if I remember correctly the reason would be because then the chain is also in n-1 dimensional space.

I could not figure out how to use simplex with Cholesky version, as now alpha~normal(0,1) is the sampled parameter...

Aki

keskiviikko, 12. maaliskuuta 2014 23.04.34 UTC+2 Herra Huu kirjoitti:

Yes, this is definitely true. Should have stated in the original post, that n_eff estimates for the parameters seems to a bit higher as well.


keskiviikko, 12. maaliskuuta 2014 22.44.02 UTC+2 Bob Carpenter kirjoitti:

[Herra Huu]

My original post had the log_softmax version of the code, but I edited it afterwards so that the Matlab code for the plots etc would work directly.

[Michael Betancourt]

Aki,

Let me be a little more detailed about the importance of the simplex.  We do indeed have a Markov chain sampling from an n-dimensional space but the important part is that the latent function, and hence the data, are only sensitive

to an (n - 1) dimensional submanifold.  That means that directions transverse to that submanifold are unconstrained, and the chain will just drift around in a random walk.  The simplex allows us to enforce the constraint analytically

so we can sample from the (n - 1) dimensional space and avoid the performance-reducing drift.

That said we have to be really careful here.  Firstly, the simplex parameters are already defined to be positive -- recall that a simplex is a collection { p_i } satisfying p_i > 0 and sum_i p_i = 1 -- so the f = log℗ is actually making your

posterior worse.  Moreover, the GP shouldn't be in the n dimensional observation space but rather the (n - 1) dimensional latent space.  So what you want is something like

f ~ GP( mean, covariance function)

p ~ simplex_map(f)

y ~ binomial℗

where simplex_map is a map from the latent space to the observation space, R^{n - 1} -> R^{n}.  Unfortunately this isn't exposed and I'm not sure if there's an easy way to do this right now in Stan. Let me start a separate thread

for that.

-Mike

On Thu, Mar 13, 2014 at 5:46 PM, Herra Huu <heikki.uu...@gmail.com> wrote:

My original post had the log_softmax version of the code, but I edited it afterwards so that the original Matlab code for the plots etc would work directly.

torstai, 13. maaliskuuta 2014 18.16.40 UTC+2 Aki Vehtari kirjoitti:

Bob: Why in Huu's version there is softmax but in your posting the quotation of Huu's code has log_softmax? Did you edit that?

[Bob Carpenter]

On Mar 14, 2014, at 12:19 PM, Michael Betancourt <betan...@gmail.com> wrote:

> Aki,
>
> Let me be a little more detailed about the importance of the simplex. We do indeed have a Markov chain sampling from an n-dimensional space but the important part is that the latent function, and hence the data, are only sensitive
> to an (n - 1) dimensional submanifold. That means that directions transverse to that submanifold are unconstrained, and the chain will just drift around in a random walk. The simplex allows us to enforce the constraint analytically
> so we can sample from the (n - 1) dimensional space and avoid the performance-reducing drift.
>
> That said we have to be really careful here. Firstly, the simplex parameters are already defined to be positive -- recall that a simplex is a collection { p_i } satisfying p_i > 0 and sum_i p_i = 1 -- so the f = log℗ is actually making your
> posterior worse. Moreover, the GP shouldn't be in the n dimensional observation space but rather the (n - 1) dimensional latent space. So what you want is something like
>
> f ~ GP( mean, covariance function)
> p ~ simplex_map(f)
> y ~ binomial℗

If p is a simplex, then you probably want a multinomial or
categorical for y.

>
> where simplex_map is a map from the latent space to the observation space, R^{n - 1} -> R^{n}. Unfortunately this isn't exposed and I'm not sure if there's an easy way to do this right now in Stan. Let me start a separate thread
> for that.

The problem is that our softmax() is defined on R^n -> R^n, so it takes
a redundant parameterization for the input. Now what you can do is use
a point in R^(n-1) and set one of the values to 0.

But we should really have a softmaxKminus1() function or something like that
that works on R^(n-1) -> R^n. And perhaps the inverse_softmaxKminus1() that
goes the other way around.

This is all related to the discussions we had with Andrew about whether to
use K or (K-1) dimensions. I believe the issue with using (K-1) is that
the prior becomes asymmetric. I seem to recall that Andrew preferred the
K-ary version, which is why it's the way it is. In retrospect, I think we
should've done it the other way around. We can just let users do it either
way using a new softmax function defined on (n-1).

If you do it with K parameters, you need stronger priors to keep the estimates
well-behaved so they don't just drift around as Michael points out they will do.

- Bob

[Herra Huu]

perjantai, 14. maaliskuuta 2014 15.49.20 UTC+2 Bob Carpenter kirjoitti:

The problem is that our softmax() is defined on R^n -> R^n, so it takes
a redundant parameterization for the input.  Now what you can do is use
a point in R^(n-1) and set one of the values to 0.  

i was trying to think how could we use this approach here with GP, but couldn't figure it out. The problem is that as in the original model we have f ~ GP( mean(x), covariance function(x)), where x is n-dimensional. And we would like to set one of the component of f to zero, say f[1]=0. But how can we take into account the covariance between f[1] and the rest of the f[i]? And basically what would be the right distribution for f[2]...f[n] corresponding to the original model?

 

[Michael Betancourt]

There is no unique way to do that, that's why you want to model the true (n -1)-dimensional latent space directly.

On Fri, Mar 14, 2014 at 4:41 PM, Herra Huu <heikki.uu...@gmail.com> wrote:

perjantai, 14. maaliskuuta 2014 15.49.20 UTC+2 Bob Carpenter kirjoitti:

The problem is that our softmax() is defined on R^n -> R^n, so it takes
a redundant parameterization for the input.  Now what you can do is use
a point in R^(n-1) and set one of the values to 0.  

i was trying to think how could we use this approach here with GP, but couldn't figure it out. The problem is that as in the original model we have f ~ GP( mean(x), covariance function(x)), where x is n-dimensional. And we would like to set one of the component of f to zero, say f[1]=0. But how can we take into account the covariance between f[1] and the rest of the f[i]? And basically what would be the right distribution for f[2]...f[n]?

 

[Herra Huu]

Is there some good references for this latent variable approach? It kinda makes sense, but I'm missing something. For example, shouldn't the mean and covariance functions be changed as well? Ie. in the original model the mean function is R^{n} -> R^{n}, in the latent version shouldn't it be R^{n} -> R^{n-1} and so on?

[Michael Betancourt]

I haven't seen many references on these things -- at least not in any general review.

In the latent approach you're modeling the latent space, so everything should be in (n - 1)-dimensions, especially if you're doing a non-parameteric model with little prior information.

That said, if mu has some physical meaning then you can define it on the n-dimensional space and then compose it with the latent mapping, mu( latent_to_simplex( x ) )

to implicitly define a R^{n - 1} -> R^{n} function.

On Fri, Mar 14, 2014 at 5:08 PM, Herra Huu <heikki.uu...@gmail.com> wrote:

Is there some good references for this latent variable approach? It kind makes sense, but I'm missing something. For example, shouldn't the mean and covariance functions be changed as well? Ie. in the original model the mean function is R^{n} -> R^{n}, in the latent version shouldn't it be R^{n} -> R^{n-1} and so on?

perjantai, 14. maaliskuuta 2014 18.45.13 UTC+2 Michael Betancourt kirjoitti: