From 8b73342a11976b89d2210cbf15d5c7dc66c03f3f Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 23 Mar 2020 12:31:38 -0400 Subject: [PATCH 1/7] Added tutorial for reduce_sum (design-docs pull request #17) --- reduce_sum_tutorial.Rmd | 284 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 284 insertions(+) create mode 100644 reduce_sum_tutorial.Rmd diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd new file mode 100644 index 0000000..8e51093 --- /dev/null +++ b/reduce_sum_tutorial.Rmd @@ -0,0 +1,284 @@ +--- +title: "Reduce Sum: A Minimal Example" +date: "23 Mar 2020" +output: html_document +--- + +```{r setup, include=FALSE} +library("rstan") +``` + +This introduction to `reduce_sum` copies directly from Richard McElreath's [Multithreading and Map-Reduce in Stan 2.18.0: A Minimal Example](https://github.com/rmcelreath/cmdstan_map_rect_tutorial). + +## Introduction + +Stan 2.23 introduced `reduce_sum`, a new way to parallelize the execution of a single Stan chain across +multiple cores. This is in addition to the already existing `map_rect` utility, and introduces a number of +features that make it easier to use parallelism. + +The main advantages to `reduce_sum` are: + +1. `reduce_sum` has a more flexible argument interface, avoiding the packing and unpacking that was necessary with `map_rect`. +2. `reduce_sum` partitions the data for parallelization automatically. + +The drawbacks to `reduce_sum` are: + +1. `reduce_sum` requires that the result of the calculation be a scalar, while `map_rect` returns a list of vectors. +2. `reduce_sum` is only parallelized across multiple cores, not across multiple computers like `map_rect`. + +This tutorial will highlight these advantages on a simple logistic regression using `cmdstan`. + +## Overview + +Modifying a Stan model for use with `reduce_sum` has three steps. + +(1) Identify the portion of the model that you want to parallelize. `reduce_sum` is design specifically for speeding up log likelihood evaluations that are composed of a number of conditionally independent terms that can be computed in parallel then added together. + +(2) Write a reduce function that can compute a given number of terms that must be added together. Pass this function (and other necessary arguments) to a call to `reduce_sum`. + +(3) Configure your compiler to enable multithreading. Then you can compile and sample as usual. + +In the rest of this tutorial, we'll first build the model frist without and then with `reduce_sum`. + +[EDIT: Does Windows work with threading yet?] + +If you use Windows, it looks like multithreading isn't working yet (see ). So flip over to a Linux partition, if you have one. If you don't, hope is that an upcoming update to RTools will make things work. + +## Prepping the data + +Let's consider a simple logistic regression, just so the model doesn't get in the way. We'll use a reasonably big data table, the football red card data set from the recent crowdsourced data analysis project (). This data table contains 146,028 player-referee dyads. For each dyad, the table records the total number of red cards the referee assigned to the player over the observed number of games. + +The ``RedcardData.csv`` file is provided in the repository here. Load the data in R, and let's take a look at the distribution of red cards: + ```R +d <- read.csv( "RedcardData.csv" , stringsAsFactors=FALSE ) +table( d$redCards ) +``` +```text +0 1 2 +144219 1784 25 +``` +The vast majority of dyads have zero red cards. Only 25 dyads show 2 red cards. These counts are our inference target. + +The motivating hypothesis behind these data is that referees are biased against darker skinned players. So we're going to try to predict these counts using the skin color ratings of each player. Not all players actually received skin color ratings in these data, so let's reduce down to dyads with ratings: + ```R +d2 <- d[ !is.na(d$rater1) , ] +out_data <- list( n_redcards=d2$redCards , n_games=d2$games , rating=d2$rater1 ) +out_data$N <- nrow(d2) +``` +This leaves us with 124,621 dyads to predict. + +At this point, you are thinking: "But there are repeat observations on players and referees! You need some cluster variables in there in order to build a proper multilevel model!" You are right. But let's start simple. Keep your partial pooling on idle for the moment. + +Now to use these data with ``cmdstan``, we need to export the table to a file that ``cmdstan`` can read in. This is made easy: +```R +library(rstan) +stan_rdump(ls(out_data), "redcard_input.R", envir = list2env(out_data)) +``` + +## Making the model + +A Stan model for this problem is just a simple logistic (binomial) GLM. I'll assume you know Stan well enough already that I can just plop the code down here. It's contained in the ``logistic0.stan`` file in the repository. It's not a complex model: + +``` + data { + int N; + int n_redcards[N]; + int n_games[N]; + vector[N] rating; + } +parameters { + vector[2] beta; +} +model { + beta ~ normal(0,1); + n_redcards ~ binomial_logit(n_games , beta[1] + beta[2] * rating); +} +``` + +## Building and Running the Basic Model + +Make a new folder called `redcard` inside your cmdstan folder folder and drop the `redcard_input.R` and `logistic0.stan` files in it. To build the Stan model, while still in the cmdstan folder, enter: + +```bash +make redcards/logistic0 +``` + +To sample from the model: + +```bash +cd redcard +time ./logistic0 sample data file=redcard_input.R +``` + +The `time` command at displays a a processor time report at the end. We'll want to compare this time to the time you get later, after you enable multithreading. On my computer, I get: + +``` +real 2m11.467s +user 1m55.844s +sys 0m15.553s +``` + +The first line is the "real" time, the one you probably care about for benchmarking. + +By default, the output file is called `output.csv`. You'll want to go back into R to analyze it. Just open R in the same folder and enter: + +```R +library(rstan) +m0 <- read_stan_csv("output.csv") +``` + +Now you can proceed as usual to work with the samples. + +## Rewriting the Model to Enable Multithreading + +The key point to getting this calculation into `reduce_sum`, is recognizing that +the statement: + +``` +n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating); +``` + +can be rewritten (up to a proportionality constant) as: +``` +for(n in 1:N) { + target += binomial_logit_lpmf(n_redcards[n] | n_games[n], beta[1] + beta[2] * rating[n]) +} +``` + +Now it is clear that the calculation is the sum of a number of conditionally +independent Bernoulli log probability statements, which is the condition where +`reduce_sum` is useful. + +To use `reduce_sum`, a function must be written that can be used to compute +arbitrary subsets of the sums. + +Using the reducer interface defined in [Reduce-Sum](#reduce-sum), such a function +can be written like: + +``` +functions { + real reduce_func(int start, int end, + int[] slice_n_redcards, + int[] n_games, + vector rating, + vector beta) { + return binomial_logit_lpmf(slice_n_redcards | + n_games[start:end], + beta[1] + beta[2] * rating[start:end] ); + } +} +``` + +And the likelihood statement in the model can now be written: + +``` +target += reduce_func(1, N, n_redcards, n_games, x, beta); // Sum terms 1 to N in the likelihood +``` + +In this example, `n_redcards` was chosen to be sliced over because there +is one term in the summation per value of `n_redcards`. Technically `n_games` +or `rating` would have worked just as well. Use whatever conceptually makes +the most sense. + +Because `n_games` and `rating` are shared arguments, they must be subset +manually with `start:end`. + +With this function, `reduce_sum` can be used to automatically parallelize the +likelihood: + +``` +int grainsize = 100; +target += reduce_sum(reduce_func, n_redcards, grainsize, + n_games, rating, beta); +``` + +`reduce_sum` automatically breaks the sum into roughly `grainsize` sized pieces +and computes them in parallel. `grainsize = 1` specifies that the grainsize should +be estimated automatically. + +Making grainsize data (this makes it convenient to find a good one), the final model +should look like: +``` +functions { + real reduce_func(int start, int end, + int[] slice_n_redcards, + int[] n_games, + vector rating, + vector beta) { + return binomial_logit_lpmf(slice_n_redcards | + n_games[start:end], + beta[1] + beta[2] * rating[start:end] ); + } +} +data { + int N; + int n_redcards[N]; + int n_games[N]; + int grainsize; + vector[N] rating; +} +parameters { + vector[2] beta; +} +model { + beta ~ normal(0,1); + + target += reduce_sum(reduce_func, n_redcards, grainsize, + n_games, rating, beta); +} +``` + +Save this as `redcard/logistic1.stan`. + +## Running the Multithreaded Model + +### Configuring Cmdstan + +Now we're almost ready to go. First, make sure threading support is enabled in cmdstan (instructions [here](https://github.com/stan-dev/math/wiki/Threading-Support)). That basically means adding: + +``` +STAN_THREADS=true +``` + +to `make/local`. If this was not set, before building the model, clean your install with a: + +``` +make clean-all +``` + +and then build it again using: + +``` +make -j4 build +``` + +(setting the argument to `j` equal to the number of cores you want to build in parallel with) + +### Build and Run + +```bash +make redcards/logistic0 +``` + +To sample from the model: + +```bash +cd redcard +time STAN_NUM_THREADS=8 ./logistic0 sample data file=redcard_input.R +``` + +The environment variable `STAN_NUM_THREADS` defines how many threads to use for this calculation. On this computer I have 8-cores, and so I choose to use 8 threads. + +The time command reports: + +``` +real 0m18.086s +user 2m15.234s +sys 0m0.783s +``` + +This gives roughly a 7x speedup. This model was particularly easy to parallelize because a lot of the arguments are data (and do not need to be autodiffed). If a model has a large number of arguments that are either defined in the parameters block, transformed parameters block, or model block, or do not do very much computation inside the reduce function, the speedup will be much more limited. + +## More + +There is a new section in the 2.23 User Manual and Function Reference that contains more details of using `reduce_sum`. From 03b3686cdc9493825d5d36503ac517dfa6f474df Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 23 Mar 2020 12:34:04 -0400 Subject: [PATCH 2/7] Slight modifications of text for reduce_sum example (design-docs pull request #17) --- reduce_sum_tutorial.Rmd | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd index 8e51093..a5c6d35 100644 --- a/reduce_sum_tutorial.Rmd +++ b/reduce_sum_tutorial.Rmd @@ -232,23 +232,21 @@ Save this as `redcard/logistic1.stan`. ## Running the Multithreaded Model -### Configuring Cmdstan - Now we're almost ready to go. First, make sure threading support is enabled in cmdstan (instructions [here](https://github.com/stan-dev/math/wiki/Threading-Support)). That basically means adding: -``` +```bash STAN_THREADS=true ``` to `make/local`. If this was not set, before building the model, clean your install with a: -``` +```bash make clean-all ``` and then build it again using: -``` +```bash make -j4 build ``` From 2dc8dcb43915de746b08b7e2be09a52b76db989b Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 23 Mar 2020 15:30:39 -0400 Subject: [PATCH 3/7] Resolved issue of windows threading (design-docs pull request #17) --- reduce_sum_tutorial.Rmd | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd index a5c6d35..7b6d65a 100644 --- a/reduce_sum_tutorial.Rmd +++ b/reduce_sum_tutorial.Rmd @@ -38,11 +38,7 @@ Modifying a Stan model for use with `reduce_sum` has three steps. (3) Configure your compiler to enable multithreading. Then you can compile and sample as usual. -In the rest of this tutorial, we'll first build the model frist without and then with `reduce_sum`. - -[EDIT: Does Windows work with threading yet?] - -If you use Windows, it looks like multithreading isn't working yet (see ). So flip over to a Linux partition, if you have one. If you don't, hope is that an upcoming update to RTools will make things work. +In the rest of this tutorial, we'll first build the model first without and then with `reduce_sum`. ## Prepping the data @@ -252,17 +248,17 @@ make -j4 build (setting the argument to `j` equal to the number of cores you want to build in parallel with) -### Build and Run +To build the new model: ```bash -make redcards/logistic0 +make redcards/logistic1 ``` -To sample from the model: +To sample from the new model: ```bash cd redcard -time STAN_NUM_THREADS=8 ./logistic0 sample data file=redcard_input.R +time STAN_NUM_THREADS=8 ./logistic1 sample data file=redcard_input.R ``` The environment variable `STAN_NUM_THREADS` defines how many threads to use for this calculation. On this computer I have 8-cores, and so I choose to use 8 threads. From 2444a95e7bae655c1fb2618acc2cd465907844fd Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 25 Mar 2020 16:12:34 +0100 Subject: [PATCH 4/7] expand a bit here and there --- reduce_sum_tutorial.Rmd | 47 +++++++++++++++++++++++++++++++---------- 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd index 7b6d65a..cf4791a 100644 --- a/reduce_sum_tutorial.Rmd +++ b/reduce_sum_tutorial.Rmd @@ -19,7 +19,10 @@ features that make it easier to use parallelism. The main advantages to `reduce_sum` are: 1. `reduce_sum` has a more flexible argument interface, avoiding the packing and unpacking that was necessary with `map_rect`. -2. `reduce_sum` partitions the data for parallelization automatically. +2. `reduce_sum` partitions the data for parallelization automatically + such that the sharding pattern is chosen automatically. +3. `reduce_sum` allows to take advantage of fast vectorized + expressions easily. The drawbacks to `reduce_sum` are: @@ -32,7 +35,7 @@ This tutorial will highlight these advantages on a simple logistic regression us Modifying a Stan model for use with `reduce_sum` has three steps. -(1) Identify the portion of the model that you want to parallelize. `reduce_sum` is design specifically for speeding up log likelihood evaluations that are composed of a number of conditionally independent terms that can be computed in parallel then added together. +(1) Identify the portion of the model that you want to parallelize. `reduce_sum` is designed specifically for speeding up log likelihood evaluations that are composed of a number of conditionally independent terms that can be computed in parallel then added together. (2) Write a reduce function that can compute a given number of terms that must be added together. Pass this function (and other necessary arguments) to a call to `reduce_sum`. @@ -141,12 +144,14 @@ for(n in 1:N) { } ``` -Now it is clear that the calculation is the sum of a number of conditionally -independent Bernoulli log probability statements, which is the condition where -`reduce_sum` is useful. +Now it is clear that the calculation is the sum of a number of +conditionally independent Bernoulli log probability statements. So +whenever we need to calculate a large sum where each term is +independent of all others and associativity holds, then `reduce_sum` +is useful. -To use `reduce_sum`, a function must be written that can be used to compute -arbitrary subsets of the sums. +To use `reduce_sum`, a function must be written that can be used to +compute arbitrary partial sums out of the entire large sum. Using the reducer interface defined in [Reduce-Sum](#reduce-sum), such a function can be written like: @@ -165,7 +170,13 @@ functions { } ``` -And the likelihood statement in the model can now be written: +The first three arguments of the reducer function must always be +present. These encode the start and end index of the partial sum as +well as a slice of the entire data over which the reduction is to +be performed. Any additional arguments are optional. These optional +arguments are shared (read-only) among all computations. + +Thus, the likelihood statement in the model can now be written: ``` target += reduce_func(1, N, n_redcards, n_games, x, beta); // Sum terms 1 to N in the likelihood @@ -190,9 +201,10 @@ target += reduce_sum(reduce_func, n_redcards, grainsize, `reduce_sum` automatically breaks the sum into roughly `grainsize` sized pieces and computes them in parallel. `grainsize = 1` specifies that the grainsize should -be estimated automatically. +be estimated automatically, which is recommended to use as a starting +point for all analyses. -Making grainsize data (this makes it convenient to find a good one), the final model +Making `grainsize` data (this makes it convenient to find a good one), the final model should look like: ``` functions { @@ -271,7 +283,20 @@ user 2m15.234s sys 0m0.783s ``` -This gives roughly a 7x speedup. This model was particularly easy to parallelize because a lot of the arguments are data (and do not need to be autodiffed). If a model has a large number of arguments that are either defined in the parameters block, transformed parameters block, or model block, or do not do very much computation inside the reduce function, the speedup will be much more limited. +This gives roughly a 7x speedup. This model was particularly easy to +parallelize because a lot of the arguments are data (and do not need +to be autodiffed). If a model has a large number of arguments that are +either defined in the parameters block, transformed parameters block, +or model block, or do not do very much computation inside the reduce +function, the speedup will be much more limited. + +In case you turn this into a multi-level model where you introduce +groups of data for which you fit group specific parameters, then it is +often helpful to perform the slicing over the groups and use as +slicing argument of `reduce_sum` the group parameters themselves. This +leads to a more efficient parallelization, since the sliced argument +only gets doubled in memory while the shared arguments will be copied +as many times as there are splits of the large sum. ## More From c16b4e3d9196c6e6197e35092b0398e596692253 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 7 Apr 2020 14:40:04 -0400 Subject: [PATCH 5/7] Updated tutorial to reflect manual changes (design-doc #17) --- reduce_sum_tutorial.Rmd | 274 ++++++++++++++++++++-------------------- 1 file changed, 139 insertions(+), 135 deletions(-) diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd index 7b6d65a..50fcb70 100644 --- a/reduce_sum_tutorial.Rmd +++ b/reduce_sum_tutorial.Rmd @@ -5,86 +5,106 @@ output: html_document --- ```{r setup, include=FALSE} -library("rstan") +library("cmdstanr") ``` -This introduction to `reduce_sum` copies directly from Richard McElreath's [Multithreading and Map-Reduce in Stan 2.18.0: A Minimal Example](https://github.com/rmcelreath/cmdstan_map_rect_tutorial). +This introduction to `reduce_sum` copies directly from Richard McElreath's +[Multithreading and Map-Reduce in Stan 2.18.0: A Minimal Example](https://github.com/rmcelreath/cmdstan_map_rect_tutorial). ## Introduction -Stan 2.23 introduced `reduce_sum`, a new way to parallelize the execution of a single Stan chain across -multiple cores. This is in addition to the already existing `map_rect` utility, and introduces a number of -features that make it easier to use parallelism. +Stan 2.23 introduced `reduce_sum`, a new way to parallelize the execution of +a single Stan chain across multiple cores. This is in addition to the already +existing `map_rect` utility, and introduces a number of features that make it +easier to use parallelism: -The main advantages to `reduce_sum` are: +1. More flexible argument interface, avoiding the packing and + unpacking that is necessary with `map_rect`. +2. Partitions data for parallelization automatically (this is done manually + in `map_rect`). +3. Is easier to use. -1. `reduce_sum` has a more flexible argument interface, avoiding the packing and unpacking that was necessary with `map_rect`. -2. `reduce_sum` partitions the data for parallelization automatically. - -The drawbacks to `reduce_sum` are: - -1. `reduce_sum` requires that the result of the calculation be a scalar, while `map_rect` returns a list of vectors. -2. `reduce_sum` is only parallelized across multiple cores, not across multiple computers like `map_rect`. - -This tutorial will highlight these advantages on a simple logistic regression using `cmdstan`. +This tutorial will highlight these new features and demonstrate how to use +`reduce_sum` on a simple logistic regression using `cmdstanr`. ## Overview Modifying a Stan model for use with `reduce_sum` has three steps. -(1) Identify the portion of the model that you want to parallelize. `reduce_sum` is design specifically for speeding up log likelihood evaluations that are composed of a number of conditionally independent terms that can be computed in parallel then added together. +(1) Identify the portion of the model that you want to parallelize. `reduce_sum` +is design specifically for speeding up log likelihood evaluations that are +composed of a number of conditionally independent terms that can be computed +in parallel then added together. -(2) Write a reduce function that can compute a given number of terms that must be added together. Pass this function (and other necessary arguments) to a call to `reduce_sum`. +(2) Write a partial sum function that can compute a given number of terms in +that sum. Pass this function (and other necessary arguments) to a call +to `reduce_sum`. -(3) Configure your compiler to enable multithreading. Then you can compile and sample as usual. +(3) Configure your compiler to enable multithreading. Then you can compile and +sample as usual. -In the rest of this tutorial, we'll first build the model first without and then with `reduce_sum`. +To demonstrate this process, we'll first build a model first without, +`reduce_sum`, identify the part to be parallelized, and then work out how to use +`reduce_sum` to parallelize it. ## Prepping the data -Let's consider a simple logistic regression, just so the model doesn't get in the way. We'll use a reasonably big data table, the football red card data set from the recent crowdsourced data analysis project (). This data table contains 146,028 player-referee dyads. For each dyad, the table records the total number of red cards the referee assigned to the player over the observed number of games. +Let's consider a simple logistic regression. We'll use a reasonably large amount +of data, the football red card data set from the recent crowdsourced data +analysis project (). This data table contains +146,028 player-referee dyads. For each dyad, the table records the total number +of red cards the referee assigned to the player over the observed number of +games. -The ``RedcardData.csv`` file is provided in the repository here. Load the data in R, and let's take a look at the distribution of red cards: - ```R -d <- read.csv( "RedcardData.csv" , stringsAsFactors=FALSE ) -table( d$redCards ) -``` -```text -0 1 2 -144219 1784 25 -``` -The vast majority of dyads have zero red cards. Only 25 dyads show 2 red cards. These counts are our inference target. +`RedcardData.csv` is provided in the repository +[here](https://github.com/bbbales2/cmdstan_map_rect_tutorial/blob/reduce_sum/RedcardData.csv). +Load the data in R, and let's take a look at the distribution of red cards: -The motivating hypothesis behind these data is that referees are biased against darker skinned players. So we're going to try to predict these counts using the skin color ratings of each player. Not all players actually received skin color ratings in these data, so let's reduce down to dyads with ratings: - ```R -d2 <- d[ !is.na(d$rater1) , ] -out_data <- list( n_redcards=d2$redCards , n_games=d2$games , rating=d2$rater1 ) -out_data$N <- nrow(d2) +```{r} +d <- read.csv("RedcardData.csv", stringsAsFactors = FALSE) +table(d$redCards) ``` -This leaves us with 124,621 dyads to predict. -At this point, you are thinking: "But there are repeat observations on players and referees! You need some cluster variables in there in order to build a proper multilevel model!" You are right. But let's start simple. Keep your partial pooling on idle for the moment. +The vast majority of dyads have zero red cards. Only 25 dyads show 2 red cards. +These counts are our inference target. -Now to use these data with ``cmdstan``, we need to export the table to a file that ``cmdstan`` can read in. This is made easy: -```R -library(rstan) -stan_rdump(ls(out_data), "redcard_input.R", envir = list2env(out_data)) +The motivating hypothesis behind these data is that referees are biased against +darker skinned players. So we're going to try to predict these counts using the +skin color ratings of each player. Not all players actually received skin color +ratings in these data, so let's reduce down to dyads with ratings: + +```{r} +d2 <- d[!is.na(d$rater1),] +redcard_data <- list(n_redcards = d2$redCards, n_games = d2$games, rating = d2$rater1) +redcard_data$N <- nrow(d2) ``` +This leaves us with 124,621 dyads to predict. + +At this point, you are thinking: "But there are repeat observations on players +and referees! You need some cluster variables in there in order to build a +proper multilevel model!" You are right. But let's start simple. Keep your +partial pooling on idle for the moment. + ## Making the model -A Stan model for this problem is just a simple logistic (binomial) GLM. I'll assume you know Stan well enough already that I can just plop the code down here. It's contained in the ``logistic0.stan`` file in the repository. It's not a complex model: +A Stan model for this problem is just a simple logistic (binomial) GLM. I'll +assume you know Stan well enough already that I can just plop the code down +here. Save this code as `logistic0.stan` where you saved `RedcardData.csv` +(we'll assume this is your current working directory): + +```{stan} +data { + int N; + int n_redcards[N]; + int n_games[N]; + vector[N] rating; +} -``` - data { - int N; - int n_redcards[N]; - int n_games[N]; - vector[N] rating; - } parameters { vector[2] beta; } + model { beta ~ normal(0,1); n_redcards ~ binomial_logit(n_games , beta[1] + beta[2] * rating); @@ -93,49 +113,37 @@ model { ## Building and Running the Basic Model -Make a new folder called `redcard` inside your cmdstan folder folder and drop the `redcard_input.R` and `logistic0.stan` files in it. To build the Stan model, while still in the cmdstan folder, enter: - -```bash -make redcards/logistic0 -``` - -To sample from the model: - -```bash -cd redcard -time ./logistic0 sample data file=redcard_input.R -``` - -The `time` command at displays a a processor time report at the end. We'll want to compare this time to the time you get later, after you enable multithreading. On my computer, I get: +To build the model in `cmdstanr` run: -``` -real 2m11.467s -user 1m55.844s -sys 0m15.553s +```{r} +logistic0 <- cmdstan_model("logistic0.stan") ``` -The first line is the "real" time, the one you probably care about for benchmarking. +To sample from the model run: -By default, the output file is called `output.csv`. You'll want to go back into R to analyze it. Just open R in the same folder and enter: - -```R -library(rstan) -m0 <- read_stan_csv("output.csv") +```{r} +system.time(fit0 <- logistic0$sample(redcard_data, + num_cores = 1, + num_chains = 1)) ``` -Now you can proceed as usual to work with the samples. +The `elapsed` time is the time that we would have recorded if we were timing +this process with a stop-watch, so that is the one relevant to understanding +performance here (`system` time is time spent on system functions like I/O, +and `user` time is for parallel calculations). ## Rewriting the Model to Enable Multithreading -The key point to getting this calculation into `reduce_sum`, is recognizing that -the statement: +The key to porting this calculation to `reduce_sum`, is recognizing that the +statement: -``` +```{stan} n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating); ``` can be rewritten (up to a proportionality constant) as: -``` + +```{stan} for(n in 1:N) { target += binomial_logit_lpmf(n_redcards[n] | n_games[n], beta[1] + beta[2] * rating[n]) } @@ -146,53 +154,63 @@ independent Bernoulli log probability statements, which is the condition where `reduce_sum` is useful. To use `reduce_sum`, a function must be written that can be used to compute -arbitrary subsets of the sums. +arbitrary sections of this sum. -Using the reducer interface defined in [Reduce-Sum](#reduce-sum), such a function -can be written like: +Using the reducer interface defined in +[Reduce-Sum](https://mc-stan.org/docs/2_23/functions-reference/functions-reduce.html): -``` +```{stan} functions { - real reduce_func(int start, int end, + real partial_sum(int start, int end, int[] slice_n_redcards, int[] n_games, vector rating, vector beta) { return binomial_logit_lpmf(slice_n_redcards | n_games[start:end], - beta[1] + beta[2] * rating[start:end] ); + beta[1] + beta[2] * rating[start:end]); } } ``` -And the likelihood statement in the model can now be written: +The likelihood statement in the model can now be written: +```{stan} +target += partial_sum(1, N, n_redcards, n_games, rating, beta); // Sum terms 1 to N in the likelihood ``` -target += reduce_func(1, N, n_redcards, n_games, x, beta); // Sum terms 1 to N in the likelihood + +Equivalently it could be broken into two pieces and written like: +```{stan} +int M = N / 2; +target += partial_sum(1, M, n_redcards[1:M], n_games, rating, beta) // Sum terms 1 to M + partial_sum(M + 1, N, n_redcards[(M + 1):N], n_games, rating, beta); // Sum terms M + 1 to N ``` -In this example, `n_redcards` was chosen to be sliced over because there -is one term in the summation per value of `n_redcards`. Technically `n_games` -or `rating` would have worked just as well. Use whatever conceptually makes -the most sense. +By passing `partial_sum` to `reduce_sum`, we allow Stan to +automatically break up these calculations and do them in parallel. -Because `n_games` and `rating` are shared arguments, they must be subset -manually with `start:end`. +Notice the difference in how `n_redcards` is split in half (to reflect +which terms of the sum are being accumulated) and the rest of the arguments +(`n_games`, `x`, and `beta`) are left alone. This distinction is important +and more fully described in the User's Guide section on +[Reduce-sum](https://mc-stan.org/docs/2_23/stan-users-guide/reduce-sum.html). -With this function, `reduce_sum` can be used to automatically parallelize the +Given the partial sum function, `reduce_sum` can be used to automatically parallelize the likelihood: ``` -int grainsize = 100; +int grainsize = 1; target += reduce_sum(reduce_func, n_redcards, grainsize, - n_games, rating, beta); + n_games, rating, beta); ``` `reduce_sum` automatically breaks the sum into roughly `grainsize` sized pieces and computes them in parallel. `grainsize = 1` specifies that the grainsize should -be estimated automatically. +be estimated automatically (`grainsize` should be left at 1 unless specific tests +are done to +[pick a different one](https://mc-stan.org/docs/2_23/stan-users-guide/reduce-sum.html#reduce-sum-grainsize)). -Making grainsize data (this makes it convenient to find a good one), the final model +Making grainsize data (this makes it convenient to experiment with), the final model should look like: ``` functions { @@ -224,55 +242,41 @@ model { } ``` -Save this as `redcard/logistic1.stan`. +Save this as `logistic1.stan`. ## Running the Multithreaded Model -Now we're almost ready to go. First, make sure threading support is enabled in cmdstan (instructions [here](https://github.com/stan-dev/math/wiki/Threading-Support)). That basically means adding: +Compile this model with threading support: -```bash -STAN_THREADS=true +```{r} +logistic1 <- cmdstan_model("logistic1.stan", threads = TRUE) ``` -to `make/local`. If this was not set, before building the model, clean your install with a: +Set the number of threads each chain will use with `set_num_threads`: -```bash -make clean-all +```{r} +set_num_threads(8) # This computer has eight cores ``` -and then build it again using: +Run and time the model with: -```bash -make -j4 build -``` - -(setting the argument to `j` equal to the number of cores you want to build in parallel with) - -To build the new model: - -```bash -make redcards/logistic1 -``` - -To sample from the new model: - -```bash -cd redcard -time STAN_NUM_THREADS=8 ./logistic1 sample data file=redcard_input.R -``` - -The environment variable `STAN_NUM_THREADS` defines how many threads to use for this calculation. On this computer I have 8-cores, and so I choose to use 8 threads. - -The time command reports: - -``` -real 0m18.086s -user 2m15.234s -sys 0m0.783s +```{r} +system.time(fit0 <- logistic1$sample(redcard_data, + num_cores = 1, + num_chains = 1)) ``` -This gives roughly a 7x speedup. This model was particularly easy to parallelize because a lot of the arguments are data (and do not need to be autodiffed). If a model has a large number of arguments that are either defined in the parameters block, transformed parameters block, or model block, or do not do very much computation inside the reduce function, the speedup will be much more limited. +Again, `real` time is the time recorded as if by a stopwatch. This gives +roughly an X speedup. This model was particularly easy to parallelize because +a lot of the arguments are data (and do not need to be autodiffed). If a +model has a large number of arguments that are either defined in the +parameters block, transformed parameters block, or model block, or do not +do very much computation inside the reduce function, the speedup will be +much more limited. -## More +## More Information -There is a new section in the 2.23 User Manual and Function Reference that contains more details of using `reduce_sum`. +For a more detailed description of how `reduce_sum` works, see the +[User's Manual](https://mc-stan.org/docs/2_23/stan-users-guide/reduce-sum.html). +For a more complete description of the interface, see the +[Function Reference](https://mc-stan.org/docs/2_23/functions-reference/functions-reduce.html). From 59002ab6910167017d4084b38c309bdb42fe0eff Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 17 Apr 2020 13:44:40 -0400 Subject: [PATCH 6/7] Updated tutorial for latest reduce_sum changes. --- logistic0.stan | 8 +- logistic1.stan | 47 ++- reduce_sum_tutorial.Rmd | 112 ++++--- reduce_sum_tutorial.html | 676 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 767 insertions(+), 76 deletions(-) create mode 100644 reduce_sum_tutorial.html diff --git a/logistic0.stan b/logistic0.stan index c930f66..bd3dcf6 100644 --- a/logistic0.stan +++ b/logistic0.stan @@ -2,12 +2,14 @@ data { int N; int n_redcards[N]; int n_games[N]; - real rating[N]; + vector[N] rating; } parameters { vector[2] beta; } model { - beta ~ normal(0,1); - n_redcards ~ binomial_logit( n_games , beta[1] + beta[2] * to_vector(rating) ); + beta[1] ~ normal(0, 10); + beta[2] ~ normal(0, 1); + + n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating); } diff --git a/logistic1.stan b/logistic1.stan index dd9c64a..e07825a 100644 --- a/logistic1.stan +++ b/logistic1.stan @@ -1,40 +1,29 @@ functions { - vector lp_reduce( vector beta , vector theta , real[] xr , int[] xi ) { - int n = size(xr); - int y[n] = xi[1:n]; - int m[n] = xi[(n+1):(2*n)]; - real lp = binomial_logit_lpmf( y | m , beta[1] + to_vector(xr) * beta[2] ); - return [lp]'; + real partial_sum(int[] slice_n_redcards, + int start, int end, + int[] n_games, + vector rating, + vector beta) { + return binomial_logit_lpmf(slice_n_redcards | + n_games[start:end], + beta[1] + beta[2] * rating[start:end]); } -} +} data { int N; int n_redcards[N]; int n_games[N]; - real rating[N]; -} -transformed data { - // 7 shards - // M = N/7 = 124621/7 = 17803 - int n_shards = 7; - int M = N/n_shards; - int xi[n_shards, 2*M]; // 2M because two variables, and they get stacked in array - real xr[n_shards, M]; - // an empty set of per-shard parameters - vector[0] theta[n_shards]; - // split into shards - for ( i in 1:n_shards ) { - int j = 1 + (i-1)*M; - int k = i*M; - xi[i,1:M] = n_redcards[ j:k ]; - xi[i,(M+1):(2*M)] = n_games[ j:k ]; - xr[i] = rating[ j:k ]; - } + vector[N] rating; } parameters { vector[2] beta; } model { - beta ~ normal(0,1); - target += sum( map_rect( lp_reduce , beta , theta , xr , xi ) ); -} + int grainsize = 1; + + beta[1] ~ normal(0, 10); + beta[2] ~ normal(0, 1); + + target += reduce_sum(partial_sum, n_redcards, grainsize, + n_games, rating, beta); +} \ No newline at end of file diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd index 3161a44..66834c7 100644 --- a/reduce_sum_tutorial.Rmd +++ b/reduce_sum_tutorial.Rmd @@ -1,6 +1,6 @@ --- title: "Reduce Sum: A Minimal Example" -date: "23 Mar 2020" +date: "17 April 2020" output: html_document --- @@ -40,10 +40,10 @@ in parallel then added together. that sum. Pass this function (and other necessary arguments) to a call to `reduce_sum`. -(3) Configure your compiler to enable multithreading. Then you can compile and +(3) Configure your compiler to enable multithreading. Then compile and sample as usual. -To demonstrate this process, we'll first build a model first without, +To demonstrate this process, we'll first build a model first without `reduce_sum`, identify the part to be parallelized, and then work out how to use `reduce_sum` to parallelize it. @@ -93,21 +93,21 @@ assume you know Stan well enough already that I can just plop the code down here. Save this code as `logistic0.stan` where you saved `RedcardData.csv` (we'll assume this is your current working directory): -```{stan} +```{stan, output.var = "", eval = FALSE} data { int N; int n_redcards[N]; int n_games[N]; vector[N] rating; } - parameters { vector[2] beta; } - model { - beta ~ normal(0,1); - n_redcards ~ binomial_logit(n_games , beta[1] + beta[2] * rating); + beta[1] ~ normal(0, 10); + beta[2] ~ normal(0, 1); + + n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating); } ``` @@ -120,30 +120,35 @@ logistic0 <- cmdstan_model("logistic0.stan") ``` To sample from the model run: - + ```{r} -system.time(fit0 <- logistic0$sample(redcard_data, - num_cores = 1, - num_chains = 1)) +time0 = system.time(fit0 <- logistic0$sample(redcard_data, + num_cores = 1, + num_chains = 4, + refresh = 1000)) + +time0 ``` -The `elapsed` time is the time that we would have recorded if we were timing -this process with a stop-watch, so that is the one relevant to understanding +In this case, we just look at the performance of four chains running in sequence +on one core so we can compare that directly to four chains chain running with +threads. The `elapsed` time is the time that we would have recorded if we were +timing this process with a stop-watch, so that is the one relevant to understanding performance here (`system` time is time spent on system functions like I/O, -and `user` time is for parallel calculations). + and `user` time is for parallel calculations). ## Rewriting the Model to Enable Multithreading The key to porting this calculation to `reduce_sum`, is recognizing that the statement: - -```{stan} + +```{stan, output.var = "", eval = FALSE} n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating); ``` can be rewritten (up to a proportionality constant) as: - -```{stan} + +```{stan, output.var = "", eval = FALSE} for(n in 1:N) { target += binomial_logit_lpmf(n_redcards[n] | n_games[n], beta[1] + beta[2] * rating[n]) } @@ -160,11 +165,11 @@ arbitrary sections of this sum. Using the reducer interface defined in [Reduce-Sum](https://mc-stan.org/docs/2_23/functions-reference/functions-reduce.html): - -```{stan} + +```{stan, output.var = "", eval = FALSE} functions { - real partial_sum(int start, int end, - int[] slice_n_redcards, + real partial_sum(int[] slice_n_redcards, + int start, int end, int[] n_games, vector rating, vector beta) { @@ -176,23 +181,24 @@ functions { ``` The likelihood statement in the model can now be written: - -```{stan} -target += partial_sum(1, N, n_redcards, n_games, rating, beta); // Sum terms 1 to N in the likelihood + +```{stan, output.var = "", eval = FALSE} +target += partial_sum(n_redcards, 1, N, n_games, rating, beta); // Sum terms 1 to N in the likelihood ``` Equivalently it could be broken into two pieces and written like: -```{stan} + +```{stan, output.var = "", eval = FALSE} int M = N / 2; -target += partial_sum(1, M, n_redcards[1:M], n_games, rating, beta) // Sum terms 1 to M - partial_sum(M + 1, N, n_redcards[(M + 1):N], n_games, rating, beta); // Sum terms M + 1 to N +target += partial_sum(n_redcards[1:M], 1, M, n_games, rating, beta) // Sum terms 1 to M + partial_sum(n_redcards[(M + 1):N], M + 1, N, n_games, rating, beta); // Sum terms M + 1 to N ``` By passing `partial_sum` to `reduce_sum`, we allow Stan to automatically break up these calculations and do them in parallel. Notice the difference in how `n_redcards` is split in half (to reflect -which terms of the sum are being accumulated) and the rest of the arguments + which terms of the sum are being accumulated) and the rest of the arguments (`n_games`, `x`, and `beta`) are left alone. This distinction is important and more fully described in the User's Guide section on [Reduce-sum](https://mc-stan.org/docs/2_23/stan-users-guide/reduce-sum.html). @@ -200,9 +206,9 @@ and more fully described in the User's Guide section on Given the partial sum function, `reduce_sum` can be used to automatically parallelize the likelihood: -``` +```{stan, output.var = "", eval = FALSE} int grainsize = 1; -target += reduce_sum(reduce_func, n_redcards, grainsize, +target += reduce_sum(partial_sum, n_redcards, grainsize, n_games, rating, beta); ``` @@ -214,32 +220,34 @@ are done to Making grainsize data (this makes it convenient to experiment with), the final model is: -``` +```{stan, output.var = "", eval = FALSE} functions { - real reduce_func(int start, int end, - int[] slice_n_redcards, + real partial_sum(int[] slice_n_redcards, + int start, int end, int[] n_games, vector rating, vector beta) { return binomial_logit_lpmf(slice_n_redcards | n_games[start:end], - beta[1] + beta[2] * rating[start:end] ); + beta[1] + beta[2] * rating[start:end]); } } data { int N; int n_redcards[N]; int n_games[N]; - int grainsize; vector[N] rating; } parameters { vector[2] beta; } model { - beta ~ normal(0,1); + int grainsize = 1; - target += reduce_sum(reduce_func, n_redcards, grainsize, + beta[1] ~ normal(0, 10); + beta[2] ~ normal(0, 1); + + target += reduce_sum(partial_sum, n_redcards, grainsize, n_games, rating, beta); } ``` @@ -263,19 +271,35 @@ set_num_threads(8) # This computer has eight cores Run and time the model with: ```{r} -system.time(fit0 <- logistic1$sample(redcard_data, - num_cores = 1, - num_chains = 1)) +time1 = system.time(fit1 <- logistic1$sample(redcard_data, + num_cores = 1, + num_chains = 4, + refresh = 1000)) + +time1 +``` + +Again, `elapsed` time is the time recorded as if by a stopwatch. Computing the +ratios of the two times gives the speedup on eight cores of: + +```{r} +time0[["elapsed"]] / time1[["elapsed"]] ``` -Again, `real` time is the time recorded as if by a stopwatch. This gives -roughly an X speedup. This model was particularly easy to parallelize because +This model was particularly easy to parallelize because a lot of the arguments are data (and do not need to be autodiffed). If a model has a large number of arguments that are either defined in the parameters block, transformed parameters block, or model block, or do not do very much computation inside the reduce function, the speedup will be much more limited. +We can always get speedup in terms of effective sample size per time +by running multiple chains on different cores. `reduce_sum` is not a +replacement for that, and it is still important to run multiple chains +to check diagnostics. `reduce_sum` is a tool for speeding up single chain +calculations, which can be useful for model development and on computers with +large numbers of cores. + ## More Information For a more detailed description of how `reduce_sum` works, see the diff --git a/reduce_sum_tutorial.html b/reduce_sum_tutorial.html new file mode 100644 index 0000000..49fb7a4 --- /dev/null +++ b/reduce_sum_tutorial.html @@ -0,0 +1,676 @@ + + + + + + + + + + + + + + + +Reduce Sum: A Minimal Example + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + +

This introduction to reduce_sum copies directly from Richard McElreath’s Multithreading and Map-Reduce in Stan 2.18.0: A Minimal Example.

+
+

Introduction

+

Stan 2.23 introduced reduce_sum, a new way to parallelize the execution of a single Stan chain across multiple cores. This is in addition to the already existing map_rect utility, and introduces a number of features that make it easier to use parallelism:

+
    +
  1. More flexible argument interface, avoiding the packing and unpacking that is necessary with map_rect.
    2. +
  2. Partitions data for parallelization automatically (this is done manually in map_rect).
    3. +
  3. Is easier to use.
    4. +
+

This tutorial will highlight these new features and demonstrate how to use reduce_sum on a simple logistic regression using cmdstanr.

+
+
+

Overview

+

Modifying a Stan model for use with reduce_sum has three steps.

+
    +

  1. Identify the portion of the model that you want to parallelize. reduce_sum is design specifically for speeding up log likelihood evaluations that are composed of a number of conditionally independent terms that can be computed in parallel then added together.

    2. +

  2. Write a partial sum function that can compute a given number of terms in that sum. Pass this function (and other necessary arguments) to a call to reduce_sum.

    3. +

  3. Configure your compiler to enable multithreading. Then compile and sample as usual.

    4. +
+

To demonstrate this process, we’ll first build a model first without reduce_sum, identify the part to be parallelized, and then work out how to use reduce_sum to parallelize it.

+
+
+

Prepping the data

+

Let’s consider a simple logistic regression. We’ll use a reasonably large amount of data, the football red card data set from the recent crowdsourced data analysis project (https://psyarxiv.com/qkwst/). This data table contains 146,028 player-referee dyads. For each dyad, the table records the total number of red cards the referee assigned to the player over the observed number of games.

+

RedcardData.csv is provided in the repository here. Load the data in R, and let’s take a look at the distribution of red cards:

+
d <- read.csv("RedcardData.csv", stringsAsFactors = FALSE)
+table(d$redCards)
+
## 
+##      0      1      2 
+## 144219   1784     25
+

The vast majority of dyads have zero red cards. Only 25 dyads show 2 red cards. These counts are our inference target.

+

The motivating hypothesis behind these data is that referees are biased against darker skinned players. So we’re going to try to predict these counts using the skin color ratings of each player. Not all players actually received skin color ratings in these data, so let’s reduce down to dyads with ratings:

+
d2 <- d[!is.na(d$rater1),]
+redcard_data <- list(n_redcards = d2$redCards, n_games = d2$games, rating = d2$rater1)
+redcard_data$N <- nrow(d2)
+

This leaves us with 124,621 dyads to predict.

+

At this point, you are thinking: “But there are repeat observations on players and referees! You need some cluster variables in there in order to build a proper multilevel model!” You are right. But let’s start simple. Keep your partial pooling on idle for the moment.

+
+
+

Making the model

+

A Stan model for this problem is just a simple logistic (binomial) GLM. I’ll assume you know Stan well enough already that I can just plop the code down here. Save this code as logistic0.stan where you saved RedcardData.csv (we’ll assume this is your current working directory):

+
data {
+  int N;
+  int n_redcards[N];
+  int n_games[N];
+  vector[N] rating;
+}
+parameters {
+  vector[2] beta;
+}
+model {
+  beta[1] ~ normal(0, 10);
+  beta[2] ~ normal(0, 1);
+
+  n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating);
+}
+
+
+

Building and Running the Basic Model

+

To build the model in cmdstanr run:

+
logistic0 <- cmdstan_model("logistic0.stan")
+
## Compiling Stan program...
+
## A change in the compiler flags was found. Forcing recompilation.
+

To sample from the model run:

+
time0 = system.time(fit0 <- logistic0$sample(redcard_data,
+                                             num_cores = 1,
+                                             num_chains = 4,
+                                             refresh = 1000))
+
## Running MCMC with 4 chain(s) on 1 core(s)...
+## 
+## Running ./logistic0 'id=1' random 'seed=1988361638' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-1-8fd5e9.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 1 finished in 139.6 seconds.
+## Running ./logistic0 'id=2' random 'seed=2082795421' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-2-8fd5e9.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 2 finished in 149.3 seconds.
+## Running ./logistic0 'id=3' random 'seed=1823299851' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-3-8fd5e9.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 3 finished in 136.1 seconds.
+## Running ./logistic0 'id=4' random 'seed=1739568412' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-4-8fd5e9.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 4 finished in 145.6 seconds.
+## 
+## All 4 chains finished succesfully.
+## Mean chain execution time: 142.6 seconds.
+## Total execution time: 570.6 seconds.
+
time0
+
##    user  system elapsed 
+## 516.498  58.746 570.813
+

In this case, we just look at the performance of four chains running in sequence on one core so we can compare that directly to four chains chain running with threads. The elapsed time is the time that we would have recorded if we were timing this process with a stop-watch, so that is the one relevant to understanding performance here (system time is time spent on system functions like I/O, and user time is for parallel calculations).

+
+
+

Rewriting the Model to Enable Multithreading

+

The key to porting this calculation to reduce_sum, is recognizing that the statement:

+
n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating);
+

