rhat and ess updates in analyze

[avehtari]

Summary:

In analyze service

• Offer interfaces a unified calculation of rhat
• Update effective sample size calculation
• Updating ess requires a small change in the autocovariance computation

Description:

#2569 and #2575 moved effective sample size computation from chains.hpp to analyze service.
Move rhat calculation in the same way. In addition make minor changes to the effective sample size and rhat calculation. These changes make the computations more robust, but in most cases affect the computed values only a little. Especially these changes don't include rank normalization, folding or quantiles (these features will be added later).

The reference implementation in R is at https://github.com/stan-dev/rstan/blob/develop/rstan/rstan/R/monitor.R

• autocovariance: use "biased" estimate as recommended by Geyer (1992)
• rhat: for constant finite values return 1, for non-finite values return NaN
• ess: 1) use split chains, 2) add half rho of the next even lag after the truncation, 3) negative values and maximum ess are set chainsn_sampleslog10(chains*n_samples)

Reproducible Steps:

Refrence values from R implementation

upath <- system.file('unitTests', package='rstan')
f1 <- file.path(upath, 'testdata', 'blocker1.csv')  
f2 <- file.path(upath, 'testdata', 'blocker2.csv')  
fit <- read_stan_csv(c(f1,f2))
monitor(fit)
sims <- as.array(fit)
rh <- rhat(sims[,,1])
checkEquals(rh, 1.007794, tolerance = 0.001); 
eb <- ess_bulk(sims[,,1])
checkEquals(eb, 465.3917, tolerance = 0.001); 
et <- ess_tail(sims[,,1])
checkEquals(et, 349.1353, tolerance = 0.001); 
em <- ess_mean(sims[,,1])
checkEquals(em, 467.8447, tolerance = 0.001); 

Same blocker1.csv and blocker2.csv files are in stan repo at src/test/unit/mcmc/test_csv_files/ and have been used in unit tests in stan repo, too.

EDIT: updated ess_mean reference value.

[avehtari]

ping @roualdes

[roualdes]

Thanks, @avehtari. I think autocovariance is in stan-dev/math. I'll start this week to work on these.

[avehtari]

I think autocovariance is in stan-dev/math

If it's used in elsewhere in math or exposed in Stan language, then it's better to re-implement in analyze

[betanalpha]

We had some discussion about this a long, long time ago. On one hand people argued that the math library should be just for the functions exposed to the autodiff library and on the other people argued that it should include _all_ mathematical functionality used in Stan. That latter side won, which is why the autocorrelation and even the Welford accumulators used in the adaptation are in the math library. That said, pushing everything into math does increase the testing and maintenance burden, where as having things contained within analyze would help distribute the burden and localize the dependences. I’d be okay going in that direction if there are no strong objections.

…

On Apr 21, 2019, at 9:51 PM, Aki Vehtari ***@***.***> wrote: I think autocovariance is in stan-dev/math If it's used in elsewhere in math or exposed in Stan language, then it's better to re-implement in analyze — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#2752 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AALU3FUUX7D4BVD6RXM6YYLPRUKY3ANCNFSM4HHKZHUA>.

[roualdes]

The reference values from the R implementation was a great idea, thank you. Just one weird issue about it. When I run the following R code,

library(rstan)
source('monitor.R')
upath <- system.file('unitTests', package='rstan')
f1 <- file.path(upath, 'testdata', 'blocker1.csv')
f2 <- file.path(upath, 'testdata', 'blocker2.csv')
fit <- read_stan_csv(c(f1,f2))
sims <- as.array(fit)
ess_mean(sims[,,1])

I get the result 467.8447. I also get this value in C++ based on my best efforts to update the function stan::analyze::compute_effective_sample_size(). However, this number does not match the number you provided, 467.3676. Would you please confirm or disconfirm the results of this R code? I'm I reading the number you provided incorrectly? Other? Thanks.

[avehtari]

I get the result 467.8447.

I now get the same result. My previous result was by mistake without splitting the chains, ie, the result you would get from ess_rfun(sims[,,1]). Sorry for the confusion.

[roualdes]

Thanks for double checking.

