9.4 Reduce-Sum Function
Stan provides a higher-order reduce function for summation. A function 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) + ...
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)
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 choose partial sums partitioning based on a dynamic scheduling algorithm.reduce_sum_static
: Compute the same sum asreduce_sum
, but partition the input in the same way for given data set (inreduce_sum
this partitioning might change depending on computer load). This should result in stable numerical evaluations.
9.4.1 Specifying the Reduce-sum Function
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 parallelization of the resultant sum.
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, ...)
Returns the equivalent of f(x, 1, size(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.x
: array ofT
, one for each term of the reduction,T
can be any type,grainsize
: Forreduce_sum
,grainsize
is the recommended size of the partial sum (grainsize = 1
means pick totally automatically). Forreduce_sum_static
,grainsize
determines the maximum size of the partial sums, typeint
,s1
: first (optional) shared argument, typeT1
, whereT1
can be any types2
: second (optional) shared argument, typeT2
, whereT2
can be any type,...
: remainder of shared arguments, each of which can be any type.
9.4.2 The Partial sum Function
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.
(T[] x_subset, int start, int end, T1 s1, T2 s2, ...):real
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:
x_subset
, the subset ofx
a given partial sum is responsible for computing, typeT[]
, whereT
matches the type ofx
inreduce_sum
(reduce_sum_static
)start
, the index of the first term of the partial sum, typeint
end
, the index of the last term of the partial sum (inclusive), typeint
s1
, first shared argument, typeT1
, matching type ofs1
inreduce_sum
(reduce_sum_static
)s2
, second shared argument, typeT2
, matching type ofs2
inreduce_sum
(reduce_sum_static
)...
, remainder of shared arguments, with types matching those inreduce_sum
(reduce_sum_static
)