can be rewritten (up to a proportionality constant) as:

+
for(n in 1:N) {
+  target += binomial_logit_lpmf(n_redcards[n] | n_games[n], beta[1] + beta[2] * rating[n])
+}
+

Now it is clear that the calculation is the sum of a number of conditionally independent Bernoulli log probability statements. So whenever we need to calculate a large sum where each term is independent of all others and associativity holds, then reduce_sum is useful.

+

To use reduce_sum, a function must be written that can be used to compute arbitrary sections of this sum.

+

Using the reducer interface defined in Reduce-Sum:

+
functions {
+  real partial_sum(int[] slice_n_redcards,
+                   int start, int end,
+                   int[] n_games,
+                   vector rating,
+                   vector beta) {
+    return binomial_logit_lpmf(slice_n_redcards |
+                               n_games[start:end],
+                               beta[1] + beta[2] * rating[start:end]);
+  }
+}
+

The likelihood statement in the model can now be written:

+
target += partial_sum(n_redcards, 1, N, n_games, rating, beta); // Sum terms 1 to N in the likelihood
+

Equivalently it could be broken into two pieces and written like:

+
int M = N / 2;
+target += partial_sum(n_redcards[1:M], 1, M, n_games, rating, beta) // Sum terms 1 to M
+  partial_sum(n_redcards[(M + 1):N], M + 1, N, n_games, rating, beta); // Sum terms M + 1 to N
+