Let me know what you guys think of the following functional spec. I propose to overload stan::analyze::compte_effective_sample_size() with

double compute_effective_sample_size(std::vector<const double*> draws, std::vector<size_t> sizes)
// and
double compute_effective_sample_size(const Eigen::MatrixBase<Derived>& draws)

which attempts to accomplish two things. The first signature aligns with our discussions on Discourse. The second signature allows for composing stan::analyze::compute_effective_sample_size() with stan::analyze::split_chains(), in the same way that monitor.R does, namely compute_effective_sample_size(split_chains(...)). The arguments to stan::analyze::split_chains() follow the same logic from our Discourse discussion:

Eigen::MatrixXd split_chains(const std::vector<const double*>& draws, const std::vector<size_t>& sizes)

To enable such a separation, without needlessly copying data around, I added an (intended to be) internal function

double effective_sample_size_impl(const Eigen::MatrixBase<DerivedA>& acov, const Eigen::MatrixBase<DerivedB>& chain_mean)

I agree that this adds much complication, to what should/could be a simple PR. However, the outcome is that the interfaces can call stan::analyze::compute_effective_sample_size() and stan::analyze::split_chains() separately without being penalized for needlessly copying data around, and more importantly, this will ease the future effort to adopt ess_bulk, which necessitates the separation of stan::analyze::compute_effective_sample_size() and stan::analyze::split_chains().

What do you guys think?

[betanalpha]

Moving towards an API that facilitates chain splitting, for both ESS and Rhat calculations, is the right step forwards. I’m wondering if internally this all can’t be handled with Eigen maps? In other words, the main API exposes functions that take in a std::vector of double arrays and the first internal step is laying out the memory into Eigen maps which can then be split/passed around without any unnecessary copying. This would have the double bonus of using the memory provided by the client.

…

On Apr 26, 2019, at 4:00 PM, Edward A. Roualdes ***@***.***> wrote: Thanks for double checking. Let me know what you guys think of the following functional spec. I propose to overload stan::analyze::compte_effective_sample_size() with double compute_effective_sample_size(std::vector<const double*> draws, std::vector<size_t> sizes) // and double compute_effective_sample_size(const Eigen::MatrixBase<Derived>& draws) which attempts to accomplish two things. The first signature aligns with our discussions on Discourse <https://discourse.mc-stan.org/t/analysis-api/3486/30>. The second signature allows for composing stan::analyze::compute_effective_sample_size() with stan::analyze::split_chains(), in the same way that monitor.R <https://github.com/stan-dev/rstan/blob/24c9962a8f425f67453bb14034d5317633830fec/rstan/rstan/R/monitor.R#L378> does, namely compute_effective_sample_size(split_chains(...)). The arguments to stan::analyze::split_chains() follow the same logic from our Discourse discussion: Eigen::MatrixXd split_chains(const std::vector<const double*>& draws, const std::vector<size_t>& sizes) To enable such a separation, without needlessly copying data around, I added an (intended to be) internal function double effective_sample_size_impl(const Eigen::MatrixBase<DerivedA>& acov, const Eigen::MatrixBase<DerivedB>& chain_mean) I agree that this adds much complication, to what should/could be a simple PR. However, the outcome is that the interfaces can call stan::analyze::compute_effective_sample_size() and stan::analyze::split_chains() separately without being penalized for needlessly copying data around, and more importantly, this will ease the future effort to adopt ess_bulk <https://github.com/stan-dev/rstan/blob/24c9962a8f425f67453bb14034d5317633830fec/rstan/rstan/R/monitor.R#L307>, which necessitates the separation of stan::analyze::compute_effective_sample_size() and stan::analyze::split_chains(). What do you guys think? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2752 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AALU3FQ6BB5LFQPXU2TR6Z3PSNNN7ANCNFSM4HHKZHUA>.

[bob-carpenter]

What do you guys think?

I'm on board because I'm the one who was suggesting something like this in the first place to have a baseline implementation. We probably need to do the same thing with our probability functions to cut down on compile times.

I’m wondering if internally this all can’t be handled with Eigen maps?

