Slow imputation of ordered variables with HMC, variational bayes?

[Guido Biele]

 

Hello,

As it cannot be said often enough: 

Thanks to the Stan-team for the time and effort

they invest to make this a great list!!

 

Now to my problem:

I am implementing fully conditional imputation of amongst

others ordinal variables. I validated the model (code

below) with simulated data, but now it turns out that

fitting is very slow when using HMC for real data.

When using a subset of my data -2500 cases, 4 variables,

two of which are ordinal, with 4 and 6 levels, respectively-

fewer than 100 iterations are done after 12h (with 

adapt_delta=.8, on a PC with a 3.4 ghz processor and 8gb RAM).

The same data and model can be fit in around a minute

using variational Bayes, but then vb is “experimental”.

 

My questions:

1) Does anybody see a way to speed up the model?

I assume the speed-bottle neck is marginalizing out multiple

(here 2) ordinal variables.

2) Is there a reasonable approach to determine that vb

results can be trusted?

I am thinking of (a) simulating data with the same

characteristics as my real data and validate vb by

comparing true and imputed data and/or (b) taking multiple

(very) small subsets of my real data and checking that

for those hmc and vb give the same results.

 

Thanks in advance for your help!

Guido

 

/* 

fully conditional imputation model for mixed normal & ordinal data.

model assumes two columns of ordinal and two of normal data

*/

data {

  int N;                       // number cases (rows)

  int K;                       // number variables (columns)

  int missings[N,K];     // indicator and index of missing variables

  int Nmiss[K];            // number missing variables

  int nD;                     // number of discrete variables

  int DCs[nD];             // number of unique values per column (discrete vars)

  matrix[N,K-nD] x;     // matrix with continuous variables

  int Dx[N,nD];            // matrix with discrete variables

}

transformed data {

  int nC;                  // number continuous variables

  int idvs1[N];          // number of possible values for 1st categorical variable

  int idvs2[N];          // number of possible values for 2nd categorical variable

  int n_lps[N];          // number of combinations of imputed categorical values over which to marginalize

  nC = K-nD;

  for (n in 1:N) {

    idvs1[n] = missings[n,1] == 0 ? 1 : DCs[1];

    idvs2[n] = missings[n,2] == 0 ? 1 : DCs[2];

    n_lps[n] = idvs1[n]*idvs2[n];

  }

}

parameters{

  matrix[K,K-1] betas;                   // regression weights for imputation  

  vector[K] alpha;                         // intercepts for imputation

  vector<lower=0>[nC] sigmas;     // error variance for imputation

  ordered[DCs[1]-1] cp1;               // cutpoints first discret variable

  ordered[DCs[2]-1] cp2;          

  

  real m1[Nmiss[1]];     // missing values

  real m2[Nmiss[2]];

  real m3[Nmiss[3]];

  real m4[Nmiss[4]];

}