By passing partial_sum to reduce_sum, we allow Stan to automatically break up these calculations and do them in parallel.

+

Notice the difference in how n_redcards is split in half (to reflect which terms of the sum are being accumulated) and the rest of the arguments (n_games, x, and beta) are left alone. This distinction is important and more fully described in the User’s Guide section on Reduce-sum.

+

Given the partial sum function, reduce_sum can be used to automatically parallelize the likelihood:

+
int grainsize = 1;
+target += reduce_sum(partial_sum, n_redcards, grainsize,
+                     n_games, rating, beta);
+

reduce_sum automatically breaks the sum into roughly grainsize sized pieces and computes them in parallel. grainsize = 1 specifies that the grainsize should be estimated automatically (grainsize should be left at 1 unless specific tests are done to pick a different one).

+

Making grainsize data (this makes it convenient to experiment with), the final model is:

+
functions {
+  real partial_sum(int[] slice_n_redcards,
+                   int start, int end,
+                   int[] n_games,
+                   vector rating,
+                   vector beta) {
+    return binomial_logit_lpmf(slice_n_redcards |
+                               n_games[start:end],
+                               beta[1] + beta[2] * rating[start:end]);
+  }
+}
+data {
+  int N;
+  int n_redcards[N];
+  int n_games[N];
+  vector[N] rating;
+}
+parameters {
+  vector[2] beta;
+}
+model {
+  int grainsize = 1;
+
+  beta[1] ~ normal(0, 10);
+  beta[2] ~ normal(0, 1);
+
+  target += reduce_sum(partial_sum, n_redcards, grainsize,
+                       n_games, rating, beta);
+}
+