Having aT* is a necessary and sufficient condition to create an Eigen::Map<Matrix<T, ...>>. One of the motivations here is the difficulty they've had with PyStan integrating Eigen data types through Cython and the desire to provide a non-Eigen interface.

This would have the double bonus of using the memory provided by the client.

There are competing issues of speed and scalability. The roadblock is that draws come in row-by-row (draw by draw, that is), but all the analysis is column-by-column (parameter by parameter). An efficient Eigen transpose involving a copy up front is probably worthwhile versus using non-memory local Maps everywhere. So it should be faster to transpose once, but it requires twice the dynamic memory overhead as now you have the original memory and a copy.

[roualdes]

I’m wondering if internally this all can’t be handled with Eigen maps?

I thought about and tried this, but I convinced myself that it won't work. An Eigen map is created with a pointer to a contiguous block of memory and the dimension(s). The trouble is, as far as I know, std::vector<const double*>& draws will not necessarily keep all arrays next to each other, despite the fact that the elements of each double* are contiguous.

This also means split_chains() necessarily does a copy.

If you've got 5 - 20 minutes to educate me about Eigen maps or convince yourself that I'm on to something, please do.

Here's my a link to my work-in-progress branch, which is on my fork of stan.

[betanalpha]

I wasn’t thinking of merging everything into a matrix via a map but rather handling the initial chains, or any splitting there of, as a collection of vector maps. I guess it all comes down to how the autocorrelations are handled. Welford accumulator estimators can run on each chain segment and hence don’t need everything to be in one big contiguous block of memory. Plus the individual estimators can be computed in parallel. On the other hand FFT estimators, as are currently employed, probably need the continuity. Anyways, to maintain flexibility moving forwards I think we want to have something that flows like input chain segments -> splitting into more chain segments -> analysis code that reads in an arbitrary number of chain segments In other words any copying/aggregation into contiguous memory shouldn’t happen until the analysis code. In that way we can rewrite the underlying analysis code to avoid copying in the future if desired.

…

On Apr 27, 2019, at 1:10 PM, Edward A. Roualdes ***@***.***> wrote: I’m wondering if internally this all can’t be handled with Eigen maps? I thought about and tried this, but I convinced myself that it won't work. An Eigen map is created with a pointer to a contiguous block of memory and the dimension(s). The trouble is, as far as I know, std::vector<const double*>& draws will not necessarily keep all arrays next to each other, despite the fact that the elements of each double* are contiguous. This also means split_chains() necessarily does a copy. If you've got 5 - 20 minutes to educate me about Eigen maps or convince yourself that I'm on to something, please do. Here's my a link to my work-in-progress branch <develop...roualdes:feature/issue-2752-update-rhat-ess>, which is on my fork of stan. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2752 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AALU3FROB3C67NZTNMSZVADPSSCILANCNFSM4HHKZHUA>.

[bob-carpenter]

On Apr 28, 2019, at 1:52 PM, Michael Betancourt ***@***.***> wrote: I wasn’t thinking of merging everything into a matrix via a map but rather handling the initial chains, or any splitting there of, as a collection of vector maps.

That seems fine. It's the locality of the individual chains that really matters here as that's the unit at which we use the data.

I guess it all comes down to how the autocorrelations are handled. Welford accumulator estimators can run on each chain segment and hence don’t need everything to be in one big contiguous block of memory.

Right. It's just that reading out of contiguous memory is faster.

Plus the individual estimators can be computed in parallel. On the other hand FFT estimators, as are currently employed, probably need the continuity.

They assume an Eigen structure as input---it may not need to be literally stored as a dense matrix, but it'll probably do an internal copy if it's not for efficiency. So if FFT's already going to do that, we might as well do it ourselves to save it the trouble and let us use it elsewhere.

Anyways, to maintain flexibility moving forwards I think we want to have something that flows like input chain segments -> splitting into more chain segments -> analysis code that reads in an arbitrary number of chain segments In other words any copying/aggregation into contiguous memory shouldn’t happen until the analysis code. In that way we can rewrite the underlying analysis code to avoid copying in the future if desired.