model{

  matrix[N,nC] xi;  // matrix with observed and imputed continuous values

  to_vector(betas) ~ normal(0,2);

  alpha ~ normal(0,2);

  sigmas ~ cauchy(0,2);

  m1 ~ normal(0,2);                  // model won't converge without priors  

  m2 ~ normal(0,2);                  // on imputed discrete values

  

  // fill in missing continuous values

  for (n in 1:N) {

    xi[n,1] = (missings[n,3] == 0 ? x[n,1] : m3[missings[n,3]]);

    xi[n,2] = (missings[n,4] == 0 ? x[n,2] : m4[missings[n,4]]);

  }

  for (n in 1:N) {

    int iDx1[n_lps[n]];      // instantiations of of imputed discrete variables

    int iDx2[n_lps[n]];

    real lp_dv[n_lps[n]];   // log-probabilies for combinations of imputed discrete variables

    int m;

    m = 0;

    for (idv1 in 1:idvs1[n]) {

      for (idv2 in 1:idvs2[n]) {

        m = m+1;

        iDx1[m] = missings[n,1] == 0 ? Dx[n,1] : idv1;

        iDx2[m] = missings[n,2] == 0 ? Dx[n,2] : idv2;

        lp_dv[m] = (missings[n,1] == 0 ? 0 : ordered_logistic_lpmf(idv1 | m1[missings[n,1]], cp1)) + 

                   (missings[n,2] == 0 ? 0 : ordered_logistic_lpmf(idv2 | m2[missings[n,2]], cp2));

       }

     }

    for (k in 1:K) {

      real lps[n_lps[n]];

      // marginalize out imputed discrete values 

      for (j in 1:n_lps[n]) {

        real eta;

        if (k == 1) {

          eta = alpha[k] + iDx2[j]*betas[k,1] + xi[n,1]*betas[k,2] + xi[n,2]*betas[k,3];

          lps[j] = lp_dv[j] + ordered_logistic_lpmf(iDx1[j] | eta, cp1);

        } else if (k == 2) {

          eta = alpha[k] + iDx1[j]*betas[k,1] + xi[n,1]*betas[k,2] + xi[n,2]*betas[k,3];

          lps[j] = lp_dv[j] + ordered_logistic_lpmf(iDx2[j] | eta, cp2);

        } else if (k == 3) {

          eta = alpha[k] + iDx1[j]*betas[k,1] + iDx2[j]*betas[k,2] + xi[n,2]*betas[k,3];

          lps[j] = lp_dv[j] + normal_lpdf(xi[n,1] | eta, sigmas[1]);

        } else if (k == 4) {

          eta = alpha[k] + iDx1[j]*betas[k,1] + iDx2[j]*betas[k,2] + xi[n,1]*betas[k,3];

          lps[j] = lp_dv[j] + normal_lpdf(xi[n,2] | eta, sigmas[2]);

        }

      }

      target += log_sum_exp(lps);

    }

  }

}

generated quantities {

  int i1[Nmiss[1]];

  int i2[Nmiss[2]];

  for (j in 1:Nmiss[1]) i1[j] = ordered_logistic_rng(m1[j], cp1);

  for (j in 1:Nmiss[2]) i2[j] = ordered_logistic_rng(m2[j], cp2);

}

[Bob Carpenter]

> On Dec 9, 2016, at 4:20 AM, Guido Biele <guido...@gmail.com> wrote:
>
>
> Hello,
>
> As it cannot be said often enough:
>
> Thanks to the Stan-team for the time and effort
>
> they invest to make this a great list!!
>
>
> Now to my problem:
>
> I am implementing fully conditional imputation of amongst
>
> others ordinal variables. I validated the model (code
>
> below) with simulated data, but now it turns out that
>
> fitting is very slow when using HMC for real data.
>
> When using a subset of my data -2500 cases, 4 variables,
>
> two of which are ordinal, with 4 and 6 levels, respectively-
>
> fewer than 100 iterations are done after 12h (with
>
> adapt_delta=.8, on a PC with a 3.4 ghz processor and 8gb RAM).

>
> The same data and model can be fit in around a minute
>
> using variational Bayes, but then vb is “experimental”.
>
>
> My questions:
>
> 1) Does anybody see a way to speed up the model?

I can help with speedig up the computations. You want to look
at other things going on to see if there's some way to parameterize
it better.

Overall, our ordered logistic isn't the fastest thing in the world,
especially if you can't vectorize it.

The biggest saving you can get in Stan speed is saving any intermediate
computations you do in order to reuse them. For that, you have to
comb over your program and find things in loops that don't depend on
the loop index and pull them out. And try to vectorize as much as possible.
So if you have a bunch of ordered logistics with the same cutpoints, save
up their arguments if you can and make one call (not sure this is possible).

And don't store things in lp_dv that are 0 and then put them
into other expressions. Just push the branching logic down and
only create and multiplyw hen needed.

For clarity, you don't ever want to compare to 0 or 1 in a conditional.
Instead, convert

missings[n,1] == 0 ? Dx[n,1] : idv1

into

missings[n, 1] ? idv1 : Dx[n, 1];

You can probably do things like this:

ta = alpha[k] + iDx2[j]*betas[k,1] + xi[n,1]*betas[k,2] + xi[n,2]*betas[k,3];

with matrix algebra, which might be faster.

>
> I assume the speed-bottle neck is marginalizing out multiple
>
> (here 2) ordinal variables.
>
>
>
> 2) Is there a reasonable approach to determine that vb
>
> results can be trusted?

Not that we know of other than comparing them to HMC or other
reference output.

> I am thinking of (a) simulating data with the same
>
> characteristics as my real data and validate vb by
>
> comparing true and imputed data and/or