Save this as logistic1.stan.

+
+
+

Running the Multithreaded Model

+

Compile this model with threading support:

+
logistic1 <- cmdstan_model("logistic1.stan", threads = TRUE)
+
## Compiling Stan program...
+
## A change in the compiler flags was found. Forcing recompilation.
+

Set the number of threads each chain will use with set_num_threads:

+
set_num_threads(8) # This computer has eight cores
+

Run and time the model with:

+
time1 = system.time(fit1 <- logistic1$sample(redcard_data,
+                                             num_cores = 1,
+                                             num_chains = 4,
+                                             refresh = 1000))
+
## Running MCMC with 4 chain(s) on 1 core(s)...
+## 
+## Running ./logistic1 'id=1' random 'seed=1263537462' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-1-ff9dd7.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 1 finished in 21.4 seconds.
+## Running ./logistic1 'id=2' random 'seed=74174033' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-2-ff9dd7.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 2 finished in 20.3 seconds.
+## Running ./logistic1 'id=3' random 'seed=183191061' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-3-ff9dd7.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 3 finished in 21.5 seconds.
+## Running ./logistic1 'id=4' random 'seed=828878209' data \
+##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
+##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-4-ff9dd7.csv' 'refresh=1000' \
+##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
+##   'engaged=1'
+## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
+## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
+## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
+## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
+## Chain 4 finished in 20.7 seconds.
+## 
+## All 4 chains finished succesfully.
+## Mean chain execution time: 21.0 seconds.
+## Total execution time: 83.9 seconds.
+
time1
+
##    user  system elapsed 
+## 667.260   0.460  84.042
+