The copying/aggreation into memory is happening in RStan, for example, on an iteration-by-iteration basis. In CmdStan, we just write to disk, but it's still local iteration-by-iteration. I think reading all the data in and transposing then calculating R-hat, ESS, etc. will be more efficient. It's just that it's also more memory intensive.

[roualdes]

I wasn’t thinking of merging everything into a matrix via a map
but rather handling the initial chains, or any splitting there of,
as a collection of vector maps.

They way I've coded it, if you don't call stan::analyze::split_chains(), then it uses vector maps and the memory provided by the client. On the other hand, I couldn't figure out how to reshape a map. So anytime stan::analyze::split_chains() is called, something closer to

reading all the data in and transposing then calculating R-hat, ESS, etc

is happening.

As it's my understanding that I'm quite close to respecting the discussion points so far, I'm going to push on. I'm going to finish up rhat, some documentation, and a couple of other things before I submit a PR.

Thank you all for your time.

[betanalpha]

They way I've coded it, if you don't call stan::analyze::split_chains(), then it uses vector maps and the memory provided by the client. On the other hand, I couldn't figure out how to reshape a map.

I’m not sure why you’d have to deal with any reshaping — splitting the chains can be accomplished by creating a new collection of std::vectors that point to segments of the input std::vector memory.

So anytime stan::analyze::split_chains() is called, something closer to reading all the data in and transposing then calculating R-hat, ESS, etc is happening

This worries me a bit — the implementation of analyze should be independent of any possible intermediate splitting. If split chains might not be contiguous, and require require copying into contiguous memory, then the initial chains might also not be contiguous and require the copy. In other words, analyze shouldn’t know whether the chains where split or not. It just consumes an array of chains segments and runs the diagnostics ignorant of where those segments came from. In any case this can be addressed in the PR if you prefer.

[roualdes]

splitting the chains can be accomplished by creating a new collection of std::vectors
that point to segments of the input std::vector memory

That's a reasonable way to split chains in order to use the memory provided by the client. But I still don't see how to use the memory provided by the client and separate split_chains() from compute_effective_sample_size(). Even if split_chains() created a new collection of std::vectors that point to segments of the input std::vector memory, then you'd still need to return a second std::vector that records the sizes of the new (halved) chains -- or put call split_chains() within the scope of compute_effective_sample_size(), which I'm trying to avoid.

I guess one could have split_chains() return an Eigen map of a new collection of std::vectors that point to segments of the input std::vector memory. But it's my understanding that a copy would be incurred as soon as an Eigen Map is returned from a function. If somebody knows something about the interplay between functions that return an Eigen map and return value optimization, please share.

Without forcing split_chains() to have multiple return values (std::pair or something), compute_effective_sample_size() still needs both the number of chains and the length of each chain; its arguments are const std::vector<const double*>& draws, const std::vector<size_t>& sizes.

I keep trying to do the calculations on Eigen data structures since an Eigen::MatrixXd carries along both the number of rows (iterations) and columns (number of chains). But reshaping an Eigen Map isn't obvious to me.

This worries me a bit — the implementation of analyze should be independent of any possible intermediate splitting.

That's fair, and I hope it was only due to my lack of clarity. Neither of the two functions I'm proposing to add, Eigen::MatrixXd split_chains(const std::vector<const double*>& draws, const std::vector<size_t>& sizes) nor double compute_effective_sample_size(const Eigen::MatrixBase<Derived>& draws) require knowledge of whether or not its arguments were/are split.

Under my proposal, if an interface wants to first split chains and then calculate ESS all in C++, they would call compute_effective_sample_size(split_chains(...)). As it stands, split_chains() copies into an Eigen::MatrixXd of double the number of chains and half the number of iterations. As far as I know, Eigen matrices are in column-major order. Thus, if the user calls split_chains() first, then all calculations would be performed on contiguous memory with respect to each chain. It's my understanding that this is essentially the transpose that Bob was speaking about.

On the other hand, if the user does not call split_chains() first, then the candidate function is (the already available) double compute_effective_sample_size(std::vector<const double*> draws, std::vector<size_t> sizes). This function, too, does not require knowledge of the number of chains nor the number of iterations whether or not the chains were/are split, but does directly use the memory provided by the client via Eigen vector maps.