You want to do this with interval posterior coverage,
not with point estimation. We do know that by construction
VB tends to underestimate uncertainty (the construction is
the direction of KL divergence being minimized).

> (b) taking multiple
>
> (very) small subsets of my real data and checking that
>
> for those hmc and vb give the same results.

Same results here means same posterior, not just posterior means.

- Bob

[Guido Biele]

Thanks for your quick repsonse Bob!

I have a 2 short follow up questions:

1) Is it possible that no vectorized version of the ordered_logistic_ ... 

    functions is available?

2) I implemented matrix multiplication, but this is a bit tricky because

    some of my variables are discrete and others are continuous.

For example, this line:

eta = alpha[k] + iDx1[j]*betas[k,1] + iDx2[j]*betas[k,2] + xi[n,2]*betas[k,3];

became this 

eta = alpha[k] + append_col(to_row_vector(iDx[,j]),xi[n,3]) * betas[,k];

Now I am not sure if the new version would be faster or slower, because

of the usage of append_col and to_row_vector.

Best - Guido

[Bob Carpenter]

> On Dec 12, 2016, at 9:53 AM, Guido Biele <guido...@gmail.com> wrote:
>
> Thanks for your quick repsonse Bob!
>
> I have a 2 short follow up questions:
>
> 1) Is it possible that no vectorized version of the ordered_logistic_ ...
> functions is available?

Not yet. It's on our to-do list. The framework we
set up to vectorize everyhting else doesn't work for
non-scalar parameters, so we need to figure out how to
generalize.

> 2) I implemented matrix multiplication, but this is a bit tricky because
> some of my variables are discrete and others are continuous.
> For example, this line:
> eta = alpha[k] + iDx1[j]*betas[k,1] + iDx2[j]*betas[k,2] + xi[n,2]*betas[k,3];
> became this
> eta = alpha[k] + append_col(to_row_vector(iDx[,j]),xi[n,3]) * betas[,k];
>
> Now I am not sure if the new version would be faster or slower, because
> of the usage of append_col and to_row_vector.

The append_col and to_row_vector just copy. The bottleneck is
almost always in the size of the autodiff graph, so usually it's
a win to aggregate operations. The original line was eta[j], right?
And I didn't follow why the k changed places in betas. But in general,
it's usually worth appending and reshaping if you can use matrix
operations. Sometimes you can redeclare variables so that the reshaping
is easy.

- Bob

[Guido Biele]

OK, then I am looking forward to vectorization for non-scalar parameters.

I also had already changed the model to use vectorization wherever possible,

and to compute the marginalization operatons only for rows with missing

discrete data. This made things noticeably faster.

eta was always just eta and k changed its place in betas because i

changed the dimensions of betas.

about eta: 

right now eta is intialized as a single real value and calculated within

the loop for the marginalization.  e.g.

// loop throug all possible combinations of discrete variable values

for (j in 1:n_lps[n]) {
        real eta;

       ...

        eta = alpha[k] + append_col(to_row_vector(iDx[,j]),xi[n,3]) * betas[,k];

        lps[j] = lp_dv[j] + normal_lpdf(xi[n,4] | eta, sigmas[2]);
}

Instead, I could again use matrices to calculate

all etas outside the loop:

real etas[n_lps[n]];
matrix[n_lps[n],K-1] tmp_x;

tmp_x = append_col(append_col(iDx1,IDX2), rep_vector(xi[n,4],[n_lps[n])
etas = alpha[k] + tmp_x *  betas[,k];

for (j in 1:n_lps[n]) {

  ....
  lps[j] = lp_dv[j] + normal_lpdf(xi[n,4] | etas[j], sigmas[2]);
}

would this be faster?

I do not understand well what goes into the construction of the autodiff graph.

Is there a section in the manual that explains this, so that I can figure

out myself what changes to the model most likely impove speed?

Guido

[Daniel Lee]

Hi Guido,

You want to read the paper on autodiff.

http://mc-stan.org/documentation/

Look under papers for "The Stan Math Library: Reverse-Mode Automatic Differentiation in C++."

Daniel

--
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.

[Guido Biele]

Thank you Daniel,

I had hoped for a shorter way to learn this,

but let me see if I can take a short cut :-).

Figure 1 in the paper shows an expression graph.

Is it correct to say that if I have column and row 

vectors x and b of length 3, then if I write anywhere

in the model (or transformed parameters) block 

eta = x[1]*b[1] + x[2]*b[2] + x[3]*b[3];

this results in 11 automatic differentiation variables,

and an expression tree with 3 levels, whereas 

eta = x * b;

results in only 3 automatic differentiation variables

and an expression tree of 2 levels?

Guido

[Daniel Lee]

Hi Guido,

On Wed, Dec 14, 2016 at 9:28 AM, Guido Biele <guido...@gmail.com> wrote:

Thank you Daniel,

I had hoped for a shorter way to learn this,

but let me see if I can take a short cut :-).

