From 24467fd7f022182afa9e805bdd154363541d3aee Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 23 Mar 2020 11:16:16 -0400 Subject: [PATCH 01/13] Added initial attempt at docs for reduce-sum (design-doc pull request #17) --- .../higher-order_functions.Rmd | 61 ++++- src/reference-manual/expressions.Rmd | 26 ++- src/stan-users-guide/_bookdown.yml | 2 +- ...{map-reduce.Rmd => parallel-computing.Rmd} | 221 ++++++++++++++++-- 4 files changed, 281 insertions(+), 29 deletions(-) rename src/stan-users-guide/{map-reduce.Rmd => parallel-computing.Rmd} (58%) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index 83cfa0049..38fea56d1 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -382,8 +382,67 @@ Internally the 1D integrator uses the double-exponential methods in the Boost 1D The gradients of the integral are computed in accordance with the Leibniz integral rule. Gradients of the integrand are computed internally with Stan's automatic differentiation. +## Parallel Reduce-Sum Function {#functions-reduce} -## Higher-Order Map {#functions-map} +Stan provides a higher-order ```reduce_sum``` function for parallelizing operations that can be represented as a reduce (by summation) over a sequence of terms. + +### Reduce-Sum Function + +The reduce sum function operates on a reducing function, a list of +sliced arguments (one for each term in the parallel-for), a recommended grainsize, +and a set of shared arguments. + + +\index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] sliced\_arg, int grainsize, T1 arg1, T2 arg2, ...): real}|hyperpage} + +`real` **`reduce_sum`**`(F f, T[] sliced_arg, int grainsize, T1 arg1, T2 arg2, ...)`
\newline + +Return the equivalent of `f(1, size(sliced_arg), arg1, arg2, ...)`, but compute +the result in parallel by breaking `sliced_arg` into pieces and computing each piece +in a different thread. `arg1, arg2, ...` are shared between all terms in the sum. + +* *`f`*: function literal referring to a function specifying the reduce operation with signature `(int, int, T[], T1, T2, ...):real` +The arguments represent + + (1) the index of the first term of the reduction, + + (2) the index of the last term of the reduction, + + (3) the subset `sliced_arg` this reduce is responsible for computing, + + (4) first shared argument, + + (5) second shared argument, + + ... the rest of the shared arguments. + +* *`sliced_args`*: array of non-shared arguments, one for each term of the reduction, array of `T`, where `T` can be any type, +* *`grainsize`*: recommented number of terms in each reduce call, set to zero to estimate automatically, type `int`, +* *`arg1`*: first (optional) shared argument, type `T1`, where `T1` can be any type +* *`arg2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, +* *`...`*: remainder of shared arguments, each of which can be any type. + +### Specifying the Reduce Function + +The reduce function must have the following signature where the types T, and the +types of all the variadic arguments (`T1`, `T2`, ...) match those of the original +`reduce_sum` call. + +``` +(int start, int end, T[] subset_sliced_arg, T1 arg1, T2 arg2, ...):real +``` + +The reduce function returns the sum of the `start` to `end` terms of the overall +calculations. The arguments to the reduce function are: + +* *`start`*, the index of the first term of the reduction, type `int` + +* *`end`*, the index of the last term of the reduction (inclusive), type `int` + +* *`subset_sliced_arg`*, the subset `sliced_arg` this reduce is responsible for computing, type `T[]`, where `T` matches the type of `sliced_arg` in `reduce_sum` + +* *`arg1`*, first shared argument, type `T1`, matching type of `arg1` in `reduce_sum` + +* *`arg2`*, second shared argument, type `T2`, matching type of `arg2` in `reduce_sum` + +* *`...`*, remainder of shared arguments, with types matching those in `reduce_sum` + + +## Parallel Map-Rect Function {#functions-map} Stan provides a higher-order map function. This allows map-reduce functionality to be coded in Stan as described in the user's guide. diff --git a/src/reference-manual/expressions.Rmd b/src/reference-manual/expressions.Rmd index bc31d081e..0b7c6f57d 100644 --- a/src/reference-manual/expressions.Rmd +++ b/src/reference-manual/expressions.Rmd @@ -1118,27 +1118,33 @@ literals.* `integrate_1d`, | `real, real, real[]` | `real[], int[]` | `real` `integrate_ode_X`, | `real, real[], real[]` | `real[], int[]` | `real[]` `map_rect` | `vector, vector` | `real[], int[]` | `vector` +`reduce_sum` | ```T[], T1, T2, ...``` | | `real` -For example, the rectangular mapping function might be used in the -following way to compute the log likelihood of a hierarchical model. +`T`, `T1`, `T2`, and the types of `...` can be any Stan type. + +For example, the `integrate_ode_rk45` function can be used to integrate +differential equations in Stan: ```stan functions { - vector foo_ll(vector phi, vector theta, real[] x_r, int[] x_i) { + real[] foo(real t, real[] y, real[] theta, real[] x_r, int[] x_i) { ... ... -vector[11] phi; -vector[2] thetas[N]; -real x_rs[N, 5]; -real x_is[N, 0]; +int T; +real y0[2]; +real t0; +real ts[T]; +real theta[1]; +real x_r[0]; +int x_i[0]; ... -target += sum(map_rect(foo_ll, phi, thetas, x_rs, x_is)); +real y_hat[T, 2] = integrate_ode_rk45(foo, y0, t0, ts, theta, x_r, x_i); ``` The function argument is `foo`, the name of the user-defined function; as shown in the [higher-order functions table](#higher-order-functions), `foo` -takes two vectors, a real array, and an integer array as arguments and -returns a vector. +takes a real, three more real arrays, and an integer array as arguments and +returns a real array. ### Functions Passed by Reference {-} diff --git a/src/stan-users-guide/_bookdown.yml b/src/stan-users-guide/_bookdown.yml index 5f30b91a0..55c0e1b6b 100644 --- a/src/stan-users-guide/_bookdown.yml +++ b/src/stan-users-guide/_bookdown.yml @@ -29,7 +29,7 @@ rmd_files: [ "problematic-posteriors.Rmd", "reparameterization.Rmd", "efficiency-tuning.Rmd", - "map-reduce.Rmd", + "parallel-computing.Rmd", "part-appendices.Rmd", "style-guide.Rmd", diff --git a/src/stan-users-guide/map-reduce.Rmd b/src/stan-users-guide/parallel-computing.Rmd similarity index 58% rename from src/stan-users-guide/map-reduce.Rmd rename to src/stan-users-guide/parallel-computing.Rmd index c557c123a..c505eebf9 100644 --- a/src/stan-users-guide/map-reduce.Rmd +++ b/src/stan-users-guide/parallel-computing.Rmd @@ -1,12 +1,201 @@ -# Map-Reduce {#map-reduce.chapter} +# Parallel Computing {#parallel-computing.chapter} + +Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and and `map_rect`. + +The main differences are: +1. `reduce_sum` requires the result of the calculation to be a scalar, while `map_rect` returns a list of vectors +2. `reduce_sum` has a more flexible interface and can accept arbitrary Stan types as arguments, `map_rect` is more restrictive on what arguments can be and how they are shaped +3. `map_rect` can parallelize work over multiple computers or a single computer, while `reduce_sum` works only on a single computer +4. `map_rect` requires work to be broken into pieces manually, while `reduce_sum` mostly automates this + +## Reduce-Sum { #reduce-sum } + +```reduce_sum``` is a tool for parallelizing operations that can be represented as a parallel-for combined with a sum (that returns a scalar). + +In terms of probabilistic models, an example of this comes up when N terms in a likelihood combined multiplicatively and can be computed independently of each other (where independence here is in the computational sense, not necessarily the statistical sense). In this case, computing the log density means computing the sum of a number of terms that can each be computed separately. + +```reduce_sum``` is not useful when there are dependencies between the terms. This can happen, for instance, if there were N terms in a Gaussian process likelihood. ```reduce_sum``` will not be useful for accelerating this. + +If for a set of input arguments, ```args0, args1, args2, ...``` and a scalar function ```f```, the log likelihood can be computed as: + +```f(args0) + f(args1) + f(args2) + ...``` + +then this calculation can be written as a reduction over the set of arguments. If this reducing function is called ```reduce```, then it would need to perform the operations: + +```reduce({ args0, args1, args2, ... }) = f(args0) + f(args1) + f(args2) + ...``` + +If the user can write a function like ```reduce```, then it is trivial for us to provide a function to automatically parallelize the calculations. + +Again, if the set of work is represented as a list of arguments ``{ args0, args1, args2, ... }```, then mathematically it is possible to rewrite this sum with any combination of partial-reduces. + +For instance, the sum can be written: + +1. ```reduce({ args0, args1, args2, ... })```, summing over all arguments, using one reduce function +2. ```reduce({ args0, ..., args(M - 1) }) + reduce({ argsM, args2, ...})```, computing the first M terms separately from the rest +3. ```reduce({ args0 }) + reduce({ args1 }) + reduce({ args2 }) + ...```, computing each term individually and summing them + +The first function call is completely serial, the second can be parallelized over two workers, and the last can be parallelized over as many workers as there are arguments. Depending on how the list is sliced up, it is possible to parallelize this calculation over many workers. + +```reduce_sum``` is the tool that will allow us to automatically parallelize these reduce operations (and sum them together). + +To implement this efficiently in Stan, the individual arguments are split into two types. The first are the arguments that are specific to each term in the reduction. These are called the sliced arguments (because we will slice these up to determine how to distribute the work). The second type of arguments are shared arguments, and are information that is shared in the computation of every term in the sum. + +Given this, the signature for reduce_sum is: + +``` +real reduce_sum(F func, T[] sliced_arg, int grainsize, T1 arg1, T2 arg2, ...) +``` + +1. ```func``` - The user-defined reduce function +2. ```sliced_arg``` - An array of any type, with each element of the array corresponding to a term of the final summation (the length of ```sliced_arg``` is the total number of terms to sum) +3. ```grainsize``` - A hint to the runtime of how many terms of the summation to compute in each reduction +4-. ```arg1, arg2, ...``` - All the arguments that are to be shared in the calculation of every term in the sum + +The user-defined reduce function is slightly different: + +``` +real func(int start, int end, T[] subset_sliced_arg, T1 arg1, T2 arg2, ...) +``` + +and takes the arguments: +1. ```start``` - An integer specifying the first element of the sequence of terms this reduce call is responsible for computing +2. ```end``` - An integer specifying the last element of the sequence of terms this reduce call is responsible for computing +3. ```subset_sliced_arg``` - The subset of sliced_arg for which this reduce is responsible (```sliced_arg[start:end]```) +4-. ```arg1, arg2, ...``` all the shared arguments -- passed on without modification from the reduce_sum call + +The user-provided function ```func``` is expect to compute the ```start``` through ```end``` terms of the overall sum, accumulate them, and return that value. The user function is only passed the subset ```sliced_arg[start:end]``` of sliced arg (as ```subset_sliced_arg```). ```start``` and ```end``` are passed so that ```func``` can index any of the ```argM``` appropriately. The trailing arguments ```argM``` are passed without modification to every call of ```func```. + +An overall call to ```reduce_sum``` can be replaced by either: + +``` +real sum = func(1, size(sliced_arg), sliced_arg, arg1, arg2, ...) +``` + +or (modulo differences due to rearrangements of summations) the code: + +``` +real sum = 0.0; +for(i in 1:size(sliced_arg)) { + sum = sum + func(i, i, { sliced_arg[i] }, arg1, arg2, ...); +} +``` + +### Example: Logistic Regression + +Logistic Regression is a useful example to clarify both the syntax +and semantics of `reduce_sum` and how it can be used to speed up a typical +model. + +A basic logistic regression can be coded in Stan as: + +``` +data { + int N; + int y[N]; + vector[N] x; +} +parameters { + vector[2] beta; +} +model { + beta ~ std_normal(); + y ~ bernoulli_logit(beta[1] + beta[2] * x); +} +``` + +In this model predictions are made about the `N` outputs `y` using the +covariate `x`. The intercept and slope of the linear equation are to be estimated. + +The key point to getting this calculation into `reduce_sum`, is recognizing that +the statement: + +``` +y ~ bernoulli_logit(beta[1] + beta[2] * x); +``` + +can be rewritten (up to a proportionality constant) as: +``` +for(n in 1:N) { + target += bernoulli_logit_lpmf(y[n] | beta[1] + beta[2] * x[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 reducer_func(int start, int end, + int[] subset_y, + vector x, + vector beta) { + return bernoulli_logit_lpmf(subset_y | beta[1] + beta[2] * x[start:end]); + } +} +``` + +And the likelihood statement in the model can now be written: + +``` +target += reducer_fun(1, N, y, x, beta); // Sum terms 1 to N in the likelihood +``` + +In this example, `y` was chosen to be sliced over because there +is one term in the summation per value of `y`. Technically `x` would have +worked as well. Use whatever conceptually makes the most sense. + +Because `x` is a shared argument, it is subset manually with `start:end`. + +With this function, `reduce_sum` can be used to automatically parallelize the +likelihood: + +``` +int grainsize = 100; +target += reduce_sum(reducer_func, y, + grainsize, + x, beta); +``` + +`reduce_sum` automatically breaks the sum into roughly `grainsize` sized pieces +and computes them in parallel. `grainsize = 0` specifies that the grainsize should +be estimated automatically. + +### Picking the Grainsize + +The `grainsize` is a recommendation on how large each piece of parallel work is +(how many terms it contains). If zero, it will be chosen automatically, but it +is probably best to choose this manually for each model. + +To figure out an appropriate grainsize, think about how many terms are in the summation +and on how many cores the model should run. If there are `N` terms and `M` cores, +run a quick test model with `grainsize` set roughly to `N / M`. Record the time, cut the +grainsize in half, and run the test again. Repeat this iteratively until the model runtime +begins to increase. This is a suitable grainsize for the model, because this ensures the +caculations can be carried out with the most parallelism without losing too much efficiency. + +For instance, in a model with `N=10000` and `M = 4`, start with `grainsize = 25000`, and +sequentially try `grainsize = 12500`, `grainsize = 6250`, etc. + +It is important to repeat this process until performance gets worse! It is possible +after many halvings nothing happens, but there might still be a small grainsize that performs better. +Even if a sum has many tens of thousands of terms, depending on the internal calculations, a `grainsize` +of thirty or forty or smaller might be the best, and it is difficult to predict this behavior. +Without doing these halvings until performance actually gets worse, it is easy to miss this. + +## Map-Rect Map-reduce allows large calculations (e.g., log likelihoods) to be broken into components which may be calculated modularly (e.g., data blocks) and combined (e.g., by summation and incrementing the target log density). -## Overview of Map-Reduce - A _map function_ is a higher-order function that applies an argument function to every member of some collection, returning a collection of the results. For example, mapping the square function, @@ -21,14 +210,14 @@ are summation (with the return being a single value) or sorting (with the return being a sorted sequence). The combination of mapping and reducing is so common it has its own name, _map-reduce_. -## Map Function +### Map Function In order to generalize the form of functions and results that are possible and accommodate both parameters (which need derivatives) and data values (which don't), Stan's map function operates on more than just a sequence of inputs. -### Map Function Signature {-} +#### Map Function Signature {-} Stan's map function has the following signature @@ -68,7 +257,7 @@ Although `f` will often return a vector of size one, the built-in flexibility allows general multivariate functions to be mapped, even raggedly. -### Map Function Semantics {-} +#### Map Function Semantics {-} Stan's map function applies the function `f` to the shared parameters along with one element each of the job parameters, real @@ -88,13 +277,13 @@ of each application of `f` is a vector, and the sequence of `N` vectors is concatenated together to return a single vector. -## Example: Mapping Logistic Regression +### Example: Logistic Regression An example should help to clarify both the syntax and semantics of the mapping operation and how it may be combined with reductions built into Stan to provide a map-reduce implementation. -### Unmapped Logistic Regression {-} +#### Unmapped Logistic Regression {-} Consider the following simple logistic regression model, which is coded unconventionally to accomodate direct translation to a mapped @@ -124,7 +313,7 @@ because the argument is on the logit scale---it implicitly applies the inverse logit function to map the argument to a probability. -### Mapped Logistic Regression {-} +#### Mapped Logistic Regression {-} The unmapped logistic regression model described in the previous subsection may be implemented using Stan's rectangular mapping @@ -196,7 +385,7 @@ shared across shards. There are no shard-specific parameters, so the array of job-specific parameters `theta` contains only empty vectors. -## Example: Hierarchical Logistic Regression +### Example: Hierarchical Logistic Regression Consider a hierarchical model of American presidential voting behavior based on state of residence.^[This example is a simplified form of the model @@ -220,7 +409,7 @@ The slopes and intercepts get hierarchical priors, \beta_k &\sim \textsf{normal}(\mu_{\beta}, \sigma_{\beta}) \end{align*} -### Unmapped Implementation {-} +#### Unmapped Implementation {-} This model can be coded up in Stan directly as follows. @@ -263,9 +452,7 @@ for (n in 1:N) y[n] ~ bernoulli_logit(beta[kk[n], 1] + beta[kk[n], 2] * x[n]); ``` - - -### Mapped Implementation {-} +#### Mapped Implementation {-} The mapped version of the model will map over the states `K`. This means the group-level parameters, real data, and integer-data @@ -368,13 +555,13 @@ parallelization with the version that leaves the prior in the model block. -## Ragged Inputs and Outputs +### Ragged Inputs and Outputs The previous examples included rectangular data structures and single outputs. Despite the name, this is not technically required by `map_rect`. -### Ragged Inputs {-} +#### Ragged Inputs {-} If each group has a different number of observations, then the rectangular data structures for predictors and outcomes will need to @@ -383,7 +570,7 @@ structure will need to be passed as integer data. This holds for shards with varying numbers of parameters as well as varying numbers of data points. -### Ragged Outputs {-} +#### Ragged Outputs {-} The output of each mapped function is concatenated in order of inputs to produce the output of `map_rect`. When every shard returns a singleton From 9c3af806dff8d9f0e4c7b304d18f757b1c34e6c2 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 23 Mar 2020 12:34:56 -0400 Subject: [PATCH 02/13] Added space to make ordered list format correctly in reduce-sum section (design-docs pull req #17) --- src/stan-users-guide/parallel-computing.Rmd | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stan-users-guide/parallel-computing.Rmd b/src/stan-users-guide/parallel-computing.Rmd index c505eebf9..705294a26 100644 --- a/src/stan-users-guide/parallel-computing.Rmd +++ b/src/stan-users-guide/parallel-computing.Rmd @@ -3,6 +3,7 @@ Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and and `map_rect`. The main differences are: + 1. `reduce_sum` requires the result of the calculation to be a scalar, while `map_rect` returns a list of vectors 2. `reduce_sum` has a more flexible interface and can accept arbitrary Stan types as arguments, `map_rect` is more restrictive on what arguments can be and how they are shaped 3. `map_rect` can parallelize work over multiple computers or a single computer, while `reduce_sum` works only on a single computer @@ -164,7 +165,7 @@ target += reduce_sum(reducer_func, y, ``` `reduce_sum` automatically breaks the sum into roughly `grainsize` sized pieces -and computes them in parallel. `grainsize = 0` specifies that the grainsize should +and computes them in parallel. `grainsize = 1` specifies that the grainsize should be estimated automatically. ### Picking the Grainsize From 3fcff92d370bcaf73d3f77bd0a0798869d069d41 Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 26 Mar 2020 15:46:08 -0400 Subject: [PATCH 03/13] Updated reduce_sum docs to reflect design-doc changes (design-doc pull request #17) --- .../higher-order_functions.Rmd | 71 +++++---- src/stan-users-guide/parallel-computing.Rmd | 148 ++++++++++++------ 2 files changed, 139 insertions(+), 80 deletions(-) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index 38fea56d1..b2c9d3845 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -10,6 +10,7 @@ if (knitr::is_html_output()) { cat(' * Algebraic Equation Solver\n') cat(' * Ordinary Differential Equation (ODE) Solvers\n') cat(' * 1D Integrator\n') +cat(' * Reduce-Sum\n') cat(' * Higher-Order Map\n') } ``` @@ -382,67 +383,75 @@ Internally the 1D integrator uses the double-exponential methods in the Boost 1D The gradients of the integral are computed in accordance with the Leibniz integral rule. Gradients of the integrand are computed internally with Stan's automatic differentiation. -## Parallel Reduce-Sum Function {#functions-reduce} +## Reduce-Sum Function {#functions-reduce} -Stan provides a higher-order ```reduce_sum``` function for parallelizing operations that can be represented as a reduce (by summation) over a sequence of terms. +Stan provides a higher-order `reduce_sum` function for parallelizing operations that can be represented as a sum of a function, `g: U -> real`, evaluated at each element of a list of type `U[]`, `{ x1, x2, ... }`. That is: -### Reduce-Sum Function +```g(x1) + g(x2) + ...``` -The reduce sum function operates on a reducing function, a list of -sliced arguments (one for each term in the parallel-for), a recommended grainsize, -and a set of shared arguments. +`reduce_sum` doesn't work on `g` itself, but takes a partial sum function, `f: U[] -> real`, where: - -\index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] sliced\_arg, int grainsize, T1 arg1, T2 arg2, ...): real}|hyperpage} +```f({ x1 }) = g(x1)``` +```f({ x1, x2 }) = g(x1) + g(x2)``` +```f({ x1, x2, ... }) = g(x1) + g(x2) + ...``` -`real` **`reduce_sum`**`(F f, T[] sliced_arg, int grainsize, T1 arg1, T2 arg2, ...)`
\newline +### The Reduce-sum Function -Return the equivalent of `f(1, size(sliced_arg), arg1, arg2, ...)`, but compute -the result in parallel by breaking `sliced_arg` into pieces and computing each piece -in a different thread. `arg1, arg2, ...` are shared between all terms in the sum. +The `reduce_sum` function takes a partial sum function, an array argument x +(with one for each term in the sum), a recommended grainsize, and a set of shared arguments and +parallelizes the resultant sum. -* *`f`*: function literal referring to a function specifying the reduce operation with signature `(int, int, T[], T1, T2, ...):real` + +\index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] x, int grainsize, T1 s1, T2 s2, ...): real}|hyperpage} + +`real` **`reduce_sum`**`(F f, T[] x, int grainsize, T1 s1, T2 s2, ...)`
\newline + +Return the equivalent of `f(1, size(x), x, s1, s2, ...)`, but compute +the result in parallel by breaking the array `x` into pieces and computing each piece +in a different thread. `s1, s2, ...` are shared between all terms in the sum. + +* *`f`*: function literal referring to a function specifying the partial sum operation with signature `(int, int, T[], T1, T2, ...):real` The arguments represent - + (1) the index of the first term of the reduction, - + (2) the index of the last term of the reduction, - + (3) the subset `sliced_arg` this reduce is responsible for computing, + + (1) the index of the first term of the partial sum, + + (2) the index of the last term of the partial sum, + + (3) the subset `x` this reduce is responsible for computing, + (4) first shared argument, + (5) second shared argument, + ... the rest of the shared arguments. -* *`sliced_args`*: array of non-shared arguments, one for each term of the reduction, array of `T`, where `T` can be any type, -* *`grainsize`*: recommented number of terms in each reduce call, set to zero to estimate automatically, type `int`, -* *`arg1`*: first (optional) shared argument, type `T1`, where `T1` can be any type -* *`arg2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, +* *`x`*: array of `T`, one for each term of the reduction, `T` can be any type, +* *`grainsize`*: recommented number of terms in each reduce call, set to one to estimate automatically, type `int`, +* *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type +* *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, * *`...`*: remainder of shared arguments, each of which can be any type. -### Specifying the Reduce Function +### The Partial-sum Function -The reduce function must have the following signature where the types T, and the -types of all the variadic arguments (`T1`, `T2`, ...) match those of the original +The partial sum function must have the following signature where the types `T`, and the +types of all the shared arguments (`T1`, `T2`, ...) match those of the original `reduce_sum` call. ``` -(int start, int end, T[] subset_sliced_arg, T1 arg1, T2 arg2, ...):real +(int start, int end, T[] x_subset, T1 s1, T2 s2, ...):real ``` -The reduce function returns the sum of the `start` to `end` terms of the overall +The reduce function returns the sum of the `start` to `end` terms (inclusive) of the overall calculations. The arguments to the reduce function are: -* *`start`*, the index of the first term of the reduction, type `int` +* *`start`*, the index of the first term of the partial sum, type `int` -* *`end`*, the index of the last term of the reduction (inclusive), type `int` +* *`end`*, the index of the last term of the partial sum (inclusive), type `int` -* *`subset_sliced_arg`*, the subset `sliced_arg` this reduce is responsible for computing, type `T[]`, where `T` matches the type of `sliced_arg` in `reduce_sum` +* *`x_subset`*, the subset of `x` this partial sum is responsible for computing, type `T[]`, where `T` matches the type of `x` in `reduce_sum` -* *`arg1`*, first shared argument, type `T1`, matching type of `arg1` in `reduce_sum` +* *`s1`*, first shared argument, type `T1`, matching type of `s1` in `reduce_sum` -* *`arg2`*, second shared argument, type `T2`, matching type of `arg2` in `reduce_sum` +* *`s2`*, second shared argument, type `T2`, matching type of `s2` in `reduce_sum` * *`...`*, remainder of shared arguments, with types matching those in `reduce_sum` -## Parallel Map-Rect Function {#functions-map} +## Map-Rect Function {#functions-map} Stan provides a higher-order map function. This allows map-reduce functionality to be coded in Stan as described in the user's guide. diff --git a/src/stan-users-guide/parallel-computing.Rmd b/src/stan-users-guide/parallel-computing.Rmd index 705294a26..2336f7f70 100644 --- a/src/stan-users-guide/parallel-computing.Rmd +++ b/src/stan-users-guide/parallel-computing.Rmd @@ -2,82 +2,106 @@ Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and and `map_rect`. -The main differences are: +The main advantages to `reduce_sum` are: -1. `reduce_sum` requires the result of the calculation to be a scalar, while `map_rect` returns a list of vectors -2. `reduce_sum` has a more flexible interface and can accept arbitrary Stan types as arguments, `map_rect` is more restrictive on what arguments can be and how they are shaped -3. `map_rect` can parallelize work over multiple computers or a single computer, while `reduce_sum` works only on a single computer -4. `map_rect` requires work to be broken into pieces manually, while `reduce_sum` mostly automates this +1. `reduce_sum` has a more flexible argument interface, avoiding the packing and unpacking that is necessary with `map_rect`. +2. `reduce_sum` partitions the data for parallelization automatically (this is done manually in `map_rect`). +3. `reduce_sum` is easier to use. -## Reduce-Sum { #reduce-sum } +while the advantages of `map_rect` are: + +1. `map_rect` returns a list of vectors, while `reduce_sum` returns only a real. +2. `map_rect` can be parallelized across multiple computers, while `reduce_sum` can only parallelized across multiple cores. -```reduce_sum``` is a tool for parallelizing operations that can be represented as a parallel-for combined with a sum (that returns a scalar). +## Reduce-Sum { #reduce-sum } -In terms of probabilistic models, an example of this comes up when N terms in a likelihood combined multiplicatively and can be computed independently of each other (where independence here is in the computational sense, not necessarily the statistical sense). In this case, computing the log density means computing the sum of a number of terms that can each be computed separately. +```reduce_sum``` is a tool for parallelizing operations that can be represented as a sum of functions, `g: U -> real`. -```reduce_sum``` is not useful when there are dependencies between the terms. This can happen, for instance, if there were N terms in a Gaussian process likelihood. ```reduce_sum``` will not be useful for accelerating this. +For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... }```, we might compute the sum: -If for a set of input arguments, ```args0, args1, args2, ...``` and a scalar function ```f```, the log likelihood can be computed as: +```g(x1) + g(x2) + ...``` -```f(args0) + f(args1) + f(args2) + ...``` +In probabilistic modeling this comes up when there are N conditionally independent terms in a likelihood. Because of the conditional independence, these terms can be computed in parallel. If dependencies exist between the terms, then this isn't possible. For instance, in evaluating the log density of a Gaussian process ```reduce_sum``` would not be very useful. -then this calculation can be written as a reduction over the set of arguments. If this reducing function is called ```reduce```, then it would need to perform the operations: +```reduce_sum``` doesn't actually take ```g: U -> real``` as an input argument. Instead it takes ```f: U[] -> real```, where ```f``` computes the partial sum corresponding to the slice of the sequence ```x``` passed in. For instance: -```reduce({ args0, args1, args2, ... }) = f(args0) + f(args1) + f(args2) + ...``` +``` +f({ x1, x2, x3 }) = g(x1) + g(x2) + g(x3) +f({ x1 }) = g(x1) +f({ x1, x2, x3 }) = f({ x1, x2 }) + f({ x3 }) +``` -If the user can write a function like ```reduce```, then it is trivial for us to provide a function to automatically parallelize the calculations. +If the user can write a function ```f: U[] -> real``` to compute the necessary partial sums in the calculation, then we can provide a function to automatically parallelize the calculations (and this is what ```reduce_sum``` is). -Again, if the set of work is represented as a list of arguments ``{ args0, args1, args2, ... }```, then mathematically it is possible to rewrite this sum with any combination of partial-reduces. +If the set of work is represented as an array ```{ x1, x2, x3, ... }```, then mathematically it is possible to rewrite this sum with any combination of partial sums. For instance, the sum can be written: -1. ```reduce({ args0, args1, args2, ... })```, summing over all arguments, using one reduce function -2. ```reduce({ args0, ..., args(M - 1) }) + reduce({ argsM, args2, ...})```, computing the first M terms separately from the rest -3. ```reduce({ args0 }) + reduce({ args1 }) + reduce({ args2 }) + ...```, computing each term individually and summing them +1. ```f({ x1, x2, x3, ... })```, summing over all arguments, using one partial sum +2. ```f({ x1, ..., xM }) + reduce({ x(M + 1), x(M + 2), ...})```, computing the first M terms separately from the rest +3. ```f({ x1 }) + f({ x2 }) + f({ x3 }) + ...```, computing each term individually and summing them -The first function call is completely serial, the second can be parallelized over two workers, and the last can be parallelized over as many workers as there are arguments. Depending on how the list is sliced up, it is possible to parallelize this calculation over many workers. +The first form uses only one partial sum and no parallelization can be done, the second uses two partial sums and so can be parallelized over two workers, and the last can be parallelized over as many workers as there are elements in the array ```x```. Depending on how the list is sliced up, it is possible to parallelize this calculation over many workers. -```reduce_sum``` is the tool that will allow us to automatically parallelize these reduce operations (and sum them together). +```reduce_sum``` is the tool that will allow us to automatically parallelize this calculation. + +For efficiency and convenience, ```reduce_sum``` allows for additional shared arguments to be passed to every term in the sum. So for the array ```{ x1, x2, ... }``` and the shared arguments ```s1, s2, ...``` the effective sum (with individual terms) looks like: + +``` +g(x1, s1, s2, ...) + g(x2, s1, s2, ...) + g(x3, s1, s2, ...) + ... +``` -To implement this efficiently in Stan, the individual arguments are split into two types. The first are the arguments that are specific to each term in the reduction. These are called the sliced arguments (because we will slice these up to determine how to distribute the work). The second type of arguments are shared arguments, and are information that is shared in the computation of every term in the sum. +which can be written equivalently with partial sums to look like: + +``` +f({ x1, x2 }, s1, s2, ...) + f({ x3 }, s1, s2, ...) +``` + +where the particular slicing of the ```x``` array can change. Given this, the signature for reduce_sum is: ``` -real reduce_sum(F func, T[] sliced_arg, int grainsize, T1 arg1, T2 arg2, ...) +real reduce_sum(F func, T[] x, int grainsize, T1 s1, T2 s2, ...) ``` -1. ```func``` - The user-defined reduce function -2. ```sliced_arg``` - An array of any type, with each element of the array corresponding to a term of the final summation (the length of ```sliced_arg``` is the total number of terms to sum) -3. ```grainsize``` - A hint to the runtime of how many terms of the summation to compute in each reduction -4-. ```arg1, arg2, ...``` - All the arguments that are to be shared in the calculation of every term in the sum +1. ```func``` - User defined function that computes partial sums +2. ```x``` - Array to slice, each element corresponds to a term in the summation +3. ```grainsize``` - Target for size of slices +4-. ```s1, s2, ...``` - Arguments shared in every term -The user-defined reduce function is slightly different: +The user-defined partial sum functions have the signature: ``` -real func(int start, int end, T[] subset_sliced_arg, T1 arg1, T2 arg2, ...) +real func(int start, int end, T[] x_subset, T1 arg1, T2 arg2, ...) ``` -and takes the arguments: -1. ```start``` - An integer specifying the first element of the sequence of terms this reduce call is responsible for computing -2. ```end``` - An integer specifying the last element of the sequence of terms this reduce call is responsible for computing -3. ```subset_sliced_arg``` - The subset of sliced_arg for which this reduce is responsible (```sliced_arg[start:end]```) -4-. ```arg1, arg2, ...``` all the shared arguments -- passed on without modification from the reduce_sum call +and take the arguments: +1. ```start``` - An integer specifying the first term in the partial sum +2. ```end``` - An integer specifying the last term in the partial sum (inclusive) +3. ```x_subset``` - The subset of ```x``` (from ```reduce_sum```) for which this partial sum is responsible (```x[start:end]```) +4-. ```arg1, arg2, ...``` Arguments shared in every term (passed on without modification from the reduce_sum call) -The user-provided function ```func``` is expect to compute the ```start``` through ```end``` terms of the overall sum, accumulate them, and return that value. The user function is only passed the subset ```sliced_arg[start:end]``` of sliced arg (as ```subset_sliced_arg```). ```start``` and ```end``` are passed so that ```func``` can index any of the ```argM``` appropriately. The trailing arguments ```argM``` are passed without modification to every call of ```func```. +The user-provided function ```func``` is expect to compute the ```start``` through ```end``` terms of the overall sum, accumulate them, and return that value. The user function is passed the subset ```x[start:end]``` as ```x_subset```. ```start``` and ```end``` are passed so that ```func``` can index any of the tailing ```sM``` arguments as necessary. The trailing ```sM``` arguments are passed without modification to every call of ```func```. + +The ```reduce_sum``` call: + +``` +real sum = reduce_sum(func, x, grainsize, s1, s2, ...) +``` -An overall call to ```reduce_sum``` can be replaced by either: +can be replaced by either: ``` -real sum = func(1, size(sliced_arg), sliced_arg, arg1, arg2, ...) +real sum = func(1, size(x), x, s1, s2, ...) ``` -or (modulo differences due to rearrangements of summations) the code: +or the code: ``` real sum = 0.0; -for(i in 1:size(sliced_arg)) { - sum = sum + func(i, i, { sliced_arg[i] }, arg1, arg2, ...); +for(i in 1:size(x)) { + sum = sum + func(i, i, { x[i] }, s1, s2, ...); } ``` @@ -126,18 +150,18 @@ 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 partial sums of the total sum. -Using the reducer interface defined in [Reduce-Sum](#reduce-sum), such a function +Using the interface defined in [Reduce-Sum](#reduce-sum), such a function can be written like: ``` functions { - real reducer_func(int start, int end, - int[] subset_y, - vector x, - vector beta) { - return bernoulli_logit_lpmf(subset_y | beta[1] + beta[2] * x[start:end]); + real partial_sum(int start, int end, + int[] y_subset, + vector x, + vector beta) { + return bernoulli_logit_lpmf(y_subset | beta[1] + beta[2] * x[start:end]); } } ``` @@ -145,7 +169,7 @@ functions { And the likelihood statement in the model can now be written: ``` -target += reducer_fun(1, N, y, x, beta); // Sum terms 1 to N in the likelihood +target += partial_sum(1, N, y, x, beta); // Sum terms 1 to N in the likelihood ``` In this example, `y` was chosen to be sliced over because there @@ -159,14 +183,40 @@ likelihood: ``` int grainsize = 100; -target += reduce_sum(reducer_func, y, +target += reduce_sum(partial_sum, y, grainsize, x, 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. The final model looks like: + +``` +functions { + real partial_sum(int start, int end, + int[] y_subset, + vector x, + vector beta) { + return bernoulli_logit_lpmf(y_subset | beta[1] + beta[2] * x[start:end]); + } +} +data { + int N; + int y[N]; + vector[N] x; +} +parameters { + vector[2] beta; +} +model { + int grainsize = 100; + beta ~ std_normal(); + target += reduce_sum(partial_sum, y, + grainsize, + x, beta); +} +``` ### Picking the Grainsize From 3b855e82733cd602402dc5ca0f6f5100136d0846 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 27 Mar 2020 13:40:01 -0400 Subject: [PATCH 04/13] Responded to review comments for reduce_sum changes (mostly use slice instead of subset, rearranging a couple sentences, etc.) (design-doc pull request #17) --- src/stan-users-guide/_bookdown.yml | 2 +- ...llel-computing.Rmd => parallelization.Rmd} | 26 +++++++++---------- 2 files changed, 14 insertions(+), 14 deletions(-) rename src/stan-users-guide/{parallel-computing.Rmd => parallelization.Rmd} (94%) diff --git a/src/stan-users-guide/_bookdown.yml b/src/stan-users-guide/_bookdown.yml index 55c0e1b6b..938fe7660 100644 --- a/src/stan-users-guide/_bookdown.yml +++ b/src/stan-users-guide/_bookdown.yml @@ -29,7 +29,7 @@ rmd_files: [ "problematic-posteriors.Rmd", "reparameterization.Rmd", "efficiency-tuning.Rmd", - "parallel-computing.Rmd", + "parallelization.Rmd", "part-appendices.Rmd", "style-guide.Rmd", diff --git a/src/stan-users-guide/parallel-computing.Rmd b/src/stan-users-guide/parallelization.Rmd similarity index 94% rename from src/stan-users-guide/parallel-computing.Rmd rename to src/stan-users-guide/parallelization.Rmd index 2336f7f70..62b841d93 100644 --- a/src/stan-users-guide/parallel-computing.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -1,4 +1,4 @@ -# Parallel Computing {#parallel-computing.chapter} +# Parallelization {#parallelization.chapter} Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and and `map_rect`. @@ -8,14 +8,14 @@ The main advantages to `reduce_sum` are: 2. `reduce_sum` partitions the data for parallelization automatically (this is done manually in `map_rect`). 3. `reduce_sum` is easier to use. -while the advantages of `map_rect` are: +The advantages of `map_rect` are: 1. `map_rect` returns a list of vectors, while `reduce_sum` returns only a real. 2. `map_rect` can be parallelized across multiple computers, while `reduce_sum` can only parallelized across multiple cores. ## Reduce-Sum { #reduce-sum } -```reduce_sum``` is a tool for parallelizing operations that can be represented as a sum of functions, `g: U -> real`. +```reduce_sum``` parallelizes operations that can be represented as a sum of functions, `g: U -> real`. For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... }```, we might compute the sum: @@ -23,7 +23,7 @@ For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... In probabilistic modeling this comes up when there are N conditionally independent terms in a likelihood. Because of the conditional independence, these terms can be computed in parallel. If dependencies exist between the terms, then this isn't possible. For instance, in evaluating the log density of a Gaussian process ```reduce_sum``` would not be very useful. -```reduce_sum``` doesn't actually take ```g: U -> real``` as an input argument. Instead it takes ```f: U[] -> real```, where ```f``` computes the partial sum corresponding to the slice of the sequence ```x``` passed in. For instance: +```reduce_sum``` takes a function ```f: U[] -> real```, where ```f``` computes the partial sum corresponding to the slice of the sequence ```x``` passed in. For instance: ``` f({ x1, x2, x3 }) = g(x1) + g(x2) + g(x3) @@ -31,7 +31,7 @@ f({ x1 }) = g(x1) f({ x1, x2, x3 }) = f({ x1, x2 }) + f({ x3 }) ``` -If the user can write a function ```f: U[] -> real``` to compute the necessary partial sums in the calculation, then we can provide a function to automatically parallelize the calculations (and this is what ```reduce_sum``` is). +If the user can write a function ```f: U[] -> real``` to compute the necessary partial sums in the calculation, then ```reduce_sum``` can automatically parallelize the calculations. If the set of work is represented as an array ```{ x1, x2, x3, ... }```, then mathematically it is possible to rewrite this sum with any combination of partial sums. @@ -73,16 +73,16 @@ real reduce_sum(F func, T[] x, int grainsize, T1 s1, T2 s2, ...) The user-defined partial sum functions have the signature: ``` -real func(int start, int end, T[] x_subset, T1 arg1, T2 arg2, ...) +real func(int start, int end, T[] x_slice, T1 arg1, T2 arg2, ...) ``` and take the arguments: 1. ```start``` - An integer specifying the first term in the partial sum 2. ```end``` - An integer specifying the last term in the partial sum (inclusive) -3. ```x_subset``` - The subset of ```x``` (from ```reduce_sum```) for which this partial sum is responsible (```x[start:end]```) +3. ```x_slice``` - The subset of ```x``` (from ```reduce_sum```) for which this partial sum is responsible (```x[start:end]```) 4-. ```arg1, arg2, ...``` Arguments shared in every term (passed on without modification from the reduce_sum call) -The user-provided function ```func``` is expect to compute the ```start``` through ```end``` terms of the overall sum, accumulate them, and return that value. The user function is passed the subset ```x[start:end]``` as ```x_subset```. ```start``` and ```end``` are passed so that ```func``` can index any of the tailing ```sM``` arguments as necessary. The trailing ```sM``` arguments are passed without modification to every call of ```func```. +The user-provided function ```func``` is expect to compute the ```start``` through ```end``` terms of the overall sum, accumulate them, and return that value. The user function is passed the subset ```x[start:end]``` as ```x_slice```. ```start``` and ```end``` are passed so that ```func``` can index any of the tailing ```sM``` arguments as necessary. The trailing ```sM``` arguments are passed without modification to every call of ```func```. The ```reduce_sum``` call: @@ -158,10 +158,10 @@ can be written like: ``` functions { real partial_sum(int start, int end, - int[] y_subset, + int[] y_slice, vector x, vector beta) { - return bernoulli_logit_lpmf(y_subset | beta[1] + beta[2] * x[start:end]); + return bernoulli_logit_lpmf(y_slice | beta[1] + beta[2] * x[start:end]); } } ``` @@ -195,10 +195,10 @@ be estimated automatically. The final model looks like: ``` functions { real partial_sum(int start, int end, - int[] y_subset, + int[] y_slice, vector x, vector beta) { - return bernoulli_logit_lpmf(y_subset | beta[1] + beta[2] * x[start:end]); + return bernoulli_logit_lpmf(y_slice | beta[1] + beta[2] * x[start:end]); } } data { @@ -221,7 +221,7 @@ model { ### Picking the Grainsize The `grainsize` is a recommendation on how large each piece of parallel work is -(how many terms it contains). If zero, it will be chosen automatically, but it +(how many terms it contains). If one, it will be chosen automatically, but it is probably best to choose this manually for each model. To figure out an appropriate grainsize, think about how many terms are in the summation From 450ec3fff227499a12bc87e440ff56bae8fe5e78 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Tue, 31 Mar 2020 14:45:49 +0200 Subject: [PATCH 05/13] refine reduce_sum function reference --- .../higher-order_functions.Rmd | 53 ++++++++++--------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index b2c9d3845..de6ea97fc 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -385,47 +385,48 @@ The gradients of the integral are computed in accordance with the Leibniz integr ## Reduce-Sum Function {#functions-reduce} -Stan provides a higher-order `reduce_sum` function for parallelizing operations that can be represented as a sum of a function, `g: U -> real`, evaluated at each element of a list of type `U[]`, `{ x1, x2, ... }`. That is: +Stan provides a higher-order `reduce_sum` function which maps +evaluation of a function `g: U -> real` to a list of type `U[]`, `{ +x1, x2, ... }`, and performs as reduction operation a sum over the +results. That is: ```g(x1) + g(x2) + ...``` -`reduce_sum` doesn't work on `g` itself, but takes a partial sum function, `f: U[] -> real`, where: +`reduce_sum` doesn't work with the element-wise evaluated +function `g` itself, but instead works through evaluating partial +sums, `f: U[] -> real`, where: -```f({ x1 }) = g(x1)``` -```f({ x1, x2 }) = g(x1) + g(x2)``` -```f({ x1, x2, ... }) = g(x1) + g(x2) + ...``` +``` +f({ x1 }) = g(x1) +f({ x1, x2 }) = g(x1) + g(x2) +f({ x1, x2, ... }) = g(x1) + g(x2) + ... +``` -### The Reduce-sum Function +### Specifying the Reduce-sum Function -The `reduce_sum` function takes a partial sum function, an array argument x -(with one for each term in the sum), a recommended grainsize, and a set of shared arguments and -parallelizes the resultant sum. +The `reduce_sum` function takes a partial sum function `f`, an array argument `x` +(with one array element for each term in the sum), a recommended +`grainsize`, and a set of shared arguments. This representation allows +to parallelize the resultant sum. \index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] x, int grainsize, T1 s1, T2 s2, ...): real}|hyperpage} `real` **`reduce_sum`**`(F f, T[] x, int grainsize, T1 s1, T2 s2, ...)`
\newline -Return the equivalent of `f(1, size(x), x, s1, s2, ...)`, but compute -the result in parallel by breaking the array `x` into pieces and computing each piece -in a different thread. `s1, s2, ...` are shared between all terms in the sum. - -* *`f`*: function literal referring to a function specifying the partial sum operation with signature `(int, int, T[], T1, T2, ...):real` -The arguments represent - + (1) the index of the first term of the partial sum, - + (2) the index of the last term of the partial sum, - + (3) the subset `x` this reduce is responsible for computing, - + (4) first shared argument, - + (5) second shared argument, - + ... the rest of the shared arguments. +Returns the equivalent of `f(1, size(x), x, s1, s2, ...)`, but computes +the result in parallel by breaking the array `x` into independent +partial sums. `s1, s2, ...` are shared between all terms in the sum. +* *`f`*: function literal referring to a function specifying the +partial sum operation. Refer to the [partial sum function](#functions-partial-sum). * *`x`*: array of `T`, one for each term of the reduction, `T` can be any type, -* *`grainsize`*: recommented number of terms in each reduce call, set to one to estimate automatically, type `int`, +* *`grainsize`*: recommended number of terms in each reduce call, set to one to estimate automatically, type `int`, * *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type * *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, * *`...`*: remainder of shared arguments, each of which can be any type. -### The Partial-sum Function +### The Partial-sum Function {#functions-partial-sum} The partial sum function must have the following signature where the types `T`, and the types of all the shared arguments (`T1`, `T2`, ...) match those of the original @@ -435,14 +436,14 @@ types of all the shared arguments (`T1`, `T2`, ...) match those of the original (int start, int end, T[] x_subset, T1 s1, T2 s2, ...):real ``` -The reduce function returns the sum of the `start` to `end` terms (inclusive) of the overall -calculations. The arguments to the reduce function are: +The partial sum function returns the sum of the `start` to `end` terms (inclusive) of the overall +calculations. The arguments to the partial sum function are: * *`start`*, the index of the first term of the partial sum, type `int` * *`end`*, the index of the last term of the partial sum (inclusive), type `int` -* *`x_subset`*, the subset of `x` this partial sum is responsible for computing, type `T[]`, where `T` matches the type of `x` in `reduce_sum` +* *`x_subset`*, the subset of `x` a given partial sum is responsible for computing, type `T[]`, where `T` matches the type of `x` in `reduce_sum` * *`s1`*, first shared argument, type `T1`, matching type of `s1` in `reduce_sum` From 8fca9292654618218debe78f01592aa9b5eaec97 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Tue, 31 Mar 2020 16:47:06 +0200 Subject: [PATCH 06/13] refine stan user manual --- src/stan-users-guide/parallelization.Rmd | 135 ++++++++++++++--------- 1 file changed, 80 insertions(+), 55 deletions(-) diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index 62b841d93..499346f6b 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -1,51 +1,62 @@ # Parallelization {#parallelization.chapter} -Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and and `map_rect`. +Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and `map_rect`. -The main advantages to `reduce_sum` are: +The advantages of `reduce_sum` are: 1. `reduce_sum` has a more flexible argument interface, avoiding the packing and unpacking that is necessary with `map_rect`. -2. `reduce_sum` partitions the data for parallelization automatically (this is done manually in `map_rect`). +2. `reduce_sum` partitions data for parallelization automatically (this is done manually in `map_rect`). 3. `reduce_sum` is easier to use. The advantages of `map_rect` are: 1. `map_rect` returns a list of vectors, while `reduce_sum` returns only a real. -2. `map_rect` can be parallelized across multiple computers, while `reduce_sum` can only parallelized across multiple cores. +2. `map_rect` can be parallelized across multiple cores and multiple + computers, while `reduce_sum` can only parallelized across multiple + cores on a single machine. + +The actual speedup gained from using these functions will depend on +many details. It is strongly recommended to only parallelize the +computationally most expensive operations in a Stan +program. Oftentimes this is the evaluation of the log likelihood for +the observed data. Since only portions of a Stan program will run in +parallel, the maximal speedup one can achieve is capped, a phenomen +described by [Amdahl's +law](https://en.wikipedia.org/wiki/Amdahl's_law). ## Reduce-Sum { #reduce-sum } -```reduce_sum``` parallelizes operations that can be represented as a sum of functions, `g: U -> real`. - -For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... }```, we might compute the sum: +```reduce_sum``` maps evaluation of a function `g: U -> real` to a list of type `U[]`, `{ +x1, x2, ... }`, and performs as reduction operation a sum over the +results. For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... }```, we might compute the sum: ```g(x1) + g(x2) + ...``` In probabilistic modeling this comes up when there are N conditionally independent terms in a likelihood. Because of the conditional independence, these terms can be computed in parallel. If dependencies exist between the terms, then this isn't possible. For instance, in evaluating the log density of a Gaussian process ```reduce_sum``` would not be very useful. -```reduce_sum``` takes a function ```f: U[] -> real```, where ```f``` computes the partial sum corresponding to the slice of the sequence ```x``` passed in. For instance: +```reduce_sum``` requires the partial sum function ```f: U[] -> +real```, where ```f``` computes the partial sum corresponding to the +slice of the sequence ```x``` passed in. ```reduce_sum``` +exploits the associativity of the sum operation as it holds that: ``` -f({ x1, x2, x3 }) = g(x1) + g(x2) + g(x3) -f({ x1 }) = g(x1) -f({ x1, x2, x3 }) = f({ x1, x2 }) + f({ x3 }) +g(x1) + g(x2) + g(x3) = f({ x1, x2, x3 }) + = f({ x1, x2 }) + f({ x3 }) + = f({ x1 }) + f({ x2, x3 }) + = f({ x1 }) + f({ x2 }) + f({ x3 }) ``` -If the user can write a function ```f: U[] -> real``` to compute the necessary partial sums in the calculation, then ```reduce_sum``` can automatically parallelize the calculations. - -If the set of work is represented as an array ```{ x1, x2, x3, ... }```, then mathematically it is possible to rewrite this sum with any combination of partial sums. - -For instance, the sum can be written: - -1. ```f({ x1, x2, x3, ... })```, summing over all arguments, using one partial sum -2. ```f({ x1, ..., xM }) + reduce({ x(M + 1), x(M + 2), ...})```, computing the first M terms separately from the rest -3. ```f({ x1 }) + f({ x2 }) + f({ x3 }) + ...```, computing each term individually and summing them - -The first form uses only one partial sum and no parallelization can be done, the second uses two partial sums and so can be parallelized over two workers, and the last can be parallelized over as many workers as there are elements in the array ```x```. Depending on how the list is sliced up, it is possible to parallelize this calculation over many workers. +If the user can write a function ```f: U[] -> real``` to compute the +necessary partial sums of the calculation, then ```reduce_sum``` can +automatically parallelize the calculations. The exact partitioning +into partial sums is arbitrary as these are mathematical equivalent to +one another. As the partitioning is flexible, it is be adapted to the +available ressources (number of concurrent threads) given to Stan. -```reduce_sum``` is the tool that will allow us to automatically parallelize this calculation. - -For efficiency and convenience, ```reduce_sum``` allows for additional shared arguments to be passed to every term in the sum. So for the array ```{ x1, x2, ... }``` and the shared arguments ```s1, s2, ...``` the effective sum (with individual terms) looks like: +For efficiency and convenience, ```reduce_sum``` allows for additional +shared arguments to be passed to every term in the sum. So for the +array ```{ x1, x2, ... }``` and the shared arguments ```s1, s2, ...``` +the effective sum (with individual terms) looks like: ``` g(x1, s1, s2, ...) + g(x2, s1, s2, ...) + g(x3, s1, s2, ...) + ... @@ -59,41 +70,49 @@ f({ x1, x2 }, s1, s2, ...) + f({ x3 }, s1, s2, ...) where the particular slicing of the ```x``` array can change. -Given this, the signature for reduce_sum is: +Given this, the signature for ```reduce_sum``` is: ``` -real reduce_sum(F func, T[] x, int grainsize, T1 s1, T2 s2, ...) +real reduce_sum(F f, T[] x, int grainsize, T1 s1, T2 s2, ...) ``` -1. ```func``` - User defined function that computes partial sums +1. ```f``` - User defined function that computes partial sums 2. ```x``` - Array to slice, each element corresponds to a term in the summation 3. ```grainsize``` - Target for size of slices -4-. ```s1, s2, ...``` - Arguments shared in every term +4. ```s1, s2, ...``` - Arguments shared in every term The user-defined partial sum functions have the signature: ``` -real func(int start, int end, T[] x_slice, T1 arg1, T2 arg2, ...) +real f(int start, int end, T[] x_slice, T1 s1, T2 s2, ...) ``` and take the arguments: + 1. ```start``` - An integer specifying the first term in the partial sum 2. ```end``` - An integer specifying the last term in the partial sum (inclusive) -3. ```x_slice``` - The subset of ```x``` (from ```reduce_sum```) for which this partial sum is responsible (```x[start:end]```) -4-. ```arg1, arg2, ...``` Arguments shared in every term (passed on without modification from the reduce_sum call) - -The user-provided function ```func``` is expect to compute the ```start``` through ```end``` terms of the overall sum, accumulate them, and return that value. The user function is passed the subset ```x[start:end]``` as ```x_slice```. ```start``` and ```end``` are passed so that ```func``` can index any of the tailing ```sM``` arguments as necessary. The trailing ```sM``` arguments are passed without modification to every call of ```func```. +3. ```x_slice``` - The subset of ```x``` (from ```reduce_sum```) for +which this partial sum is responsible (```x_slice = x[start:end]```) +4. ```s1, s2, ...``` - Arguments shared in every term (passed on without modification from the ```reduce_sum``` call) + +The user-provided function ```f``` is expected to compute the partial +sum with the terms ```start``` through ```end``` of the overall +sum. The user function is passed the subset ```x[start:end]``` as +```x_slice```. ```start``` and ```end``` are passed so that ```f``` +can index any of the tailing ```sM``` arguments as necessary. The +trailing ```sM``` arguments are passed without modification to every +call of ```f```. The ```reduce_sum``` call: ``` -real sum = reduce_sum(func, x, grainsize, s1, s2, ...) +real sum = reduce_sum(f, x, grainsize, s1, s2, ...); ``` can be replaced by either: ``` -real sum = func(1, size(x), x, s1, s2, ...) +real sum = f(1, size(x), x, s1, s2, ...); ``` or the code: @@ -101,14 +120,14 @@ or the code: ``` real sum = 0.0; for(i in 1:size(x)) { - sum = sum + func(i, i, { x[i] }, s1, s2, ...); + sum += f(i, i, { x[i] }, s1, s2, ...); } ``` ### Example: Logistic Regression -Logistic Regression is a useful example to clarify both the syntax -and semantics of `reduce_sum` and how it can be used to speed up a typical +Logistic regression is a useful example to clarify both the syntax +and semantics of ```reduce_sum``` and how it can be used to speed up a typical model. A basic logistic regression can be coded in Stan as: @@ -147,9 +166,9 @@ 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. +```reduce_sum``` is useful. -To use `reduce_sum`, a function must be written that can be used to compute +To use ```reduce_sum```, a function must be written that can be used to compute arbitrary partial sums of the total sum. Using the interface defined in [Reduce-Sum](#reduce-sum), such a function @@ -169,14 +188,14 @@ functions { And the likelihood statement in the model can now be written: ``` -target += partial_sum(1, N, y, x, beta); // Sum terms 1 to N in the likelihood +target += partial_sum(1, N, y, x, beta); // Sum terms 1 to N of the likelihood ``` In this example, `y` was chosen to be sliced over because there is one term in the summation per value of `y`. Technically `x` would have worked as well. Use whatever conceptually makes the most sense. -Because `x` is a shared argument, it is subset manually with `start:end`. +Because `x` is a shared argument, it is subset accordingly with `start:end`. With this function, `reduce_sum` can be used to automatically parallelize the likelihood: @@ -220,25 +239,31 @@ model { ### Picking the Grainsize -The `grainsize` is a recommendation on how large each piece of parallel work is -(how many terms it contains). If one, it will be chosen automatically, but it -is probably best to choose this manually for each model. - -To figure out an appropriate grainsize, think about how many terms are in the summation -and on how many cores the model should run. If there are `N` terms and `M` cores, -run a quick test model with `grainsize` set roughly to `N / M`. Record the time, cut the -grainsize in half, and run the test again. Repeat this iteratively until the model runtime -begins to increase. This is a suitable grainsize for the model, because this ensures the -caculations can be carried out with the most parallelism without losing too much efficiency. +The `grainsize` is a recommendation on how large each piece of +parallel work is (how many terms it contains). It is recommended to +choose one as a starting point which will select an appropiate value +automatically. + +From empirical experience, the automatic grainsize determination works +well and no further tuning is required in most cases. In order to +figure out an optimal grainsize, think about how many terms are in the +summation and on how many cores the model should run. If there are `N` +terms and `M` cores, run a quick test model with `grainsize` set +roughly to `N / M`. Record the time, cut the grainsize in half, and +run the test again. Repeat this iteratively until the model runtime +begins to increase. This is a suitable grainsize for the model, +because this ensures the caculations can be carried out with the most +parallelism without losing too much efficiency. For instance, in a model with `N=10000` and `M = 4`, start with `grainsize = 25000`, and sequentially try `grainsize = 12500`, `grainsize = 6250`, etc. It is important to repeat this process until performance gets worse! It is possible -after many halvings nothing happens, but there might still be a small grainsize that performs better. +after many halvings nothing happens, but there might still be a smaller grainsize that performs better. Even if a sum has many tens of thousands of terms, depending on the internal calculations, a `grainsize` of thirty or forty or smaller might be the best, and it is difficult to predict this behavior. -Without doing these halvings until performance actually gets worse, it is easy to miss this. +Without doing these halvings until performance actually gets worse, it +is easy to miss this. ## Map-Rect From f4e1bce5254332ac74ea469c7bbbaa223f487b96 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 1 Apr 2020 15:54:45 +0200 Subject: [PATCH 07/13] add to reference the reduce_sum_static function ref --- .../higher-order_functions.Rmd | 46 +++++++++++++------ 1 file changed, 32 insertions(+), 14 deletions(-) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index de6ea97fc..a4decad53 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -385,16 +385,16 @@ The gradients of the integral are computed in accordance with the Leibniz integr ## Reduce-Sum Function {#functions-reduce} -Stan provides a higher-order `reduce_sum` function which maps -evaluation of a function `g: U -> real` to a list of type `U[]`, `{ -x1, x2, ... }`, and performs as reduction operation a sum over the -results. That is: +Stan provides a higher-order reduce function for summation. A function +`g: U -> real`, which returns a scalar, is mapped to every element of +a list of type `U[]`, `{ x1, x2, ... }` and all the results are +accumulated, ```g(x1) + g(x2) + ...``` -`reduce_sum` doesn't work with the element-wise evaluated -function `g` itself, but instead works through evaluating partial -sums, `f: U[] -> real`, where: +For efficiency reasons the reduce function doesn't work with the +element-wise evaluated function `g` itself, but instead works through +evaluating partial sums, `f: U[] -> real`, where: ``` f({ x1 }) = g(x1) @@ -402,9 +402,24 @@ f({ x1, x2 }) = g(x1) + g(x2) f({ x1, x2, ... }) = g(x1) + g(x2) + ... ``` +Mathematically the summation reduction is associative and forming +arbitrary partial sums in an aribitrary order will not change the +result. However, floating point numerics on computers only have +a limited precision such that associativity does not hold +exactly. This implies that the order of summation determines the exact +numerical result. For this reason, the higher-order reduce function is +available in two variants: + +* `reduce_sum`: Automatically forms partial sums resulting usually in good + performance without further tuning. +* `reduce_sum_static`: Creates for the same input always the same +call graph resulting in stable numerical evaluation. This version +requires setting a tuning parameter which controls the maximal size of partial +sums formed. + ### Specifying the Reduce-sum Function -The `reduce_sum` function takes a partial sum function `f`, an array argument `x` +The higher-order reduce function takes a partial sum function `f`, an array argument `x` (with one array element for each term in the sum), a recommended `grainsize`, and a set of shared arguments. This representation allows to parallelize the resultant sum. @@ -413,6 +428,7 @@ to parallelize the resultant sum. \index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] x, int grainsize, T1 s1, T2 s2, ...): real}|hyperpage} `real` **`reduce_sum`**`(F f, T[] x, int grainsize, T1 s1, T2 s2, ...)`
\newline +`real` **`reduce_sum_static`**`(F f, T[] x, int grainsize, T1 s1, T2 s2, ...)`
\newline Returns the equivalent of `f(1, size(x), x, s1, s2, ...)`, but computes the result in parallel by breaking the array `x` into independent @@ -421,7 +437,9 @@ partial sums. `s1, s2, ...` are shared between all terms in the sum. * *`f`*: function literal referring to a function specifying the partial sum operation. Refer to the [partial sum function](#functions-partial-sum). * *`x`*: array of `T`, one for each term of the reduction, `T` can be any type, -* *`grainsize`*: recommended number of terms in each reduce call, set to one to estimate automatically, type `int`, +* *`grainsize`*: recommended number of terms in each reduce call, set +to one to estimate automatically for `reduce_sum` while for +`reduce_sum_static` this determines the maximal size of the partial sums, type `int`, * *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type * *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, * *`...`*: remainder of shared arguments, each of which can be any type. @@ -430,7 +448,7 @@ partial sum operation. Refer to the [partial sum function](#functions-partial-su The partial sum function must have the following signature where the types `T`, and the types of all the shared arguments (`T1`, `T2`, ...) match those of the original -`reduce_sum` call. +`reduce_sum` (`reduce_sum_static`) call. ``` (int start, int end, T[] x_subset, T1 s1, T2 s2, ...):real @@ -443,13 +461,13 @@ calculations. The arguments to the partial sum function are: * *`end`*, the index of the last term of the partial sum (inclusive), type `int` -* *`x_subset`*, the subset of `x` a given partial sum is responsible for computing, type `T[]`, where `T` matches the type of `x` in `reduce_sum` +* *`x_subset`*, the subset of `x` a given partial sum is responsible for computing, type `T[]`, where `T` matches the type of `x` in `reduce_sum` (`reduce_sum_static`) -* *`s1`*, first shared argument, type `T1`, matching type of `s1` in `reduce_sum` +* *`s1`*, first shared argument, type `T1`, matching type of `s1` in `reduce_sum` (`reduce_sum_static`) -* *`s2`*, second shared argument, type `T2`, matching type of `s2` in `reduce_sum` +* *`s2`*, second shared argument, type `T2`, matching type of `s2` in `reduce_sum` (`reduce_sum_static`) -* *`...`*, remainder of shared arguments, with types matching those in `reduce_sum` +* *`...`*, remainder of shared arguments, with types matching those in `reduce_sum` (`reduce_sum_static`) ## Map-Rect Function {#functions-map} From df152e4694b430692557fcc908c2ef62c5cc317e Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 1 Apr 2020 16:57:30 +0200 Subject: [PATCH 08/13] add static version to user guide --- src/stan-users-guide/parallelization.Rmd | 129 +++++++++++++++-------- 1 file changed, 83 insertions(+), 46 deletions(-) diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index 499346f6b..2692340d9 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -1,18 +1,20 @@ # Parallelization {#parallelization.chapter} -Stan has two mechanisms for parallelizing calculations used in a model: `reduce_sum` and `map_rect`. +Stan has two mechanisms for parallelizing calculations used in a +model: reduce with summation and rectangular map. -The advantages of `reduce_sum` are: +The advantages of reduce with summation are: -1. `reduce_sum` has a more flexible argument interface, avoiding the packing and unpacking that is necessary with `map_rect`. -2. `reduce_sum` partitions data for parallelization automatically (this is done manually in `map_rect`). -3. `reduce_sum` is easier to use. +1. More flexible argument interface, avoiding the packing and + unpacking that is necessary with rectanguar map. +2. Partitions data for parallelization automatically (this is done manually in rectanguar map). +3. Is easier to use. -The advantages of `map_rect` are: +The advantages of rectangular map are: -1. `map_rect` returns a list of vectors, while `reduce_sum` returns only a real. -2. `map_rect` can be parallelized across multiple cores and multiple - computers, while `reduce_sum` can only parallelized across multiple +1. Returns a list of vectors, while the reduce summation returns only a scalar. +2. Can be parallelized across multiple cores and multiple + computers, while reduce summation can only parallelized across multiple cores on a single machine. The actual speedup gained from using these functions will depend on @@ -26,18 +28,26 @@ law](https://en.wikipedia.org/wiki/Amdahl's_law). ## Reduce-Sum { #reduce-sum } -```reduce_sum``` maps evaluation of a function `g: U -> real` to a list of type `U[]`, `{ -x1, x2, ... }`, and performs as reduction operation a sum over the -results. For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... }```, we might compute the sum: +The higher-order reduce with summation facility maps evaluation of a +function `g: U -> real`, which returns a scalar, to a list of type +`U[]`, `{ x1, x2, ... }`, and performs as reduction operation a sum +over the results. For instance, for a sequence of ```x``` values of +type ```U```, ```{ x1, x2, ... }```, we might compute the sum: ```g(x1) + g(x2) + ...``` -In probabilistic modeling this comes up when there are N conditionally independent terms in a likelihood. Because of the conditional independence, these terms can be computed in parallel. If dependencies exist between the terms, then this isn't possible. For instance, in evaluating the log density of a Gaussian process ```reduce_sum``` would not be very useful. +In probabilistic modeling this comes up when there are $N$ +conditionally independent terms in a likelihood. Because of the +conditional independence, these terms can be computed in parallel. If +dependencies exist between the terms, then this isn't possible. For +instance, in evaluating the log density of a Gaussian process then +summation of independent terms isn't applicable. -```reduce_sum``` requires the partial sum function ```f: U[] -> -real```, where ```f``` computes the partial sum corresponding to the -slice of the sequence ```x``` passed in. ```reduce_sum``` -exploits the associativity of the sum operation as it holds that: +For efficiency reasons the reduce function doesn’t work with the +element-wise evaluated function `g`, but instead requires the partial +sum function ```f: U[] -> real```, where ```f``` computes the partial +sum corresponding to a slice of the sequence ```x``` passed in. Due to the +the associativity of the sum reduction it holds that: ``` g(x1) + g(x2) + g(x3) = f({ x1, x2, x3 }) @@ -46,15 +56,31 @@ g(x1) + g(x2) + g(x3) = f({ x1, x2, x3 }) = f({ x1 }) + f({ x2 }) + f({ x3 }) ``` -If the user can write a function ```f: U[] -> real``` to compute the -necessary partial sums of the calculation, then ```reduce_sum``` can -automatically parallelize the calculations. The exact partitioning -into partial sums is arbitrary as these are mathematical equivalent to -one another. As the partitioning is flexible, it is be adapted to the -available ressources (number of concurrent threads) given to Stan. - -For efficiency and convenience, ```reduce_sum``` allows for additional -shared arguments to be passed to every term in the sum. So for the +With the partial sum function ```f: U[] -> real``` reduction of a +large number of terms can be evaluated in parallel automatically, since the +overall sum can be partitioned into arbitrary smaller partial +sums. The exact partitioning into the partial sums is not under the +control of the user. However, since the exact numerical result will +depend on the order of summation, Stan provides two versions of the +reduce summation facility: + +* `reduce_sum`: Automatically forms partial sums resulting usually in good + performance without further tuning. +* `reduce_sum_static`: Creates for the same input always the same +call graph resulting in stable numerical evaluation. This version +requires setting a sensible tuning parameter for good performance. + +The tuning parameter is the so-called `grainsize`. For the +`reduce_sum` version the `grainsize` is merely a suggested partial sum +size while for the `reduce_sum_static` version the `grainsize` +specifies the maximal partial sum size. While for `reduce_sum` a +`grainsize` of one commonly leads to good performance already (since +automatic aggregation is performed), the `reduce_sum_static` variant +requires setting a sensible `grainsize` for good performance as +explained in [more detail below](#reduce-sum-grainsize). + +For efficiency and convenience additional +shared arguments can be passed to every term in the sum. So for the array ```{ x1, x2, ... }``` and the shared arguments ```s1, s2, ...``` the effective sum (with individual terms) looks like: @@ -70,10 +96,11 @@ f({ x1, x2 }, s1, s2, ...) + f({ x3 }, s1, s2, ...) where the particular slicing of the ```x``` array can change. -Given this, the signature for ```reduce_sum``` is: +Given this, the signatures are: ``` real reduce_sum(F f, T[] x, int grainsize, T1 s1, T2 s2, ...) +real reduce_sum_static(F f, T[] x, int grainsize, T1 s1, T2 s2, ...) ``` 1. ```f``` - User defined function that computes partial sums @@ -91,9 +118,10 @@ and take the arguments: 1. ```start``` - An integer specifying the first term in the partial sum 2. ```end``` - An integer specifying the last term in the partial sum (inclusive) -3. ```x_slice``` - The subset of ```x``` (from ```reduce_sum```) for +3. ```x_slice``` - The subset of ```x``` (from ```reduce_sum``` / `reduce_sum_static`) for which this partial sum is responsible (```x_slice = x[start:end]```) -4. ```s1, s2, ...``` - Arguments shared in every term (passed on without modification from the ```reduce_sum``` call) +4. ```s1, s2, ...``` - Arguments shared in every term (passed on +without modification from the ```reduce_sum``` / `reduce_sum_static` call) The user-provided function ```f``` is expected to compute the partial sum with the terms ```start``` through ```end``` of the overall @@ -103,7 +131,7 @@ can index any of the tailing ```sM``` arguments as necessary. The trailing ```sM``` arguments are passed without modification to every call of ```f```. -The ```reduce_sum``` call: +A ```reduce_sum``` (or `reduce_sum_static`) call: ``` real sum = reduce_sum(f, x, grainsize, s1, s2, ...); @@ -127,7 +155,7 @@ for(i in 1:size(x)) { ### Example: Logistic Regression Logistic regression is a useful example to clarify both the syntax -and semantics of ```reduce_sum``` and how it can be used to speed up a typical +and semantics of reduce summation and how it can be used to speed up a typical model. A basic logistic regression can be coded in Stan as: @@ -150,7 +178,7 @@ model { In this model predictions are made about the `N` outputs `y` using the covariate `x`. The intercept and slope of the linear equation are to be estimated. -The key point to getting this calculation into `reduce_sum`, is recognizing that +The key point to getting this calculation to use reduce summation, is recognizing that the statement: ``` @@ -166,9 +194,9 @@ 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. +reduce summation is useful. -To use ```reduce_sum```, a function must be written that can be used to compute +To use the reduce summation, a function must be written that can be used to compute arbitrary partial sums of the total sum. Using the interface defined in [Reduce-Sum](#reduce-sum), such a function @@ -197,7 +225,7 @@ worked as well. Use whatever conceptually makes the most sense. Because `x` is a shared argument, it is subset accordingly with `start:end`. -With this function, `reduce_sum` can be used to automatically parallelize the +With this function, reduce summation can be used to automatically parallelize the likelihood: ``` @@ -207,7 +235,7 @@ target += reduce_sum(partial_sum, y, x, beta); ``` -`reduce_sum` automatically breaks the sum into roughly `grainsize` sized pieces +The reduce summation facility automatically breaks the sum into roughly `grainsize` sized pieces and computes them in parallel. `grainsize = 1` specifies that the grainsize should be estimated automatically. The final model looks like: @@ -237,12 +265,19 @@ model { } ``` -### Picking the Grainsize +### Picking the Grainsize {#reduce-sum-grainsize} The `grainsize` is a recommendation on how large each piece of -parallel work is (how many terms it contains). It is recommended to -choose one as a starting point which will select an appropiate value -automatically. +parallel work is (how many terms it contains). When using the +non-static version, it is recommended to choose one as a starting +point as automatic aggregation of partial sums are performed. However, +for the static version the `grainsize` defines the maximal size of the +partial sums, e.g. the static variant will split the input sequence +until all partial sums are just smaller than `grainsize`. Therefore, +for the static version it is more important to select a sensible +value. The rational for choosing a sensible `grainsize` is based on +balancing the overhead implied by creating many small tasks versus +creating fewer large tasks which limits the potential parallelism. From empirical experience, the automatic grainsize determination works well and no further tuning is required in most cases. In order to @@ -258,12 +293,14 @@ parallelism without losing too much efficiency. For instance, in a model with `N=10000` and `M = 4`, start with `grainsize = 25000`, and sequentially try `grainsize = 12500`, `grainsize = 6250`, etc. -It is important to repeat this process until performance gets worse! It is possible -after many halvings nothing happens, but there might still be a smaller grainsize that performs better. -Even if a sum has many tens of thousands of terms, depending on the internal calculations, a `grainsize` -of thirty or forty or smaller might be the best, and it is difficult to predict this behavior. -Without doing these halvings until performance actually gets worse, it -is easy to miss this. +It is important to repeat this process until performance gets worse! +It is possible after many halvings nothing happens, but there might +still be a smaller grainsize that performs better. Even if a sum has +many tens of thousands of terms, depending on the internal +calculations, a `grainsize` of thirty or forty or smaller might be the +best, and it is difficult to predict this behavior. Without doing +these halvings until performance actually gets worse, it is easy to +miss this. ## Map-Rect From 88311633c43ea0b23600d887fb431cecf8fc0f66 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Tue, 7 Apr 2020 09:26:01 +0200 Subject: [PATCH 09/13] address review comments --- .../higher-order_functions.Rmd | 22 +++++----- src/stan-users-guide/parallelization.Rmd | 42 ++++++++----------- 2 files changed, 29 insertions(+), 35 deletions(-) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index a4decad53..ea1c1d612 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -11,7 +11,7 @@ cat(' * Algebraic Equation Solver\ cat(' * Ordinary Differential Equation (ODE) Solvers\n') cat(' * 1D Integrator\n') cat(' * Reduce-Sum\n') -cat(' * Higher-Order Map\n') +cat(' * Map-Rect\n') } ``` @@ -386,11 +386,11 @@ The gradients of the integral are computed in accordance with the Leibniz integr ## Reduce-Sum Function {#functions-reduce} Stan provides a higher-order reduce function for summation. A function -`g: U -> real`, which returns a scalar, is mapped to every element of -a list of type `U[]`, `{ x1, x2, ... }` and all the results are +which returns a scalar `g: U -> real` is mapped to every element of a +list of type `U[]`, `{ x1, x2, ... }` and all the results are accumulated, -```g(x1) + g(x2) + ...``` +`g(x1) + g(x2) + ...` For efficiency reasons the reduce function doesn't work with the element-wise evaluated function `g` itself, but instead works through @@ -410,12 +410,12 @@ exactly. This implies that the order of summation determines the exact numerical result. For this reason, the higher-order reduce function is available in two variants: -* `reduce_sum`: Automatically forms partial sums resulting usually in good - performance without further tuning. -* `reduce_sum_static`: Creates for the same input always the same -call graph resulting in stable numerical evaluation. This version -requires setting a tuning parameter which controls the maximal size of partial -sums formed. +* `reduce_sum`: Compute partial sums automatically. This usually + results in good performance without further tuning. +* `reduce_sum_static`: For the same input, always create the same call +graph. This results in stable numerical evaluation. This version +requires setting a tuning parameter which controls the maximal size of +partial sums formed. ### Specifying the Reduce-sum Function @@ -438,7 +438,7 @@ partial sums. `s1, s2, ...` are shared between all terms in the sum. partial sum operation. Refer to the [partial sum function](#functions-partial-sum). * *`x`*: array of `T`, one for each term of the reduction, `T` can be any type, * *`grainsize`*: recommended number of terms in each reduce call, set -to one to estimate automatically for `reduce_sum` while for +to 1 to estimate automatically for `reduce_sum` while for `reduce_sum_static` this determines the maximal size of the partial sums, type `int`, * *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type * *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index 2692340d9..f2372f83d 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -34,7 +34,7 @@ function `g: U -> real`, which returns a scalar, to a list of type over the results. For instance, for a sequence of ```x``` values of type ```U```, ```{ x1, x2, ... }```, we might compute the sum: -```g(x1) + g(x2) + ...``` +`g(x1) + g(x2) + ...` In probabilistic modeling this comes up when there are $N$ conditionally independent terms in a likelihood. Because of the @@ -71,10 +71,10 @@ call graph resulting in stable numerical evaluation. This version requires setting a sensible tuning parameter for good performance. The tuning parameter is the so-called `grainsize`. For the -`reduce_sum` version the `grainsize` is merely a suggested partial sum -size while for the `reduce_sum_static` version the `grainsize` +`reduce_sum` version, the `grainsize` is merely a suggested partial sum +size, while for the `reduce_sum_static` version the `grainsize` specifies the maximal partial sum size. While for `reduce_sum` a -`grainsize` of one commonly leads to good performance already (since +`grainsize` of 1 commonly leads to good performance already (since automatic aggregation is performed), the `reduce_sum_static` variant requires setting a sensible `grainsize` for good performance as explained in [more detail below](#reduce-sum-grainsize). @@ -156,9 +156,7 @@ for(i in 1:size(x)) { Logistic regression is a useful example to clarify both the syntax and semantics of reduce summation and how it can be used to speed up a typical -model. - -A basic logistic regression can be coded in Stan as: +model. A basic logistic regression can be coded in Stan as: ``` data { @@ -177,7 +175,6 @@ model { In this model predictions are made about the `N` outputs `y` using the covariate `x`. The intercept and slope of the linear equation are to be estimated. - The key point to getting this calculation to use reduce summation, is recognizing that the statement: @@ -194,13 +191,10 @@ 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 summation is useful. - -To use the reduce summation, a function must be written that can be used to compute -arbitrary partial sums of the total sum. - -Using the interface defined in [Reduce-Sum](#reduce-sum), such a function -can be written like: +reduce summation is useful. To use the reduce summation, a function +must be written that can be used to compute arbitrary partial sums of +the total sum. Using the interface defined in +[Reduce-Sum](#reduce-sum), such a function can be written like: ``` functions { @@ -213,7 +207,7 @@ functions { } ``` -And the likelihood statement in the model can now be written: +The likelihood statement in the model can now be written: ``` target += partial_sum(1, N, y, x, beta); // Sum terms 1 to N of the likelihood @@ -221,12 +215,12 @@ target += partial_sum(1, N, y, x, beta); // Sum terms 1 to N of the likelihood In this example, `y` was chosen to be sliced over because there is one term in the summation per value of `y`. Technically `x` would have -worked as well. Use whatever conceptually makes the most sense. - -Because `x` is a shared argument, it is subset accordingly with `start:end`. - -With this function, reduce summation can be used to automatically parallelize the -likelihood: +worked as well. Use whatever conceptually makes the most +sense for a given model, e.g. slice over independent terms like +conditionally independent observations or groups of observations as in +hierarchical models. Because `x` is a shared argument, it is subset +accordingly with `start:end`. With this function, reduce summation can +be used to automatically parallelize the likelihood: ``` int grainsize = 100; @@ -237,7 +231,7 @@ target += reduce_sum(partial_sum, y, The reduce summation facility automatically breaks the sum into roughly `grainsize` sized pieces and computes them in parallel. `grainsize = 1` specifies that the grainsize should -be estimated automatically. The final model looks like: +be estimated automatically. The final model looks as: ``` functions { @@ -269,7 +263,7 @@ model { The `grainsize` is a recommendation on how large each piece of parallel work is (how many terms it contains). When using the -non-static version, it is recommended to choose one as a starting +non-static version, it is recommended to choose 1 as a starting point as automatic aggregation of partial sums are performed. However, for the static version the `grainsize` defines the maximal size of the partial sums, e.g. the static variant will split the input sequence From 31d12d75b0913bc24a25bbc4385e9474187be4c8 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Tue, 7 Apr 2020 16:19:13 +0200 Subject: [PATCH 10/13] more review comments --- src/stan-users-guide/parallelization.Rmd | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index f2372f83d..f9b03a00e 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -32,9 +32,8 @@ The higher-order reduce with summation facility maps evaluation of a function `g: U -> real`, which returns a scalar, to a list of type `U[]`, `{ x1, x2, ... }`, and performs as reduction operation a sum over the results. For instance, for a sequence of ```x``` values of -type ```U```, ```{ x1, x2, ... }```, we might compute the sum: - -`g(x1) + g(x2) + ...` +type ```U```, ```{ x1, x2, ... }```, we might compute the sum as +`g(x1) + g(x2) + ...`. In probabilistic modeling this comes up when there are $N$ conditionally independent terms in a likelihood. Because of the From 021d3640fc0caf94780ee33f877b09a42e55c310 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 7 Apr 2020 10:27:13 -0400 Subject: [PATCH 11/13] Reworeded user's guide introduction to reduce_sum and lowered expectations about what `grainsize = 1` can do. (design-doc #17) --- .../higher-order_functions.Rmd | 20 ++--- src/stan-users-guide/parallelization.Rmd | 86 ++++++++++--------- 2 files changed, 53 insertions(+), 53 deletions(-) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index ea1c1d612..ef4c18c00 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -410,12 +410,12 @@ exactly. This implies that the order of summation determines the exact numerical result. For this reason, the higher-order reduce function is available in two variants: -* `reduce_sum`: Compute partial sums automatically. This usually - results in good performance without further tuning. -* `reduce_sum_static`: For the same input, always create the same call -graph. This results in stable numerical evaluation. This version -requires setting a tuning parameter which controls the maximal size of -partial sums formed. +* `reduce_sum`: Automatically choose partial sums partitioning based on a dynamic + scheduling algorithm. +* `reduce_sum_static`: Compute the same sum as `reduce_sum`, but partition + the input in the same way for given data set (in `reduce_sum` this partitioning + might change depending on computer load). This should result in stable + numerical evaluations. ### Specifying the Reduce-sum Function @@ -437,16 +437,14 @@ partial sums. `s1, s2, ...` are shared between all terms in the sum. * *`f`*: function literal referring to a function specifying the partial sum operation. Refer to the [partial sum function](#functions-partial-sum). * *`x`*: array of `T`, one for each term of the reduction, `T` can be any type, -* *`grainsize`*: recommended number of terms in each reduce call, set -to 1 to estimate automatically for `reduce_sum` while for -`reduce_sum_static` this determines the maximal size of the partial sums, type `int`, +* *`grainsize`*: For `reduce_sum`, `grainsize` is the recommended size of the partial sum. For `reduce_sum_static`, `grainsize` determinse the maximum size of the partial sums, type `int`, * *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type * *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, * *`...`*: remainder of shared arguments, each of which can be any type. -### The Partial-sum Function {#functions-partial-sum} +### The Partial sum Function {#functions-partial-sum} -The partial sum function must have the following signature where the types `T`, and the +The partial sum function must have the following signature where the type `T`, and the types of all the shared arguments (`T1`, `T2`, ...) match those of the original `reduce_sum` (`reduce_sum_static`) call. diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index f2372f83d..3585cefce 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -28,25 +28,21 @@ law](https://en.wikipedia.org/wiki/Amdahl's_law). ## Reduce-Sum { #reduce-sum } -The higher-order reduce with summation facility maps evaluation of a -function `g: U -> real`, which returns a scalar, to a list of type -`U[]`, `{ x1, x2, ... }`, and performs as reduction operation a sum -over the results. For instance, for a sequence of ```x``` values of -type ```U```, ```{ x1, x2, ... }```, we might compute the sum: +It is often necessary in probabilistic modeling to compute the sum of +a number of independent function evaluations. This occurs, for instance, when +evaluating a number of conditionally independent terms in a log-likelihood. +If `g: U -> real` is the function and `{ x1, x2, ... }` is an array of +inputs, then that sum looks like: `g(x1) + g(x2) + ...` -In probabilistic modeling this comes up when there are $N$ -conditionally independent terms in a likelihood. Because of the -conditional independence, these terms can be computed in parallel. If -dependencies exist between the terms, then this isn't possible. For -instance, in evaluating the log density of a Gaussian process then -summation of independent terms isn't applicable. +The `reduce_sum` function is a tool for automatically parallelizing these +calculations. For efficiency reasons the reduce function doesn’t work with the -element-wise evaluated function `g`, but instead requires the partial -sum function ```f: U[] -> real```, where ```f``` computes the partial -sum corresponding to a slice of the sequence ```x``` passed in. Due to the +element-wise evaluated function `g`, but instead the partial +sum function `f: U[] -> real`, where `f` computes the partial +sum corresponding to a slice of the sequence `x` passed in. Due to the the associativity of the sum reduction it holds that: ``` @@ -64,20 +60,20 @@ control of the user. However, since the exact numerical result will depend on the order of summation, Stan provides two versions of the reduce summation facility: -* `reduce_sum`: Automatically forms partial sums resulting usually in good - performance without further tuning. -* `reduce_sum_static`: Creates for the same input always the same -call graph resulting in stable numerical evaluation. This version -requires setting a sensible tuning parameter for good performance. - -The tuning parameter is the so-called `grainsize`. For the -`reduce_sum` version, the `grainsize` is merely a suggested partial sum -size, while for the `reduce_sum_static` version the `grainsize` -specifies the maximal partial sum size. While for `reduce_sum` a -`grainsize` of 1 commonly leads to good performance already (since -automatic aggregation is performed), the `reduce_sum_static` variant -requires setting a sensible `grainsize` for good performance as -explained in [more detail below](#reduce-sum-grainsize). +* `reduce_sum`: Automatically choose partial sums partitioning based on a dynamic + scheduling algorithm. +* `reduce_sum_static`: Compute the same sum as `reduce_sum`, but partition + the input in the same way for given data set (in `reduce_sum` this partitioning + might change depending on computer load). + +`grainsize` is the one tuning parameter. For `reduce_sum`, `grainsize` is +a suggested partial sum size. A `grainsize` of 1 leaves the partitioning +entirely up to the scheduler. + +For `reduce_sum_static`, `grainsize` specifies the maximal partial sum size. +With `reduce_sum_static` it is more important to choose `grainsize` +carefully since it entirely determines the partitioning of work. +See details in [more detail below](#reduce-sum-grainsize). For efficiency and convenience additional shared arguments can be passed to every term in the sum. So for the @@ -230,7 +226,7 @@ target += reduce_sum(partial_sum, y, ``` The reduce summation facility automatically breaks the sum into roughly `grainsize` sized pieces -and computes them in parallel. `grainsize = 1` specifies that the grainsize should +and computes them in parallel. `grainsize = 1` specifies that the `grainsize` should be estimated automatically. The final model looks as: ``` @@ -261,35 +257,41 @@ model { ### Picking the Grainsize {#reduce-sum-grainsize} -The `grainsize` is a recommendation on how large each piece of +For `grainsize` is a recommendation on how large each piece of parallel work is (how many terms it contains). When using the non-static version, it is recommended to choose 1 as a starting point as automatic aggregation of partial sums are performed. However, for the static version the `grainsize` defines the maximal size of the -partial sums, e.g. the static variant will split the input sequence -until all partial sums are just smaller than `grainsize`. Therefore, -for the static version it is more important to select a sensible -value. The rational for choosing a sensible `grainsize` is based on +partial sums, e.g. + +The rational for choosing a sensible `grainsize` is based on balancing the overhead implied by creating many small tasks versus creating fewer large tasks which limits the potential parallelism. -From empirical experience, the automatic grainsize determination works -well and no further tuning is required in most cases. In order to -figure out an optimal grainsize, think about how many terms are in the -summation and on how many cores the model should run. If there are `N` +In `reduce_sum`, `grainsize` is a recommendation on how to partition +the work in the partial sum into smaller pieces. A `grainsize` of 1 +leaves this entirely up to the internal scheduler. Ideally this will be +efficient, but there are no guarantees. + +In `reduce_sum_static`, `grainsize` is an upper limit on the worksize. +Work will be split until all partial sums are just smaller than `grainsize` +(and the split will happen the same way every time for the same data). +For the static version it is more important to select a sensible `grainsize`. + +In order to figure out an optimal `grainsize`, if there are `N` terms and `M` cores, run a quick test model with `grainsize` set -roughly to `N / M`. Record the time, cut the grainsize in half, and +roughly to `N / M`. Record the time, cut the `grainsize` in half, and run the test again. Repeat this iteratively until the model runtime -begins to increase. This is a suitable grainsize for the model, +begins to increase. This is a suitable `grainsize` for the model, because this ensures the caculations can be carried out with the most parallelism without losing too much efficiency. For instance, in a model with `N=10000` and `M = 4`, start with `grainsize = 25000`, and sequentially try `grainsize = 12500`, `grainsize = 6250`, etc. -It is important to repeat this process until performance gets worse! +It is important to repeat this process until performance gets worse. It is possible after many halvings nothing happens, but there might -still be a smaller grainsize that performs better. Even if a sum has +still be a smaller `grainsize` that performs better. Even if a sum has many tens of thousands of terms, depending on the internal calculations, a `grainsize` of thirty or forty or smaller might be the best, and it is difficult to predict this behavior. Without doing From 3edafecf088defd7af8715cd164412dc002a8a01 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 7 Apr 2020 11:55:37 -0400 Subject: [PATCH 12/13] Responded to reviews (design-doc #17) --- .../higher-order_functions.Rmd | 6 ++-- src/stan-users-guide/parallelization.Rmd | 31 ++++++++----------- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/src/functions-reference/higher-order_functions.Rmd b/src/functions-reference/higher-order_functions.Rmd index ef4c18c00..04da9692e 100644 --- a/src/functions-reference/higher-order_functions.Rmd +++ b/src/functions-reference/higher-order_functions.Rmd @@ -135,7 +135,7 @@ package MINPACK-1 [@minpack:1980]. The Jacobian of the solution with respect to auxiliary parameters is computed using the implicit function theorem. Intermediate Jacobians -(of the the algebraic function's output with respect to the unknowns y +(of the algebraic function's output with respect to the unknowns y and with respect to the auxiliary parameters theta) are computed using Stan's automatic differentiation. @@ -422,7 +422,7 @@ available in two variants: The higher-order reduce function takes a partial sum function `f`, an array argument `x` (with one array element for each term in the sum), a recommended `grainsize`, and a set of shared arguments. This representation allows -to parallelize the resultant sum. +parallelization of the resultant sum. \index{{\tt \bfseries reduce\_sum }!{\tt (F f, T[] x, int grainsize, T1 s1, T2 s2, ...): real}|hyperpage} @@ -437,7 +437,7 @@ partial sums. `s1, s2, ...` are shared between all terms in the sum. * *`f`*: function literal referring to a function specifying the partial sum operation. Refer to the [partial sum function](#functions-partial-sum). * *`x`*: array of `T`, one for each term of the reduction, `T` can be any type, -* *`grainsize`*: For `reduce_sum`, `grainsize` is the recommended size of the partial sum. For `reduce_sum_static`, `grainsize` determinse the maximum size of the partial sums, type `int`, +* *`grainsize`*: For `reduce_sum`, `grainsize` is the recommended size of the partial sum (`grainsize = 1` means pick totally automatically). For `reduce_sum_static`, `grainsize` determines the maximum size of the partial sums, type `int`, * *`s1`*: first (optional) shared argument, type `T1`, where `T1` can be any type * *`s2`*: second (optional) shared argument, type `T2`, where `T2` can be any type, * *`...`*: remainder of shared arguments, each of which can be any type. diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index 3585cefce..4ad110b71 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -36,14 +36,14 @@ inputs, then that sum looks like: `g(x1) + g(x2) + ...` -The `reduce_sum` function is a tool for automatically parallelizing these +`reduce_sum` and `reduce_sum_static` are tools for parallelizing these calculations. For efficiency reasons the reduce function doesn’t work with the element-wise evaluated function `g`, but instead the partial sum function `f: U[] -> real`, where `f` computes the partial sum corresponding to a slice of the sequence `x` passed in. Due to the -the associativity of the sum reduction it holds that: +associativity of the sum reduction it holds that: ``` g(x1) + g(x2) + g(x3) = f({ x1, x2, x3 }) @@ -68,12 +68,13 @@ reduce summation facility: `grainsize` is the one tuning parameter. For `reduce_sum`, `grainsize` is a suggested partial sum size. A `grainsize` of 1 leaves the partitioning -entirely up to the scheduler. +entirely up to the scheduler. This should be the default way of using +`reduce_sum` unless time is spent carefully picking `grainsize`. For picking a `grainsize`, see details [below](#reduce-sum-grainsize). For `reduce_sum_static`, `grainsize` specifies the maximal partial sum size. With `reduce_sum_static` it is more important to choose `grainsize` carefully since it entirely determines the partitioning of work. -See details in [more detail below](#reduce-sum-grainsize). +See details [below](#reduce-sum-grainsize). For efficiency and convenience additional shared arguments can be passed to every term in the sum. So for the @@ -219,15 +220,15 @@ accordingly with `start:end`. With this function, reduce summation can be used to automatically parallelize the likelihood: ``` -int grainsize = 100; +int grainsize = 1; target += reduce_sum(partial_sum, y, grainsize, x, beta); ``` -The reduce summation facility automatically breaks the sum into roughly `grainsize` sized pieces -and computes them in parallel. `grainsize = 1` specifies that the `grainsize` should -be estimated automatically. The final model looks as: +The reduce summation facility automatically breaks the sum into pieces +and computes them in parallel. `grainsize = 1` specifies that the +`grainsize` should be estimated automatically. The final model looks as: ``` functions { @@ -247,7 +248,7 @@ parameters { vector[2] beta; } model { - int grainsize = 100; + int grainsize = 1; beta ~ std_normal(); target += reduce_sum(partial_sum, y, grainsize, @@ -257,25 +258,19 @@ model { ### Picking the Grainsize {#reduce-sum-grainsize} -For `grainsize` is a recommendation on how large each piece of -parallel work is (how many terms it contains). When using the -non-static version, it is recommended to choose 1 as a starting -point as automatic aggregation of partial sums are performed. However, -for the static version the `grainsize` defines the maximal size of the -partial sums, e.g. - The rational for choosing a sensible `grainsize` is based on balancing the overhead implied by creating many small tasks versus creating fewer large tasks which limits the potential parallelism. In `reduce_sum`, `grainsize` is a recommendation on how to partition the work in the partial sum into smaller pieces. A `grainsize` of 1 -leaves this entirely up to the internal scheduler. Ideally this will be +leaves this entirely up to the internal scheduler and should be chosen +if no benchmarking of other grainsizes is done. Ideally this will be efficient, but there are no guarantees. In `reduce_sum_static`, `grainsize` is an upper limit on the worksize. Work will be split until all partial sums are just smaller than `grainsize` -(and the split will happen the same way every time for the same data). +(and the split will happen the same way every time for the same inputs). For the static version it is more important to select a sensible `grainsize`. In order to figure out an optimal `grainsize`, if there are `N` From bd93bf115f466c72bb4fa05634872f865c7a0f95 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 7 Apr 2020 12:34:52 -0400 Subject: [PATCH 13/13] Parallelization edits (design-doc #17) --- src/stan-users-guide/parallelization.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stan-users-guide/parallelization.Rmd b/src/stan-users-guide/parallelization.Rmd index 4ad110b71..a10fa609a 100644 --- a/src/stan-users-guide/parallelization.Rmd +++ b/src/stan-users-guide/parallelization.Rmd @@ -228,7 +228,7 @@ target += reduce_sum(partial_sum, y, The reduce summation facility automatically breaks the sum into pieces and computes them in parallel. `grainsize = 1` specifies that the -`grainsize` should be estimated automatically. The final model looks as: +`grainsize` should be estimated automatically. The final model is: ``` functions {