command options & adaptation for Metropolis

[Bob Carpenter]

Peter's got the basic Metropolis working end to end, so
it's time to think about adaptation and configuration.

1. How do we want Metropolis to be configurable?
At very least, I'm imagining:

-- adapted single step size * unit covariance

-- adapted diagonal covariance matrix

-- adapted full covariance matrix

Do we need to include user-specified step sizes and user-specified
covar matrices in addition?

If so, any suggestions on how to set up the classes?
Specifically, do we have three classes,

metropolis_unit_covar
metropolis_diag_covar
metropolis_full_covar

that let us configure adaptation vs. user input, or six
classes, the above plus user-specified matrix versions?

2. How should we be doing adaptation for Metropolis?

Does it make sense to independently adapt a global
scale (step size in HMC) and a standardized
covariance matrix (inverse mass matrix in HMC)?
We can't just let both be free due to identifiability
issues, right?

Here's what I'm imagining:

- parameters: scale epsilon, standardized covariance Sigma_std;
proposal covariance matrix: Sigma = epsilon * Sigma_std

- initialize covariance matrix Sigma to unit diagonal

- find a scale epsilon that gives us roughly 23.4% acceptance rates
[or make acceptance rate target configurable]

- run warmup iterations after 100 (200, 400, ...) iterations:
+ use last half of iterations' acceptance rate to adjust
epsilon up or down
+ use last half of iterations to estimate standardized
covariance matrix (???)

How do we standardize the covariance matrix estimate?

- Bob

[Marcus Brubaker]

It depends on how fancy we want things to get. I think a basic
version which estimates a covariance matrix given the previous samples
(full, diagonal or scaled identity) and then uses a mixture proposal
which uses the rescaled (see
http://probability.ca/jeff/ftpdir/galinart.pdf) estimated covariance
most of the time and sometimes uses a conservative proposal to avoid
excess sticking would be good. Other possibilities are to use a
mixtures of multivariate Student T distributions which takes into
account the uncertainty in estimating the covariance matrix. User
visible parameters would then be things like probability of using
adapted covariance, whether to use student t or gaussian proposals,
etc.

Also, there isn't really a stepsize in Metropolis, that's just an HMC
thing, so no, I don't think it really makes sense to adapt a global
scale separately.

Cheers,
Marcus

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

[Andrew Gelman]

No point in getting too fancy with Metropolis. After all, we have HMC. The one place Metropolis would be useful, I suppose, is with probability densities whose gradients can't easily be calculated (for example, in a toxiology model where the expected values of data given parameters are spit out from a black-box differential-equation solver).
A

[Michael Betancourt]

> Peter's got the basic Metropolis working end to end, so
> it's time to think about adaptation and configuration.
>
> 1. How do we want Metropolis to be configurable?
> At very least, I'm imagining:
>
> -- adapted single step size * unit covariance
>
> -- adapted diagonal covariance matrix
>
> -- adapted full covariance matrix
>
> Do we need to include user-specified step sizes and user-specified
> covar matrices in addition?
>
> If so, any suggestions on how to set up the classes?
> Specifically, do we have three classes,
>
> metropolis_unit_covar
> metropolis_diag_covar
> metropolis_full_covar
>
> that let us configure adaptation vs. user input, or six
> classes, the above plus user-specified matrix versions?

I was picturing a complete analogy to the HMC code,

metro_unit_covar
metro_diag_covar
metro_dense_covar

and then use multiple inheritance with the above and the adaptive sampler to yield

adapt_metro_unit_covar
adapt_metro_diag_covar
adapt_metro_dense_covar

> 2. How should we be doing adaptation for Metropolis?
>
> Does it make sense to independently adapt a global
> scale (step size in HMC) and a standardized
> covariance matrix (inverse mass matrix in HMC)?
> We can't just let both be free due to identifiability
> issues, right?
>
> Here's what I'm imagining:
>
> - parameters: scale epsilon, standardized covariance Sigma_std;
> proposal covariance matrix: Sigma = epsilon * Sigma_std
>
> - initialize covariance matrix Sigma to unit diagonal
>
> - find a scale epsilon that gives us roughly 23.4% acceptance rates
> [or make acceptance rate target configurable]
>
> - run warmup iterations after 100 (200, 400, ...) iterations:
> + use last half of iterations' acceptance rate to adjust
> epsilon up or down
> + use last half of iterations to estimate standardized
> covariance matrix (???)
>
> How do we standardize the covariance matrix estimate?

I don't think there will be identifiability issues here. You have d(d -1) / 2 correlations for estimating the covariance and the Metropolis acceptance rate for estimating the step size. Remember that the optimal scaling in Metropolis is not the covariance but some scaling of the covariance that depends on dimension (at least in IID gaussian models, etc). Having a free step size will avoid our having to determine that a priori.

Again, I think we can reuse all the HMC code (which is why the reusable parts are in mcmc and not mcmc/hmc!).

1) Initialize covariance matrix to unit diagonal - the regularized estimator takes care of this for us.