Ha... there isn't really a short cut here.

At some point, you either have to spend the time to understand autodiff if you're serious about trying to optimize runtime. Or just follow along with best practices and stop trying to squeeze out maximal performance. (Honestly, I'd suggest the second option any day.)

 

Figure 1 in the paper shows an expression graph.

Is it correct to say that if I have column and row 

vectors x and b of length 3, then if I write anywhere

in the model (or transformed parameters) block 

eta = x[1]*b[1] + x[2]*b[2] + x[3]*b[3];

this results in 11 automatic differentiation variables,

and an expression tree with 3 levels, whereas 

eta = x * b;

results in only 3 automatic differentiation variables

and an expression tree of 2 levels?

No. It's much more complicated than that. And this is why Stan is fast. It depends whether eta, x, and b are parameters or data. Your original Stan program was too many posts before, so I don't remember what they are.

There's a large section in the paper about how we actually implement the C++ so that we cut out lots of the expression graph. So... no, you're not right here.

If you were to write the expressions in C++ directly without using built in functions, then yes, you're right. (in the second expression, eta = x * b, you'd have to write eta = stan::math::multiply(x, b), but yes).

Daniel

 

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Guido Biele]

Thanks Daniel, 

this helps (not as in "now I fully understand" but more

as in "now I know a bit better what I don't know").

Given my limited math and programming background, 

its probably wise I stick to the "follow best paractice" option.

[about what eta, x and b are: eta and b are parameters, and

x is a mix of parameters (used to imputed missing data) and

data.]

Guido

[Daniel Lee]

Hi Guido,

On Wed, Dec 14, 2016 at 10:28 AM, Guido Biele <guido...@gmail.com> wrote:

Thanks Daniel, 

this helps (not as in "now I fully understand" but more

as in "now I know a bit better what I don't know").

Given my limited math and programming background, 

its probably wise I stick to the "follow best paractice" option.

[about what eta, x and b are: eta and b are parameters, and

x is a mix of parameters (used to imputed missing data) and

data.]

If that's what it is, the fastest you can do it is to split x into the parameters and the data and do those operations separately. Then you lose the generality of jamming everything into one vector, but that's really going to be faster.

I wouldn't suggest it in general, though. It makes it complicated and hard to code. There's a tradeoff between fast computation and fast / robust coding. I'd rather code things that are slower than I'm less likely to make a mistake than to code something super fast and wrong.

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

Yes.

> I do not understand well what goes into the construction of the autodiff graph.
> Is there a section in the manual that explains this, so that I can figure
> out myself what changes to the model most likely impove speed?

Only hints. The short story is that when the log density is
evaluated, the expression graph is saved. This stores a data
structure for the result and arguments of every mathematical
operation (evaluation order gives you a nice topologically sorted
stack and proof the graph is directed and acyclic).

The graph is then used to calculate derivatives by passing
the chain rule through it backwards from teh result (hence
the "reversse" in our reverse-mode autodiff).

The derivative passing usually consumes well over 90% of the
compute time in all of Stan's algorithms.

Derivative passing is very slow as it involves walking the
stack of subexpressions and calling a virtual function for
each one. Memory locality is difficult to achieve for anything
other than matrix operations.

The smaller the expression graph is, the faster Stan is.

Almost all of the optimization tips we provide have the sole aim
of making the expression graph smaller.

Every time you execute an operation, it goes on the expression
graph. So the best thing you can do is save intermediate
results and reuse them. It automatically translates into
efficient dynamic programming on the back end (the stack and
its topological sort show it works).