Again, elapsed time is the time recorded as if by a stopwatch. Computing the ratios of the two times gives the speedup on eight cores of:

+
time0[["elapsed"]] / time1[["elapsed"]]
+
## [1] 6.791997
+

This model was particularly easy to parallelize because a lot of the arguments are data (and do not need to be autodiffed). If a model has a large number of arguments that are either defined in the parameters block, transformed parameters block, or model block, or do not do very much computation inside the reduce function, the speedup will be much more limited.

+

We can always get speedup in terms of effective sample size per time by running multiple chains on different cores. reduce_sum is not a replacement for that, and it is still important to run multiple chains to check diagnostics. reduce_sum is a tool for speeding up single chain calculations, which can be useful for model development and on computers with large numbers of cores.

+
+
+

More Information

+

For a more detailed description of how reduce_sum works, see the User’s Manual. For a more complete description of the interface, see the Function Reference.

+
+ + + + +
+ + + + + + + + + + + + + + + From b24e005a654e5eca32eb34f77c565f87e4735582 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 17 Apr 2020 15:28:52 -0400 Subject: [PATCH 7/7] Added section on checking posterior samples --- reduce_sum_tutorial.Rmd | 14 ++++++ reduce_sum_tutorial.html | 93 +++++++++++++++++++++++----------------- 2 files changed, 68 insertions(+), 39 deletions(-) diff --git a/reduce_sum_tutorial.Rmd b/reduce_sum_tutorial.Rmd index 66834c7..8c700e0 100644 --- a/reduce_sum_tutorial.Rmd +++ b/reduce_sum_tutorial.Rmd @@ -300,6 +300,20 @@ to check diagnostics. `reduce_sum` is a tool for speeding up single chain calculations, which can be useful for model development and on computers with large numbers of cores. +We can do a quick check that these two methods are mixing with posterior. +When parallelizing a model is a good thing to do to make sure something is not +breaking: +```{r} +remotes::install_github("jgabry/posterior") +library(posterior) +summarise_draws(bind_draws(fit0$draws(), fit1$draws(), along = "chain")) +``` + +Ignoring `lp__` (remember the change to `bernoulli_logit_lpmf` means the +likelihoods between the two models are only proportional up to a constant), the Rhats +are less than 1.01, so our basic convergence checks pass. Unless something strange +is going on, we probably parallelized our model correctly! + ## More Information For a more detailed description of how `reduce_sum` works, see the diff --git a/reduce_sum_tutorial.html b/reduce_sum_tutorial.html index 49fb7a4..274eed4 100644 --- a/reduce_sum_tutorial.html +++ b/reduce_sum_tutorial.html @@ -438,53 +438,53 @@