2) Initialize stepsize to about 0.234 - 0.5 should be close enough for an initial guess. Which means we can copy the init_stepsize code in HMC only using propose instead of evolve.

3) We should be able to adapt the step size as we go and then update the covariance after the usual windows.

[Bob Carpenter]

On 4/25/13 4:27 PM, Michael Betancourt wrote:

> Again, I think we can reuse all the HMC code (which is why the r
> reusable parts are in mcmc and not mcmc/hmc!).

Sounds good. Reuse base code for MCMC in both HMC and in Metropolis
and whatever else we do.

> 1) Initialize covariance matrix to unit diagonal - the regularized estimator takes care of this for us.

OK.

> 2) Initialize stepsize to about 0.234 - 0.5 should be close enough for an initial guess.
> Which means we can copy the init_stepsize code in HMC only using propose instead of evolve.

Why do we have to use init_stepsize() code in HMC if
we're just initializing to a fixed value?

And shouldn't that code go in stan::mcmc if we're going
to reuse it in Metropolis? Or did you just mean
cut-and-paste it with modifications?

We used to have a find_reasonable_stepsize() method that
Matt wrote in the constructor for HMC/NUTS that tried to
find a step size with a decent acceptance rate before
starting. Is this what init_stepsize() does? And if so,
why does it keep getting called during adaptation?

> 3) We should be able to adapt the step size as we go and then update the covariance after the usual windows.

How do we adapt the step size as we go? It looks like what
it's doing now is calling the equivalent of find_reasonable_stepsize()
at each iteration, which multiplies step size by 2 or cuts it
in half repeatedly based on proposals (from an initial position or
what -- here's where the code loses me?).

(Aside question: wouldn't a finer-grained factor like
like sqrt(2) or sqrt(sqrt(2)) be better? Adaptation is slower
by a factor of 2 or 4, but we'd get closer
than to within a factor of 2 of the "right" answer.)

I think I understand the covar matrix estimation, but let
me make sure. At intervals during warmup (say after 100,
after 200, after 400 iterations, etc.), we'll take the
second half of the iterations (so after 100, we'd take
iterations 51--100, after 200 we take iterations 101-200) and
use the "sampled" parameter values to estimate a covariance matrix.
We can do this in a fixed amount of memory and reasonable numerical
stability with Welford's algorithm.

Until we get to the end. If someone specifieds 875 warmup
iterations, then what do we do? We'll have buffered
400--800, then we'll start again, but don't just want to use
the last 75. And this'll be a huge mess in the time-based
approach when we don't even know how many warmup iterations
we'll get.

Also, someone mentioned only using the warmup iterations that
were not Metropolis rejected. Should we do that?

- Bob

[Bob Carpenter]

The main reason to have Metropolis, I think, is for
comparison. So we need a "state of the art" version of
Metropolis for that.

- Bob

[Michael Betancourt]

> Why do we have to use init_stepsize() code in HMC if
> we're just initializing to a fixed value?
>
> And shouldn't that code go in stan::mcmc if we're going
> to reuse it in Metropolis? Or did you just mean
> cut-and-paste it with modifications?

> We used to have a find_reasonable_stepsize() method that
> Matt wrote in the constructor for HMC/NUTS that tried to
> find a step size with a decent acceptance rate before
> starting. Is this what init_stepsize() does? And if so,
> why does it keep getting called during adaptation?

I meant copy and paste with modifications. The version in base_hmc is just Matt's algorithm tweaked slightly to use the proper (q, p) acceptance probability instead of just a (q) acceptance probability (which is really undefined since it's missing a volume correction).