Vectorized probability functions address both points. It
creates a single result node, so fewer virtual function calls,
and it automatically saves the results of intermediate
computations, for more efficiency. For example, when we
do y ~ normal(mu, sigma); for vector y and scalars mu and sigma,
we only need to evaluate log(sigma) once, which is hugely
expensive compared to simple arithmetic.

Matrix operations do the same thing. With scalar * vector,
you can collaps all the individual scalar updates to a single
update, saving virtual function calls and space again. Matrix
times matrix is even more cost saving.

This is all just computational efficiency. Reparameterization,
tighter priors, etc., address statistical efficiency, which can
be even more important for overall speed.

The manual chapter on efficiency and optimization addresses
both of these kinds of optimization.

- Bob

[Bob Carpenter]

Better yet, here's what we do and what we recommend everyone
else do:

1. code the model as straightforwardly as possible

2. generate fake data

3. make sure the program properly codes the model

4. run the program on real data

5. *If* the model is too slow, optimize *one step at a time*
and for each step, go back to (3).

The optimizations can be of either the statistical or the
computational variety. Slow iterations can be due to
computational statistical problems with parameterization (requiring
too many iterations) or due to slow code (each log density and
derivative evaluation being too slow).

We'd like to do a better job of helping people figure out how to
do the statistical optimization. It's usually more important and
also more variable based on the model. And it's harder to diagnose
because you need to look within iterations and at posterior
draws.

- Bob

[Guido Biele]

I agree with what you are writing Bob.

I went through steps 1-4, and I am also confident that the 

parameterization and priors are OK.

The core problem is just that I have to marginalize over the 

ordered logistic regression, and this just involves lots of 

computations. I see if I find other places where I can improve

the model and will post it again than (In case this could be useful

for somebody else).

Guido

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/zoelNcwfFYk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Guido Biele]

just to make sure that I understand what your are saying:

lets say we have a design matrix for a regression model x_imputed, which consists mostly of observed data, but some entries of which are parameters (i.e. imputed).

If we now focus on one row of x and want multiplicate x_imputed[n,] with regressions weights beta, are you saying that then it would be fastest to do straight forward matrix multiplication for the rows that consist only of data, and to separate data from parameters in a rows of x_imputed that contains parameters before multiplying each predictor with its weight beta?

guido

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/zoelNcwfFYk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Daniel Lee]

On Wed, Dec 14, 2016 at 3:46 PM, Guido Biele <guido...@gmail.com> wrote:

just to make sure that I understand what your are saying:

lets say we have a design matrix for a regression model x_imputed, which consists mostly of observed data, but some entries of which are parameters (i.e. imputed).

If we now focus on one row of x and want multiplicate x_imputed[n,] with regressions weights beta, are you saying that then it would be fastest to do straight forward matrix multiplication for the rows that consist only of data, and to separate data from parameters in a rows of x_imputed that contains parameters before multiplying each predictor with its weight beta?

Yes.

 

guido

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

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/zoelNcwfFYk/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
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+unsubscribe@googlegroups.com.

[Guido Biele]

After more experimentation I came to the conclusion that this approach (marginalization over multiple imputeted ordered variables) is to slow for my data set.

An alternative approach that seems to work is to rescale the original data to lie between 0 and 1, apply a logit transform to the rescaled data, and to impute data using a multivariate normal distribution (for the rescaled and logit-transformed data). This approach is reasonably faithful to the (my) data and orders of magnitude faster than the original approach.

On Wednesday, December 14, 2016 at 9:46:32 PM UTC+1, Guido Biele wrote:

just to make sure that I understand what your are saying:

lets say we have a design matrix for a regression model x_imputed, which consists mostly of observed data, but some entries of which are parameters (i.e. imputed).

If we now focus on one row of x and want multiplicate x_imputed[n,] with regressions weights beta, are you saying that then it would be fastest to do straight forward matrix multiplication for the rows that consist only of data, and to separate data from parameters in a rows of x_imputed that contains parameters before multiplying each predictor with its weight beta?

guido

On Wed, Dec 14, 2016 at 5:20 PM Daniel Lee   wrote:

Hi Guido,

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
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/zoelNcwfFYk/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

Neat. Thanks for reporting back.

- Bob

> 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/zoelNcwfFYk/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.