Building and Running the Basic Model

refresh = 1000)) 
## Running MCMC with 4 chain(s) on 1 core(s)...
 ## 
-## Running ./logistic0 'id=1' random 'seed=1988361638' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-1-8fd5e9.csv' 'refresh=1000' \
+## Running ./logistic0 'id=1' random 'seed=209877525' data \
+##   'file=/tmp/RtmpRR2oFF/standata-741316b6f077.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic0-202004171418-1-753111.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 1 finished in 139.6 seconds.
-## Running ./logistic0 'id=2' random 'seed=2082795421' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-2-8fd5e9.csv' 'refresh=1000' \
+## Chain 1 finished in 135.1 seconds.
+## Running ./logistic0 'id=2' random 'seed=971698889' data \
+##   'file=/tmp/RtmpRR2oFF/standata-741316b6f077.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic0-202004171418-2-753111.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 2 finished in 149.3 seconds.
-## Running ./logistic0 'id=3' random 'seed=1823299851' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-3-8fd5e9.csv' 'refresh=1000' \
+## Chain 2 finished in 147.9 seconds.
+## Running ./logistic0 'id=3' random 'seed=1212530602' data \
+##   'file=/tmp/RtmpRR2oFF/standata-741316b6f077.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic0-202004171418-3-753111.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 3 finished in 136.1 seconds.
-## Running ./logistic0 'id=4' random 'seed=1739568412' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f75dbc488.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic0-202004171327-4-8fd5e9.csv' 'refresh=1000' \
+## Chain 3 finished in 156.0 seconds.
+## Running ./logistic0 'id=4' random 'seed=84501368' data \
+##   'file=/tmp/RtmpRR2oFF/standata-741316b6f077.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic0-202004171418-4-753111.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 4 finished in 145.6 seconds.
+## Chain 4 finished in 140.7 seconds.
 ## 
 ## All 4 chains finished succesfully.