> > 3) We should be able to adapt the step size as we go and then update the covariance after the usual windows.
>
> How do we adapt the step size as we go? It looks like what
> it's doing now is calling the equivalent of find_reasonable_stepsize()
> at each iteration, which multiplies step size by 2 or cuts it
> in half repeatedly based on proposals (from an initial position or
> what -- here's where the code loses me?).

Init_stepsize is called only after --the (co)variance is updated--, not at every iteration. The problem is that when you update the (co)variance the optimal step size changes, and initializing the step size again avoids unnecessary burden on the dual-averaging (remember we're not using a "standardized" metric anymore).

> (Aside question: wouldn't a finer-grained factor like
> like sqrt(2) or sqrt(sqrt(2)) be better? Adaptation is slower
> by a factor of 2 or 4, but we'd get closer
> than to within a factor of 2 of the "right" answer.)

Sure, I didn't consider any improvements but that would be easy.

> I think I understand the covar matrix estimation, but let
> me make sure. At intervals during warmup (say after 100,
> after 200, after 400 iterations, etc.), we'll take the
> second half of the iterations (so after 100, we'd take
> iterations 51--100, after 200 we take iterations 101-200) and
> use the "sampled" parameter values to estimate a covariance matrix.
> We can do this in a fixed amount of memory and reasonable numerical
> stability with Welford's algorithm.

Yes.

> Until we get to the end. If someone specifieds 875 warmup
> iterations, then what do we do? We'll have buffered
> 400--800, then we'll start again, but don't just want to use
> the last 75. And this'll be a huge mess in the time-based
> approach when we don't even know how many warmup iterations
> we'll get.

Yeah. In fact, I think I've seen behavior to this effect already (it shows up as unoptimal acceptance rates in the final samples).
Should we not make the windows arbitrary but scale them to the total number of warm-up iterations? i.e. assume five windows and just take the first window to be of size (warm-up iterations) / 2^{5}, second window (warm-up) iterations) / 2^{4}, etc.

> Also, someone mentioned only using the warmup iterations that
> were not Metropolis rejected. Should we do that?

Only concern would be getting stuck in local minima.

[Peter Li]

Well, I got adapt_unit_metro and adapt_diag_metro working under Michael's framework, or at least I think they're working (no actual testing done yet). But when I try and compile with adapt_dense_metro, I get this output

dyn-160-39-0-93:stan Peter$ make models/basic_distributions/binormal CC=clang++ O=0 -j2
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -Wno-unused-function -Wno-tautological-compare -c -O0 -o models/basic_distributions/binormal.o models/basic_distributions/binormal.cpp
clang++ -I src -I lib/eigen_3.1.2 -I lib/boost_1.53.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -Wno-unused-function -Wno-tautological-compare -O0 -o models/basic_distributions/binormal models/basic_distributions/binormal.o -Lbin -lstan
dyn-160-39-0-93:stan Peter$ ./models/basic_distributions/binormal --unit_metroSTAN SAMPLING COMMAND
data = (specified model requires no data)
init = random initialization
init tries = 1
samples = samples.csv
append_samples = 0
save_warmup = 0
seed = 1258117860 (randomly generated)
chain_id = 1 (default)
iter = 2000
warmup = 1000
thin = 1
equal_step_sizes = 0
nondiag_mass = 0
leapfrog_steps = -1
max_treedepth = 10
epsilon = -1
epsilon_pm = 0
delta = 0.5
gamma = 0.05
run metropolis
Iteration: 1 / 2000 [ 0%] (Sampling)
Assertion failed: (other.rows() == 1 || other.cols() == 1), function resizeLike, file lib/eigen_3.1.2/Eigen/src/Core/PlainObjectBase.h, line 306.
ohiheythurAbort trap: 6

I'm not sure what that last error msg is referring to. I tried using print statements to isolate what was causing the error, but I couldn't find it. I just know that it was able to instantiate adapt_dense_metro.


Thanks,
Peter

The code is in the metropolis branch on github under mcmc/metro.

[Peter Li]

On Apr 25, 2013, at 8:59 PM, Peter Li wrote:

> Well, I got adapt_unit_metro and adapt_diag_metro working under Michael's framework, or at least I think they're working (no actual testing done yet). But when I try and compile with adapt_dense_metro, I get this output

It compiles fine. I meant that when I try to run the model, I get the following error/output:

[Ben Goodrich]

On Thursday, April 25, 2013 9:03:22 PM UTC-4, Peter Li wrote:

On Apr 25, 2013, at 8:59 PM, Peter Li wrote:

> Well, I got adapt_unit_metro and adapt_diag_metro working under Michael's framework, or at least I think they're working (no actual testing done yet). But when I try and compile with adapt_dense_metro, I get this output

It compiles fine. I meant that when I try to run the model, I get the following error/output:

Runs here (59251b0).

Ben

[Michael Betancourt]

It also works for me.

Firstly, watch out out for the camelCase -- Bob's going to find you! Move to prop_covar, etc.

Secondly, you've got some nastiness in your constructors

adapt_dense_metro(M &m, BaseRNG& rng, std::ostream* error_msg)
: dense_metro<M, BaseRNG>(m, rng, this->_propCov, error_msg),
stepsize_covar_adapter(m.num_params_r())
{
this->_propCov.resize(m.num_params_r(), m.num_params_r());
this->_propCov.setIdentity();
};

When the dense_metro constructor is called this->_propCov hasn't been constructed yet: you have a circular constructor! Personally I wouldn't even pass var/covar as an argument (the Welford adaptation is pretty quick) but just a dimension to be used for constructing _propCov, especially if you're going to hide the option in the adapt_ classes!

[Peter Li]

Alright I fixed those, but I still don't understand what's going wrong. Are you guys calling adapt_dense_metro with something like the following?

dyn-209-2-236-117:stan Peter$ ./models/basic_distributions/binormal --dense_metro


Thanks,
Peter

[Michael Betancourt]

I was.

Here's the issue in detail. You've defined the constructor as

adapt_dense_metro(M &m, BaseRNG& rng, std::ostream* error_msg)
: dense_metro<M, BaseRNG>(m, rng, this->_propCov, error_msg),
stepsize_covar_adapter(m.num_params_r())
{
this->_propCov.resize(m.num_params_r(), m.num_params_r());
this->_propCov.setIdentity();
};

When you try to construct adapt_dense_metro multiple constructors are called, from base to derived. So in this case you call

base_mcmc()
base_metro()

dense_metro<M, BaseRNG>(m, rng, this->_propCov, error_msg)

base_adapter()
stepsize_covar_adapter(m.num_params_r())

and then finally

{
this->_propCov.resize(m.num_params_r(), m.num_params_r());
this->_propCov.setIdentity();
};

Let's focus on the dense_metro constructor,

dense_metro<M, BaseRNG>(m, rng, this->_propCov, error_msg)

You've given this->_propCov as an argument, but it doesn't exist yet! It only gets created once the arguments have been passed and the constructor executes. It's like writing

double x = x;

I'm not sure what the spec says, but I wouldn't be surprised if it gives undefined behavior which would explain why Ben and I see different results than you.

-Mike

[Peter Li]

I pulled prop_cov out (all on github). Now, it's just this

adapt_dense_metro(M &m, BaseRNG& rng, std::ostream* error_msg)

: dense_metro<M, BaseRNG>(m, rng, error_msg),
stepsize_covar_adapter(m.num_params_r())
{
this->_prop_cov.resize(m.num_params_r(), m.num_params_r());
this->_prop_cov.setIdentity();
};

which is still giving me the same error.

Originally, I had included it so users could call just dense_metro, providing their own covariance matrix.


Peter

[Michael Betancourt]

Be sure to initialize the covariance matrix in the dense_metro constructor,

dense_metro(M& m,

BaseRNG& rng,
std::ostream* error_msg)

: base_metro<M, BaseRNG>(m, rng, error_msg),
_prop_covar(Eigen::MatrixXd::Identity(m.num_params_r(), m.num_params_r()) {
this->_name = "Metropolis with a dense Euclidean metric";
this->_nom_epsilon = 1;
}

You don't need to mess with _prop_covar in the adapt_dense_metro constructor.

Eventually we should gives users the option to input their own covariance matrix, but I think that's more appropriate as a seed or some function that reads in a file and builds the matrix. In any case, let's get the basic version working first.

[Peter Li]

Done. Adaptive unit/diag/dense metro should all be working now (and are in the metropolis branch).

The error before was that I was passing a square matrix in for the mu argument, which should've only been a vector.

[Bob Carpenter]

Super cool. Now we need to get all the doc in.

But let's wait until the dust settles on Michael's
refactoring of the commands to do that.

- Bob

[Michael Betancourt]

Almost done with the refactoring -- the code works but I'm just trying to figure out the best way to minimize code duplication in setting all of the sampler parameters (turns out the usual OO tricks interact terribly with templated classes :-) ). It would also be easier to do the updated I/O and diagnostics before making the argument changes, if we're not in any hurry.

[Bob Carpenter]

Just about everything seems to interact badly with
templates.

The last error in the model that wouldn't compile is because
of namespace and template specialization interactions.
It tried to cast an int to a var rather than to a double.
The fix of course is more templating :-)

This should all be done by end of tomorrow, but I'm going to
have to go back through every function that specifies
a double and template it out so it matches int better
than var matches int.

- Bob

[Michael Betancourt]

> Until we get to the end. If someone specifieds 875 warmup
> iterations, then what do we do? We'll have buffered
> 400--800, then we'll start again, but don't just want to use
> the last 75. And this'll be a huge mess in the time-based
> approach when we don't even know how many warmup iterations
> we'll get.

To be precise, any extra samples wouldn't have affected anything because the (co)variance is updated only when a window closes, not when it opens. In any case, I updated the adapters to take in the total number of adaptation iterations (just the number of warm up samples in our case) and then snap the final adaptation window to the total warmup window so that no samples are wasted.

It should be straightforward to make this work with time-based approaches, as well, if we can assume constant samples/time.

[Jiqiang Guo]

On Sun, Apr 28, 2013 at 12:47 PM, Michael Betancourt <betan...@gmail.com> wrote:

It should be straightforward to make this work with time-based approaches, as well, if we can assume constant samples/time.

Could you elaborate this a bit more?  I am working on implementing the time-bases approaches in rstan.  I do not know how to do it now as the total number of warmup iterations cannot be determined a priori. 

--

Jiqiang 

[Michael Betancourt]

The idea would be to estimate samples/time in vivo. So instead of sampling until you hit the max time, like

sample
if(time > max_time) break

you would time each sample, use that to estimate the number of samples remaining in the given time window, and then set the total number of samples to equal to the samples already run plus those estimated to be remaining.

start clock
sample
stop clock
sample_per_time = (n * sample_per_time + 1 / (delta t) ) / ( n + 1)
samples_remaining = remaining_time * sample_per_time
n_samples = n + samples_remaining

The final timing would be less accurate, but you get to have both features. There might be easier ways to accomplish this.

[Bob Carpenter]

What's the issue with not knowing number of
warmup iterations?

It'd be nice to have time per iteration estimators
even for the fixed number of iterations so we
can estimate time remaining. I think it'd be nice
to be able to see iteration number, percentage
iterations complete, log probability, and
estimated time remaining.

As far as estimating time per iteration,
it might vary pretty widely during warmup, so
I'm not sure how stable it will be and when we
can do this. Especially in the regime Andrew
imagines where you intentially start out at a point
that's likely to be too little time to converge.

Even after convergence, there can be a lot of
variation in tree depth explored (and hence time
taken per iteration).

- Bob

[Michael Betancourt]

> What's the issue with not knowing number of
> warmup iterations?

Then you don't know when to stop the adaptation, and you get the adaptation window misalignment we had discussed earlier.

> It'd be nice to have time per iteration estimators
> even for the fixed number of iterations so we
> can estimate time remaining. I think it'd be nice
> to be able to see iteration number, percentage
> iterations complete, log probability, and
> estimated time remaining.

Wouldn't be too hard -- add to the to-do?

> As far as estimating time per iteration,
> it might vary pretty widely during warmup, so
> I'm not sure how stable it will be and when we
> can do this. Especially in the regime Andrew
> imagines where you intentially start out at a point
> that's likely to be too little time to converge.
>
> Even after convergence, there can be a lot of
> variation in tree depth explored (and hence time
> taken per iteration).

Sure, I don't disagree. But the total time to convergence can also vary dramatically, which is at odds with sampling for only a preset amount of time. Just thinking of ideas to try and do as best as possible given somewhat dissonant requirements.

[Bob Carpenter]

On 4/28/13 3:05 PM, Michael Betancourt wrote:
>> What's the issue with not knowing number of
>> warmup iterations?
>
> Then you don't know when to stop the adaptation, and you get the adaptation window misalignment we had discussed earlier.

Can we do something like a moving average type adaptation where
old information dies out in strength over time? Then there'd be
no discretization of the epochs for training covariance during
warmup.

I'd think it'd be easy to add to Welford's algorithm by adjusting
the ratio in the updates to do some past discounting.

- Bob

[Michael Betancourt]

>
> Can we do something like a moving average type adaptation where
> old information dies out in strength over time? Then there'd be
> no discretization of the epochs for training covariance during
> warmup.
>
> I'd think it'd be easy to add to Welford's algorithm by adjusting
> the ratio in the updates to do some past discounting.

Technically that's not hard to do. I think the problem with a moving average is that it can still be too sensitive to the past. In particular, convergence of the Markov chain (and, perhaps more importantly, the disparity in values between unconverted and converged) can often outpace the forgetfulness of the moving average estimator. Isn't that why we went to sequential windows in the first place?