-## Mean chain execution time: 142.6 seconds.
-## Total execution time: 570.6 seconds.
+## Mean chain execution time: 144.9 seconds. +## Total execution time: 579.7 seconds. 
time0
##    user  system elapsed 
-## 516.498  58.746 570.813
+## 525.802 58.636 579.945

In this case, we just look at the performance of four chains running in sequence on one core so we can compare that directly to four chains chain running with threads. The elapsed time is the time that we would have recorded if we were timing this process with a stop-watch, so that is the one relevant to understanding performance here (system time is time spent on system functions like I/O, and user time is for parallel calculations).

@@ -569,58 +569,73 @@

Running the Multithreaded Model

refresh = 1000)) 
## Running MCMC with 4 chain(s) on 1 core(s)...
 ## 
-## Running ./logistic1 'id=1' random 'seed=1263537462' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-1-ff9dd7.csv' 'refresh=1000' \
+## Running ./logistic1 'id=1' random 'seed=1000187990' data \
+##   'file=/tmp/RtmpRR2oFF/standata-7413274f09ae.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic1-202004171428-1-c0ef1a.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 1 finished in 21.4 seconds.
-## Running ./logistic1 'id=2' random 'seed=74174033' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-2-ff9dd7.csv' 'refresh=1000' \
+## Chain 1 finished in 22.0 seconds.
+## Running ./logistic1 'id=2' random 'seed=1751892397' data \
+##   'file=/tmp/RtmpRR2oFF/standata-7413274f09ae.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic1-202004171428-2-c0ef1a.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 2 finished in 20.3 seconds.
-## Running ./logistic1 'id=3' random 'seed=183191061' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-3-ff9dd7.csv' 'refresh=1000' \
+## Chain 2 finished in 24.1 seconds.
+## Running ./logistic1 'id=3' random 'seed=1061813517' data \
+##   'file=/tmp/RtmpRR2oFF/standata-7413274f09ae.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic1-202004171428-3-c0ef1a.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 3 finished in 21.5 seconds.
-## Running ./logistic1 'id=4' random 'seed=828878209' data \
-##   'file=/tmp/RtmpVB5mWC/standata-6f3f7f443b18.json' output \
-##   'file=/tmp/RtmpVB5mWC/logistic1-202004171337-4-ff9dd7.csv' 'refresh=1000' \
+## Chain 3 finished in 22.3 seconds.
+## Running ./logistic1 'id=4' random 'seed=756761795' data \
+##   'file=/tmp/RtmpRR2oFF/standata-7413274f09ae.json' output \
+##   'file=/tmp/RtmpRR2oFF/logistic1-202004171428-4-c0ef1a.csv' 'refresh=1000' \
 ##   'method=sample' 'save_warmup=0' 'algorithm=hmc' 'engine=nuts' adapt \
 ##   'engaged=1'
 ## Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
 ## Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
 ## Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
 ## Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
-## Chain 4 finished in 20.7 seconds.
+## Chain 4 finished in 22.1 seconds.
 ## 
 ## All 4 chains finished succesfully.
-## Mean chain execution time: 21.0 seconds.
-## Total execution time: 83.9 seconds.
+## Mean chain execution time: 22.6 seconds. +## Total execution time: 90.5 seconds. 
time1
##    user  system elapsed 
-## 667.260   0.460  84.042
+## 721.053 0.675 90.662

Again, elapsed time is the time recorded as if by a stopwatch. Computing the ratios of the two times gives the speedup on eight cores of:

time0[["elapsed"]] / time1[["elapsed"]]
-
## [1] 6.791997
+
## [1] 6.396781

This model was particularly easy to parallelize because a lot of the arguments are data (and do not need to be autodiffed). If a model has a large number of arguments that are either defined in the parameters block, transformed parameters block, or model block, or do not do very much computation inside the reduce function, the speedup will be much more limited.

We can always get speedup in terms of effective sample size per time by running multiple chains on different cores. reduce_sum is not a replacement for that, and it is still important to run multiple chains to check diagnostics. reduce_sum is a tool for speeding up single chain calculations, which can be useful for model development and on computers with large numbers of cores.

+

We can do a quick check that these two methods are mixing with posterior. When parallelizing a model is a good thing to do to make sure something is not breaking:

+
remotes::install_github("jgabry/posterior")
+
## Skipping install of 'posterior' from a github remote, the SHA1 (2f4ba9e3) has not changed since last install.
+##   Use `force = TRUE` to force installation
+
library(posterior)
+
## This is posterior version 0.0.2
+
summarise_draws(bind_draws(fit0$draws(), fit1$draws(), along = "chain"))
+
## # A tibble: 3 x 10
+##   variable     mean   median      sd     mad       q5      q95  rhat ess_bulk
+##   <chr>       <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
+## 1 lp__     -9.05e+3 -9.06e+3 1.20e+3 1.78e+3 -1.03e+4 -7.85e+3  1.71     12.5
+## 2 beta[1]  -5.53e+0 -5.53e+0 3.37e-2 3.37e-2 -5.59e+0 -5.48e+0  1.00   2566. 
+## 3 beta[2]   2.89e-1  2.88e-1 8.10e-2 8.05e-2  1.55e-1  4.21e-1  1.00   2610. 
+## # … with 1 more variable: ess_tail <dbl>
+

Ignoring lp__ (remember the change to bernoulli_logit_lpmf means the likelihoods between the two models are only proportional up to a constant), the Rhats are less than 1.01, so our basic convergence checks pass. Unless something strange is going on, we probably parallelized our model correctly!

More Information