Folding in Parallel

Playing with monoids and discovering for oneself how frequently they turn up in various guises and how much is expressed by them is a delight. Not only it brightens the day of important but frankly boring data analytics. Not only it challenges to discover efficient (with the constant time and space operation) monoids, which demands ingenuity. In the words of Guy Steele, inventing efficient associative combining operators is a very big deal. An efficient monoid opens up parallel implementations -- not just parallel but embarrassingly parallel: The input sequence arbitrarily partitioned among workers, all running in parallel with no races, dependencies or even memory bank conflicts. This is the ideal for multi-core, GPU or distributed processing.

We show why monoids are a big deal, and also a big fun -- even more than expected. Horner rule and Boyer-Moore majority voting -- commonly regarded as prototypically inherently sequential -- turn out monoids, and hence embarrassingly parallelizable. As we generalize monoids to observable monoids, we discover a new, vertical way to combine them, enabling efficient iterated grouping and aggregation directly on raw deserialized (big) data, on multicore or multi-processor. Multicore OCaml and OpenMP examples bear it out.

When asked to sum an array, what first comes to mind?

For some readers, it may be

    int sum = 0;
    for(int i=0; i<N; i++) sum += arr[i];

For others, something like

    fold_left (+) 0 arr

(which, despite very different notation, ends up, with a good compiler, in the same machine code.)

The ideal, advocated by Guy Steele in his ICFP 2009 keynote, is

    Σaᵢ

This notation tells what without micromanaging how. Furthermore, there are many hows. For example,

    Σⁿ⁻¹₀ aᵢ = (((0 + a₀) + a₁) + ... + aₙ₋₁)

which is the left accumulation above. One can also accumulate from the right: fold_right in OCaml and similarly in other languages. One can also do
        Σⁿ⁻¹₀ aᵢ = (Σm-10 aᵢ) + (Σn-1m aᵢ)       where   0 < m < n
In other words, one may split the array in two (or more, and not necessarily equal) parts, sum them separately, and then tally up. The partial sums may be computed by the left- or right- accumulation, or by splitting again. Furthermore, they may be computed in parallel: not just in parallel, but embarrassingly in parallel. If the partial sum computations are assigned to separate cores, they can do their work without any synchronization or locking, without snooping at their neighbors, without even read memory conflicts, taking the full advantage of the memory bandwidth. Embarrassing parallelism especially suits GPU.

How would you compute the big sum depends on your circumstances (how big is the array, and how it is layed out in memory or disk or across nodes) and how many cores/nodes are available. All these decisions can be made later, to tuning -- or even to the run-time and or a compiler.

Summation is not unique in permitting such a variety of implementations, including the embarrassingly parallel ones. The same applies to any accumulation with an associative operation, and with a zero-like value to initialize the accumulator -- or, as we would say, to any monoid reduction.

Reduction is quite old (contemporaneous with the left and right folds), appearing already in APL in early 1960s. In APL, the summation of an array X is written, with characteristic brevity, as +/X: one should see it as `wedging in' the plus operation between the consecutive elements of the array. The array summation, contrasting the accumulation in the imperative language of the day with reduction has been a common example in APL tutorials since at least 1970, as far as I can see. (Looking a bit ahead, we will write the summation, in OCaml, as reduce (monoid Sum) -- not as brief as in APL, but more pronounceable, and also extensible, to `richer' monoids.) Modern array languages like Futhark or SaC indeed let the programmer write reductions in a similar fashion, compiling it to the efficient GPU, multicore etc. code, choosing the implementation strategy depending on the problem and the hardware. Monoid reduction got quite a boost in mid-2000s, thanks to Google's MapReduce. (We shall soon see why reduce is often accompanied by map.)

Implementing/compiling monoid reductions is well-researched and supported, especially in modern array languages or MapReduce frameworks. The problem is designing/finding the efficient monoid -- which, as Guy Steele said, is non-trivial and requires creativity. That's where the real fun starts.

References

Guy Steele: Organizing Functional Code for Parallel execution or, foldl and foldr considered slightly harmful
August 2009 (ICFP 2009, Keynote)

Olivier Danvy: ``Folding left and right matters: direct style, accumulators, and continuations''
Journal of Functional Programming, vol 33, e2. Functional Pearl, February 2023
Appendix A: A brief history of folding left and right over lists

According to Danvy, the first instances of fold_left and fold_right were investigated by Christopher Strachey in 1961(!). Reduce is part of APL (Iverson, 1962), where it is called `/': hence +/x sums the array x. Goedel recursor R in System T is a more general version of fold (called now para-fold). Church numerals are folds.

Monoid reduce is sometimes called `big operator' (by analogy with the capital Σ and Π for summation and multiplication, resp). Another name for them is `Eindhoven quantifiers'.

One may think that foldMap from the Haskell Foldable class is also related: it also talks about monoid reductions. However, foldMap is specified as an explicitly right-associative operator, and hence lacks the main property of map-reduce emphasized by Guy Steele: decoupling of specification from implementation. One may argue the Foldable type class itself is to blame: according to it, a collection has a single foldMap, whereas in fact there may be many reduction implementations for the same collection.

Before getting to monoid reductions, lets review folds -- for ease of reference and because they will appear frequently, contrasted with monoids. Left fold over a sequence is particularly important as a notation for a state machine, a general pattern of a sequential, or streaming computation.

Folding, or sequential accumulation over a sequence (a list, for concreteness) is defined as

    fold_left  : ('z -> 'a -> 'z) -> 'z -> 'a list -> 'z
    fold_right : ('a -> 'z -> 'z) -> 'a list -> 'z -> 'z

Actually, there are two operations: the left and the right fold. Their meaning should be clear from the following example:

    fold_left (+) 0 [1;2;3;4]  ≡ (((0 + 1) + 2) + 3) + 4
    
    fold_right (+) [1;2;3;4] 0 ≡ 1 + (2 + (3 + (4 + 0)))

In this case, of the folding function being addition, the results are identical: both expressions sum the list. Generally, left and right folds produce different results: try, for example, subtraction as the folding function. The types of the list elements 'a and of the accumulator 'z need not be the same: for example,

    fold_left (fun z _ -> z + 1) 0 l

computes the length of any list,

    fold_left (fun z x -> x :: z) [] l

reverses the list l and

    fold_right (fun x z -> if p x then x::z else z) l []

filters the list: omits the elements for which the predicate p returns false. Many other operations on lists (in fact, all of them) can be expressed as folds. Fold is indeed the general pattern of sequential stateful processing of a sequence.

The are alternative definitions of fold: sometimes the last two arguments of the right fold are swapped. If the arguments of the right folding function are likewise swapped, then the left and the right fold have the same signature. Their behavior, the association pattern, is still different. Olivier Danvy, see below, traces the history of list folds and the argument order.

For concreteness we showed folds over lists, but similar operations exist over arrays, streams, files, trees, dictionaries and any other collections.

In this article, by reduce we always mean the reduce over a monoid. A monoid is a set (called `carrier set') with an associative binary operation which has a neutral (also called unit, or zero) element. Concretely, in OCaml

    module type monoid = sig
      type t                                (* carrier set *)
      val zero : t
      val op : t -> t -> t
    end
    type 'a monoid = (module monoid with type t = 'a)

The operation op must be associative (and total: defined for all t values) and

    op zero x = op x zero = x

must hold for every element x of the monoid. In Google MapReduce, the operation op is also taken to be commutative. We do not impose such requirement. We will, however, emphasize efficiency: the monoids we aim for should have a constant time and space operation op. The obvious example is

    module Sum = struct
      type t = int
      let zero = 0
      let op = (+)
    end

Reduce over a sequence (here, a list) is the operation

    reduce : 'a monoid -> 'a list -> 'a

with the behavior that can be illustrated as

    reduce monoid []  ≡ monoid.zero
    reduce monoid [x] ≡ x             (* for any x and monoid *)
    
    reduce (module Sum) [1;2;3;4] ≡ 1 + 2 + 3 + 4

One may say that reduce `wedges in' the monoid operation between the consecutive elements of the sequence. Since op is associative, the parentheses are not necessary. Therefore, unlike the left and the right fold, there is only one reduce.

Folding over a sequence is executing a finite state machine whose input is the sequence: viz., fold_left f z l runs the finite state machine with the initial state z and the transition function f on the input l, and returns the final state. It is an inherently sequential computation. As said already, the left fold is just a different notation for an accumulating for-loop. (The right fold is similar, but accessing the sequence right to left.)

Reduce, however, can be implemented in various ways: for example,

    let reduce (module M) l = List.fold_left M.op M.zero l
    let reduce (module M) l = List.fold_right M.op l M.zero 

or

    let reduce m l = split l |> map reduce' m |> reduce'' m

That is, split the sequence in two or more parts (using an appropriately defined split : 'a list -> 'a list list), reduce each part (using some implementation of the reduce, denoted here reduce') and then tally up the results, using perhaps a different implementation of reduce. (Using this code as is, with regular OCaml lists, unlikely yields any benefits; however, there are better sequences, e.g., chunk lists, arrays or seekable file handles; intermediate lists can be eliminated by fusion. We shall see examples later.) The attraction of this implementation (schema) is that part reductions can all be done in parallel -- embarrassingly in parallel. For example, reduce' above can be running in its own thread (multicore OCaml Domain or an OpenMP task) as we shall see later. To stress again, no locking or synchronization is needed. Reduce thus offers quite a lot of flexibility: whether to reduce sequentially or in parallel, how many parts to split the sequence in, what implementation to use for part reduction and the final tallying up. The decisions depend on the (estimated) size of the sequence and the available hardware (cores), and can be left to later tuning or the run-time.

Since the (left or right) fold commits us to sequential evaluation, Guy Steele called it `slightly harmful'. He urged using reduce, as far as possible, because it is flexible and decouples the algorithm from the execution strategy. The strategy (sequential, parallel, distributed) and data partitioning can be chosen later, depending on circumstances and available resources.

Thus the main question is: can we convert fold into reduce? The rest of the article answers it.

References

Guy Steele: Organizing Functional Code for Parallel execution or, foldl and foldr considered slightly harmful
August 2009 (ICFP 2009, Keynote)

monoid_reduce.ml [42K]
Complete code for the article

How to zip folds
A complete library of fold-represented lists, demonstrating that all list processing operations can be expressed as folds

Accumulating tree traversals, a better tree fold with applications to XML parsing

Recall, monoid is a set with two operations (satisfying particular properties) -- in other words, an Algebra, of the following signature.

    module type monoid = sig
      type t
      val zero : t
      val op : t -> t -> t
    end

Different monoids -- monoids algebras -- are implementations of the signature. The operation op must be associative with zero as its left and right unit. Furthermore, we place special emphasis on the efficient, that is, constant space and time, operation op.

The obvious examples are

    module Sum = struct
      type t = int
      let zero = 0
      let op = (+)
    end

and Prod -- and also Min (with min as the operation) and Max; All (carrier bool, zero is true, op is conjunction) and Some (carrier bool, zero is false, op is disjunction). Lists (or sequences, in general) are prototypical monoids, with empty list as zero and append (or, reverse append) as the operation. This monoid is not efficient, however: concatenating two non-empty sequences produces a bigger sequence; furthermore, for lists, concatenation is not constant time.

One must also note the EndoF monoid: its carrier is the set of t->t functions (set-theoretical maps in general) for some type t; zero is the identity and op is the composition (for related EndoF', op is the left-to-right composition). EndoF is of great theoretical importance, giving rise to `monoid actions'. Practically however (if we are talking about programming language functions), it is wasteful; we see an example later. Although function composition is usually constant time, it creates a new closure (incorporating, by reference) its argument functions, and, therefore, not constant-space.

Homomorphism

That monoids appeared in pairs in our list, e.g., All and Some, hinted at a connection among them. The most important connection is called homomorphism. Homomorphism h from monoid A to monoid B is a map from the carrier of A to the carrier of B (in programming language terms, A.t -> B.t function) that satisfies

    h A.zero = B.zero         h (A.op x y) = B.op (h x) (b y)

for any elements x and y of the monoid A. One says, h `commutes with operations'. For example, boolean negation is a homomorphism from Some to All or vice versa, justified by the de Morgan laws. Exponentiation is the homomorphism from Sum to Prod. Try to find interesting homomorphism, e.g., from Sum to Some. A homomorphism from a monoid A to EndoF has a special name: monoid action.

There is another important homomorphism. Consider integer lists, which are monoid, as we saw already. A homomorphism h from integer lists to Sum should satisfy

    h [] = 0     h (x @ y) = h x + h y

List summation satisfies exactly these conditions: reduce is a homomorphism. List.length satisfies too: the homomorphism conditions, hence, do not define reduction uniquely. We need to extend the notion of monoids, so to precisely specify reduce as a sequence homomorphism.

Map-reduce and Collection Monoids

There are many hints that the monoid interface is too minimalistic. First, compare the types of the list fold

    fold_left  : ('z -> 'a -> 'z) -> 'z -> 'a list -> 'z

and reduce

    reduce : 'a monoid -> 'a list -> 'a

which, after inlining 'a monoid, becomes

    reduce : ('a -> 'a -> 'a) -> 'a -> 'a list -> 'a

The fold_left type is more general. This generality lets fold_left express, for example, the length of a list of an arbitrary type:

    let length 'a list -> int = List.fold_left (fun z _ -> z + 1) 0

Although list length can also be computed as a reduction, it is not a mere reduction:

    let length : 'a list -> int = map (Fun.const 1) >> reduce (module Sum)

(where >> is the left-to-right function composition). That is, we need map, to first convert element types (and values) to the monoid to accumulate. That's why reduce is often preceded by map, to do this adjustment.

If the type t in the monoid signature is taken to be abstract, the only monoid elements we can get hold of are zero, op zero zero, etc. -- all zeros. We need a way to obtain other elements: we need generators.

In abstract Algebra, generators are often introduced as zero-arity operations, or constants -- perhaps infinitely many constants. In the programming language context, perhaps its best to think of generators as being indexed by some set, or, as a mapping (function) from the index set to monoid elements:

    module type coll_monoid = sig
      include monoid
      type ix
      val gen : ix -> t
    end
    type ('a,'e) coll_monoid = (module coll_monoid with type t='a and type ix='e)

The index set may be empty, unitary (e.g., unit), finite (e.g., bool) or infinite. Whereas the type t of monoid elements may be abstract, the index type ix should be concrete. Fegaras and Maier call (a simpler version of) this interface a collection monoid.

The earlier Sum etc. monoids are trivially extended to collection monoids:

    module Sum = struct
      type t = int
      let zero = 0
      let op = (+)
      type ix = int
      let gen = Fun.id
    end

In other words, if the monoid type is concrete, we can always find/construct its values `out-of-band', so to speak, without using the monoid interface. A non-trivial example of a collection monoid is

    let list_monoid : type a e. (a list, a) coll_monoid = 
        (module struct
          type ix = a
          type t = a list
          let zero = []
          let op = (@)
          let gen x = [x]
        end)

For this monoid, all of its elements can be obtained by via monoid operations alone: as zero, gen i for some i, or applying op. In other words, it is a generated monoid.

A homomorphism over collection monoids must also commute with the operation gen: in other words, a homomorphism from collection monoid A to collection monoid B is the function h : A.t -> B.t satisfying

    h A.zero = B.zero     h (A.op x y) = B.op (h x) (b y)
    h (A.gen i) = B.gen i

In particular, the homomorphism h from the int list collection monoid to Sum must satisfy

    h [] = 0     h (x @ y) = (h x) + (b y)      h [i] = i

which is exactly the list summation. The conditions now specify the homomorphism uniquely. List length is also a homomorphism, but to a different, Sum-like monoid

    struct include Sum let gen = Fun.const 1 end

With collection monoids, list length is expressed directly as reduction (homomorphism), without the need for map. One may say, that with collection monoids, a homomorphism from a sequence is map-reduce.

Thus, a collection monoid uniquely specifies the homomorphism from a sequence. In other words, one may define map_reduce : ('a,'e) coll_monoid -> 'e list -> 'a. Before we do this, however, let's consider an alternative operation for generators:

    val opgen : t -> ix -> t  

It is equivalent, or inter-convertible, with gen:

    opgen t x = op t (gen x)
    gen x = opgen zero x

However, often opgen can be implemented more efficiently. There is another reason.

If we are to implement the collection monoid reduction as left-accumulation, it should look like

    let map_reduce : type t e. (t,e) coll_monoid -> e list -> t = 
      fun (module M) -> List.fold_left M.opgen M.zero 

In other words, opgen is exactly the folding function.

Sequential and parallel map-reduce over arrays

As stressed before, monoid reduction permits many other implementations (all giving the same result). Here we show just two, using an array slice as a collection. (OCaml does not provide array slices in the standard library; therefore, we have to define it ourselves.)

    type 'a array_slice = {arr:'a array; from:int; upe:int} 
    let whole_slice (arr:'a array) : 'a array_slice = 
      {arr; from=0; upe=Array.length arr}

The sequential, left-associative implementation:

    let map_reduce_arr : type t e. (t,e) coll_monoid -> e array_slice -> t = 
      fun (module M) {arr;from;upe} -> 
      let acc = ref M.zero in
      for i=from to upe-1 do
        acc := M.opgen !acc arr.(i)
      done;
      !acc

and the parallel implementation, over ncores cores.

    let map_reduce_par_n : 
     type t e. int -> (t,e) coll_monoid -> e array_slice -> t = 
      fun ncores (module M) {arr;from;upe} -> 
      let len = upe - from in
      let chunk = len / ncores in
      List.init ncores Fun.id |>
      List.map (fun i -> 
        Domain.spawn (fun () -> 
            map_reduce_arr (module M) 
             {arr; from = i*chunk; upe = min len ((i+1)*chunk)}))
        |>
        map_reduce (module struct include M 
            type ix = M.t Domain.t
            let gen = Domain.join
            let opgen x y = op x (Domain.join y) end)

It is a bit naive: We should be checking first if the length is big enough: if not, just do the sequential map_reduce_arr, since the ever-present overhead of parallel evaluation would otherwise dominate. We should also make chunk to be a multiple of the page size. Also, we should use Domain 0 to do some work, too. See the end of Nested grouping-aggregation: Vertical monoid composition for the real example.

Products

So far, we built monoids from scratch. The is another way: combining, or composing, existing monoids. The most obvious and straightforward composition is pairing, or a (collection) monoid product:

    module Product(M1:coll_monoid)(M2:coll_monoid with type ix = M1.ix) :
        (coll_monoid with type t = M1.t * M2.t and type ix = M1.ix)

It is literally just a pair of monoids, operated on elementwise. The two monoids should have the same index set, to generate elements of both monoids. An example is computing sequence average:

    module Av = 
      Product(Sum)(struct include Sum let opgen x _ = x+1 end)

The map-reduce over Av yields both the sum and the length of the sequence in one pass. (How would you modify Av to avoid overflows in sum computation? After all, the average is well-defined for any non-empty sequence, no matter how long and how big its elements are.)

Observable Monoids and the Conjugate Transform

If the type t in the (collection) monoid interface is indeed abstract, the result of a monoid reduction is a value of the abstract type t, which we cannot use outside the monoid. We need one more operation: observing the monoid, as a value of some concrete type obs.

    module type obs_monoid = sig
      include coll_monoid
      type obs
      val obs : t -> obs
    end

The other motivation comes from computing the average of a sequence: monoid reduction with the Av monoid above gives a pair of numbers: the sum and the length. To get the desired average, we have to `observe' the pair: do the final division. The operation obs looks like the complement of gen, which is also gratifying.

Any collection monoid can be observed as is:

    module AsIs(M:coll_monoid) : 
           (obs_monoid with type ix = M.ix and type obs = M.t) = struct
      include M
      type obs = M.t
      let obs = Fun.id
    end

A non-trivial observable monoid is

    module Average = struct
      include Av
      type obs = float
      let obs (s,l) = float_of_int s /. float_of_int l
    end

Map-reduce with the observable monoid takes the form

    let map_reduce_conj (type i) (type o)
        (module M:obs_monoid with type ix = i and type obs = o) (l:i list) : o = 
      map_reduce (module M) l |> M.obs

performing map-reduction and then observing the result. The overall process -- injecting sequence elements into the monoid using gen, reducing, and observing the result -- is called conjugate transform. It is the general monoid reduction.

References

monoids.ml [6K]
OCaml code for various monoids and monoid compositions

Leonidas Fegaras, David Maier: Towards an Effective Calculus for Object Query Languages. SIGMOD 1995, pp. 47-58.

Algebra
Introduction to the mathematical concept of Algebras for computer scientists, specifically for tagless-final programmers

Returning to Guy Steele's exhortation of preferring reduce over fold, one may wonder if one can always do that. In that section we see that fold can indeed always be expressed in terms of reduce -- but in a way that is not useful or interesting. Indeed,

    fold_right f [1;2;3;4] z
    = f 1 (f 2 (f 3 (f 4 z)))
    = (f 1 · f 2 · f 3 · f 4) z
    = map_reduce (endof_monoid f) |> (fun g -> g z)

for any f, z and the sequence of appropriate types. The very similar expression can be written for the left fold (left as an exercise to the reader). The key idea is that function composition is associative. Here:

    let endof_monoid : type a e. (e -> a -> a) -> (a->a,e) coll_monoid =
     fun f -> 
        (module struct
          type t = a -> a
          let zero = Fun.id
          let op g h = fun x -> g (h x)     (* may also consider >> *)
          type ix = e
          let gen = f
          let opgen m i = fun x -> m (f i x)
        end)

Since homomorphism to endof_monoid is called monoid action, fold is hence called `iterative action' -- the construction that is quite known in mathematics.

Unfortunately, this trivial answer is of little practical use. As we see from the example, the map-reduction creates a closure (f 1 · f 2 · f 3 · f 4), which is finally applied to z. The composed closure (f 1 · f 2 · f 3 · f 4) has the same structure as the original list -- and hence the same size (actually, a few times larger: not only we have to store the list elements themselves, we have to store the function pointers f, plus the overhead). In effect we have built an intermediate data structure of unbounded size. Although composing the closures can be parallelized, the useful work is done only when the big closure composition is finally applied to the initial accumulator z -- at which point the folding function f is applied step-by-step, sequentially, just like in the original fold_right f l z. This trivial reduction of fold only wastes time and space.

The problem hence is not merely to express fold as reduce, but do so efficiently, without undue and unbounded overhead. Only then we can profitably fold in parallel.

References

Jeremy Gibbons: Origami Programming for Fun and Profit
<https://www.cs.ox.ac.uk/publications/publication16637-abstract.html>
Exposition of a fold as endof_monoid

Representing a sequential algorithm (fold) as monoid reduce, opening up embarrassing parallelization, only makes sense if we can do so efficiently, without intermediary data of unbound size and without taking undue time. In other words, if folding works in constant time per sequence element and in constant (working) space, so should reduce.

The general approach, the conjugate transform introduced above, is a principle, or a schema. It does not tell how to actually find the efficient monoid, or if it even exists. In fact, it is often non-trivial and requires ingenuity. To appreciate this, here a few exercises for the reader. To stress, in all cases we are interested in the efficient monoid.

Iterated subtraction

Consider subtractive folding, mentioned in passing earlier. Subtraction is not associative and therefore fold_left (-) is not per se an instance of reduce -- and neither is fold_right (-) (and they obviously give different results). Yet both are expressible as map-reduce over efficient monoids. Try to find them. Hint: one of the two cases is quite easier.

ArgMax

Finding the largest element of a sequence is a reduction, over the Max monoid. What if we want to find not just a largest element, but, at the same time, also its position in the sequence -- the left-most position, if there are several occurrences. What about the right-most position? In other words, compute max and argmax as an efficient monoid reduction, in one pass.

Parentheses matching

Finding if in a string of open and close round parentheses, all parentheses match can be done as an efficient map-reduction. Try to find the corresponding monoid. As the next question, assume the string may contain other, non-parentheses characters. Can you modify the monoid that the reduction would tell not only if the parentheses match, but also, in the case of a mismatch, the position of the left-most unopened close parenthesis and the position of the right-most unclosed open parenthesis. What about the right-most unopened close parenthesis and the left-most unclosed open parenthesis?

Horner rule, or schema, is a widely-used sequential algorithm to efficiently evaluate a polynomial. A similar accumulating pattern also frequently occurs in parsing. Although the algorithm is inherently sequential, on the face of it, it turns out representable as a monoid reduction. Our reduction is more general and flexible than the other known ways to parallelize the Horner rule.

As an illustrating example, we take the conversion of a sequence of digits (most-significant first) to the corresponding number -- a simple parsing operation. It is an accumulating sequential operation and can be written as fold:

    let digits_to_num = 
      List.fold_left (fun z x -> 10*z + (Char.code x - Char.code '0')) 0

For example, digits_to_num ['1'; '7'; '5'] gives 175. One can easily factor it as a composition of a map and

    List.fold_left (fun z x -> 10*z + x) 0

which is essentially the Horner rule: evaluating the polynomial Σ d_i bⁱ at b=10. The folding function fun z x -> 10*z + x is not associative, however.

Still, there is a way to build a monoid:

    module Horner : (obs_monoid with type ix = int and type obs = int) =
    struct 
      type t = int * int
      let zero = (0,1)
      let op (x,b1) (y,b2) = (x*b2+y, b1*b2)
      type ix = int
      let gen x = (x,10)
      type obs = int
      let obs = fst
      end

and digits_to_num is a reduction

    let digits_to_num = 
      map_reduce_conj (module struct
          include Horner
          type ix = char
          let opgen m x = opgen m (Char.code x - Char.code '0') end)

Thus, the Horner rule can be expressed as map-reduce and can hence be evaluated embarrassingly in parallel.

The Horner rule has the feel of the parallel prefix problem and the famous parallel scan of Guy Blelloch (PPoPP 2009). Horner rule does not require reporting of intermediate results, however. Our monoid reduction differs from Blelloch's even-odd interleaving. It also differs from monoid-cached trees from Guy Steele's ICFP09 keynote. We do not rely on any special data structures: We can work with plain arrays, partitioning them across available cores.

Since Horner method is pervasive in polynomial evaluation, there are approaches to parallelize it. The most common is even/odd splitting of polynomial coefficients. Although it may be well suitable for SIMD processing, the even/odd splitting is bad for multicore since it creates read contention (bank conflicts) and wastes cache space. Our monoid reduce lets us assign different memory banks to different cores, for their exclusive, conflict-free access.

On the other hand, Estrin's scheme may be seen as a very particular instance of our monoid construction: binomial grouping. Our monoid reduction allows for arbitrary grouping, hierarchically if desired, and of not necessarily of the same size.

References

monoid_reduce.ml [42K]
Complete code for the article

Boyer-Moore majority voting is the algorithm for finding the majority element in a sequence -- that is, the element that occurs more than half of the time. The algorithm requires one pass over the sequence and the constant amount of working space. It is important to keep in mind that the return value is the majority element if it exists. The algorithm in general cannot tell if the sequence has majority: if majority does not exist, the return value is an arbitrary element, not necessarily the most frequently occurring. In general, therefore, one needs the second pass over the sequence to count the occurrences of the returned element and verify it is indeed the majority. The second pass may be avoided or be unnecessary in some cases.

The Boyer-Moore algorithm is called a prototypical streaming algorithm and is inherently sequential. It may be surprising, therefore, that it may be presented as monoid reduction and hence performed just as efficiently in parallel -- in fact, embarrassingly in parallel. The monoid reduction implementation is actually simple. What is complex is convincing ourselves that it indeed works and does the right thing, in all cases.

The Boyer-Moore algorithm is a state machine processing, ingesting input element-by-element. It can hence be implemented as fold (assuming the input is given as a non-empty list):

    let bm : 'a list -> 'a = function h::t -> 
      let fsm (c,m) x =
          if c = 0 then (1,x)
          else if m = x then (c+1,m)
          else (c-1,m)
      in List.fold_left fsm (1,h) t >> snd

The state is the counter c and the candidate majority m. At the end of the algorithm, when the entire sequence is scanned, m is the majority element (if the sequence has one). If the counter c at the end is more than half of the sequence length, the sequence has majority and m is the majority element. Otherwise, we have to make another pass through the sequence to verify that m is indeed the majority. The verification pass may be unnecessary if we have a prior knowledge that majority exists -- or if we are satisfied with any element if it does not.

Although the algorithm seems inherently sequential, it can in fact be represented as monoid reduction. The monoid is simple and efficient:

    type 'a bm = int * 'a
    type t = a bm
    let zero = (0,z)
    let op (c1,m1) (c2,m2) =
     if c1 = 0 then (c2,m2) else
     if c2 = 0 then (c1,m1) else
     if m1 = m2 then (c1+c2,m1) else
     if c1 > c2 then (c1-c2, m1) else
     (c2-c1, m2)
    type ix = a
    let gen x = (1,x)
    let opgen m x = op m (gen x)

Monoid elements (the carrier set) are pairs: of a natural number c and of a sequence element m. When c is zero, m may be arbitrary; in the above code we use z for such an arbitrary element. We could have avoided specifying it by defining a more sophisticated data type for monoid carrier, or using 'a option. The reasoning below would get messier; therefore, we keep the above definition for the sake of clarity.

One can see that op is commutative and that zero is indeed the zero of op. However, is op associative? That is, is bm_monoid really a monoid? It is easy to check:

    let m1 = (2,10) and m2 = (3,9) and m3 = (2,8)
    op (op m1 m2) m3  ⇝ (1, 8) 
    op m1 (op m2 m3)  ⇝ (1, 10)

Unfortunately, bm_monoid is not a monoid.

To cut the suspense, it turns out that bm_monoid is `morally' a monoid, and reduction with it indeed gives the right result, in all cases. Proving it however is not that straightforward.

For the sake of proof, we extend bm_monoid to:

    type 'a ebm = int * 'a * ('a*'a) list
    type t = a ebm
    let zero = (0,z,[])
    let op (c1,m1,l1) (c2,m2,l2) =
     if c1 = 0 then (c2,m2,l1@l2) else
     if c2 = 0 then (c1,m1,l1@l2) else
     if m1 = m2 then (c1+c2,m1,l1@l2) else
     if c1 > c2 then (c1-c2, m1, repeat c2 (m1,m2) @ l1 @ l2) else
     (c2-c1, m2, repeat c1 (m1,m2) @ l1 @ l2)
    type ix = a
    let gen x = (1,x,[])
    let opgen m x = op m (gen x)

Here, the infix operation @ is list concatenation and repeat (n:int) (x:'a) : 'a list produces the list of length n with elements all equal to x.

The carrier set of bm_monoid_ext is 'a ebm: triples, of a natural number c, a sequence element m, and a list of pairs with distinct components (that is, fst of any pair is not equal to its snd). We note that op preserves the invariant that all pairs in l have distinct components.

The operation op of thus extended bm_monoid_ext is not associative. Not only bm_monoid_ext is not a monoid, it is also not efficient: the list concatenation operations in op are neither constant space nor constant time.

Let us introduce the relation ≈ on that set, defined as follows: x ≈ y just in case flatten x is equal to flatten y, where

    let flatten : 'a ebm -> 'a multiset = fun (c,m,l) ->
      repeat c m @ List.map fst l @ List.map snd l

That is, flatten x collects all elements mentioned within x in a multiset. Therefore, x ≈ y iff the multiset of elements of x is equal to the multiset of elements of y. One easily sees that the relation ≈ is reflexive, commutative and associative. That is, it is an equivalence relation. It also has useful for us properties, as follows.

Proposition: bm_monoid_ext is a monoid modulo ≈. That is,

    op (0,z,[]) (c,m,l) ≈ op (c,m,l) (0,z,[]) ≈ (c,m,l)
    op (c1,m1,l1) (op (c2,m2,l2) (c3,m3,l3)) ≈ op (op (c1,m1,l1) (c2,m2,l2)) (c3,m3,l3)

The associativity property follows from the fact op preserves elements occurring in (c,m,l), including their multiplicities (and the multiset equality is associative).

Moreover, ≈ is adequate for the task of finding the majority.

Proposition (Adequacy): if (c1,m1,l1) ≈ (c2,m2,l2) and the multiset flatten (c1,m1,l1) has majority, then c1>0, c2>0, and m1=m2. Informally, ≈ preserves the found majority (if the majority exists). The proposition is the consequence of the fact that majority is invariant to sequence permutation, and the following representation lemma.

Lemma (Representation): if a sequence (multiset) flatten (c,m,l) has majority, then c>0 and m is the majority element.

The key idea is that (c,m,l) is a representation of a partition of the multiset flatten m: specifically, a partition of the multiset into pairs of distinct elements and the remainder (which, if exists, is comprised of one or more copies of a particular element m). Such a partition is generally not unique: consider the sequence [1; 2; 3; 3; 3] and its two partitions (3,1,[(1,3);(2,3)]) and (3,3,[(1,2)]). It has however the property specified in the lemma.

For the proof we note that if the list of distinct pairs l has n elements, it has 2*n of values. Not all of them are distinct. However, there is no value that occurs more than n times; otherwise, by the pigeonhole principle, some pair in l would have had equal components. Therefore, no value occurring in (c,m,l) that is distinct from m (if c>0) can be the majority element. If the majority exists, it has to be m, and c>0.

Our proofs are inspired by Boyer and Moore's original proofs. They however used induction (which is appropriate for the sequential algorithm). We, in contrast, used equational reasoning (and the pigeonhole principle) -- but no induction.

As an aside, it may seem that induction is the only proof method used in computer science -- many widely used textbooks present nothing but induction. In my opinion, induction is overused. Mathematics has many other proof methods.

The final step is to observe that the part of bm_monoid_ext that we actually care about -- the element m and the count c -- is computed with no regard to the list l. The list is a `ghost' list, so to speak, needed only for proving correctness but not for computing the majority. Therefore, we may omit it and recover the bm_monoid.

In the upshot, the majority voting may hence be performed as map-reduce, over the monoid bm_monoid. It is just as efficient as the original Boyer-Moore algorithm: one pass over the sequence using constant and small working space. However, the input may be split arbitrarily among available workers, who would do the search in parallel, without any interference or dependencies. Majority voting can be done embarrassingly in parallel.

References

<https://en.wikipedia.org/wiki/Boyer%E2%80%93Moore_majority_vote_algorithm>

monoid_reduce.ml [42K]
Complete code for the article

Lambert Meertens: Reducing hopeful majority. The Squiggolist, 1(1):5, 1989
<https://www.kestrel.edu/people/meertens/publications/papers/Reducing_hopeful_majority.pdf>
The paper also noted the lack of associativity and informally argued (invoking indeterminism) that it is not required for consistency. We have shown how associativity can be formally regained without invoking indeterminism, using the `ghost' list of distinct pairs. This `ghost' list is the crucial ingredient to be able to define the equivalence relation `≈'.

Another instructive example of the monoid reduction is nested grouping and aggregation -- a typical problem in data analysis. That grouping and aggregation per se are expressible as a monoid comprehension -- map-reduce -- was insightfully articulated by Fegaras and Maier in their pioneering paper on monoid comprehensions. In our case, however, the result of a grouping-aggregation is grouped-aggregated again -- and again. An earlier article described a sequential optimal solution to the problem. Curiously, it also featured monoids. However, their composition was not explicated. Here we show that `piling up' particular monoids still gives a monoid -- and hence the embarrassingly parallel processing.

The running example, as in the earlier article, is an Advent of Code 2022 problem, summarized as:

The input is a sequence of numbers separated by pipes into chunks where each chunk contains multiple comma-separated numbers. The task is to compute the sum of each chunk and find the maximum across all chunks.

A sample input is as follows:
100,200,300|400|500,600|700,800,900|1000

The earlier article showed an efficient -- single-pass with no buffering -- sequential solution. We now present an embarrassingly parallel one. Emphatically, the input can be split arbitrarily, for example

    100,20  0,300|40    0|500,600   |700,800  ,900|1000

haphazardly cutting through the numbers and delimiters. We may mmap the input file and assign a range of pages to a worker (a CPU core), not worrying about how or if the data are aligned with page boundaries. The workers may then all run in parallel, in constant memory, with no dependencies or races -- with the linear speed-up. Such a processing is also very fitting for GPU.

A more realistic variation of the problem is, say, finding the 10 most frequent words in a UTF-8--encoded text file.

Warm-up

To build intuition, let's first consider a simple version of the problem. The input is a sequence of items of type α separated by delimiters. For example, in the case α is char, or digit, to be precise:

    let adv22_input = 
    ['1'; '0'; '0'; ','; '2'; '0'; '0'; ','; '3'; '0'; '0'; '|'; '4'; '0'; '0';...]

We take the input here to be a list for simplicity; in the full version, the input is a byte array or a file. The task is to group the items between the delimiters (commas and bars).

Formally, the task is to convert (α+δ)^* to (α^*↑)^* where α^*↑ is a maximal subsequence not containing the delimiter. Here, δ is the type (set) of delimiters. We use `+' for a union (of two types or sequences) and multiplication for concatenation; the superscripted star is Kleene star. Observe:

(α+δ)^*  =  (α^*δ^*)^*  =  α^* + α^*δ(α^*δ)^*α^*  =  α^* + α^*δ(α^*↑δ)^*α^*

using the almost standard algebra (only the multiplication, being sequence concatenation, is not commutative). In words: A sequence (α+δ)^* either:

• has no delimiters: α^*
• has one: α^*δα^*
• has more than one:α^*δ(α^*↑δ)⁺α^*

In the latter case, the α^* sub-sequences between the delimiters are maximal -- by definition. This representation is closed under concatenation (`multiplication'):

(α+δ)^*  =  (α+δ)^* (α+δ)^*
 =  (α^* + α^*δ(α^*↑δ)^*α^*) (α^* + α^*δ(α^*↑δ)^*α^*
 =  α^* + α^*δ(α^*↑δ)^*α^* + α^*δ(α^*↑δ)^*α^* α^*δ(α^*↑δ)^*α^*  =  α^* + α^*δ(α^*↑δ)^*α^*

The original (α+δ)^* is a monoid: the operation is concatenation, zero is the empty sequence. As we have just shown, the representation α^* + α^*δ(α^*↑δ)^*α^* is also a monoid; the monoid operation can be read directly from the above equations.

Concretely, in OCaml, α^* + α^*δ(α^*↑δ)^*α^* is the sum data type:

    type delim = char
    type 'a chunk = A of 'a list 
               | C of 'a list * delim * (('a list) * delim) list * 'a list

The monoid is then

    let wap_monoid =
    {zero = A [];
     op = fun x y -> match (x,y) with
     | (A lx,A ly) -> A (lx @ ly)
     | (A lx,C (l,d,c,r)) -> C (lx @ l, d,c,r)
     | (C (l,d,c,r),A ly) -> C (l,d,c, r @ ly)
     | (C (lx,dx,cx,rx),C (ly,dy,cy,ry)) -> 
                  C (lx,dx, (cx @ [(rx@ly),dy] @ cy), ry)}

It is a complicated version of sequence concatenation, taking into account that α^* subsequences flanked by delimiters must be maximal.

The monoid represents an (α+δ)^* subsequence: the empty subsequence is represented by A [], the singleton subsequence by A [c] or C ([],d,[],[]) (depending on if its element is a digit c or a delimiter d). The monoid representing the whole input is found by map-reduce:

    adv22_input |>
    map_reduce (function '0'..'9' as c -> A [c] | d -> C ([],d,[],[]))
      wap_monoid

resulting in

    C (['1'; '0'; '0'], ',',
       [(['2'; '0'; '0'], ','); (['3'; '0'; '0'], '|'); ...],
       ['1'; '0'; '0'; '0'])

We may then extract the the maximal delimiter-free subsequences:

    let of_chunk : 'a chunk -> ('a list) list = function
        | A x -> [x]
        | C (l,d,c,r) -> [l] @ (List.map fst c) @ [r]

The warm-up task is hence solved by map-reduce, followed by of_chunk to extract the result.

The warm-up example converted the input list of digits and delimiters to ((char list)) list, with two occurrences of list. We now replace them with arbitrary monoids.

Vertical monoid composition

The just introduced algebraic background lets us easily generalize the warm-up example to solve arbitrary nested grouping-aggregation problems. The warm-up example, taken broadly, is converting a sequence (α+δ)^* -- of α items interspersed with delimiters δ -- to a sequence (α^*↑ + δ')^*, that is, a sequence of maximal delimiter-free subsequences α^*↑ also interspersed with, generally fewer delimiters δ'. (Group delimiters are preserved but the delimiters that only serve to separate items within a group disappear after aggregation.)

We now generalize from sequences to arbitrary monoids. First, suppose there is an observable monoid M1 which reduces a sequence α^*↑ to some value ρ (e.g., a sequence of digits to an integer). The first generalization is the reduction from (α+δ)^* to (ρ+δ')^*. Now, suppose there is a monoid M2 that reduces (ρ+δ')^*. This monoid should accept not just ρ items but also delimiters. It is worth then introducing delimiter-accepting monoids

    module type obs_del_monoid = sig
      include obs_monoid
      type del                           (* type of delimiters *)
      val del : del -> t
    end

Given such M1 and M2 we can build an obs_del_monoid that reduces (α+δ)^* in one pass:

    module ChunkMonoid(S:sig
       module M1: obs_monoid
       module M2: obs_del_monoid with type ix = M1.obs
       val del_ign : M2.del -> bool
      end) : 
        (obs_del_monoid with 
         type del = S.M2.del and type ix = S.M1.ix and type obs = S.M2.obs) = struct
      open S
      type t = A of M1.t | C of M1.t * M2.t * M1.t    (* chunk *)
      let zero = A M1.zero
      let m1_to_m2 : M1.t -> M2.t = fun x -> M2.gen (M1.obs x)
      let op x y = match (x,y) with
      | (A lx,A ly) -> A (M1.op lx ly)
      | (A lx,C (l,c,r)) -> C (M1.op lx l, c, r)
      | (C (l,c,r),A ly) -> C (l,c, M1.op r ly)
      | (C (lx,cx,rx),C (ly,cy,ry)) -> 
          C (lx, M2.op cx (M2.op (m1_to_m2 (M1.op rx ly)) cy), ry)
      type ix = M1.ix
      let gen x = A (M1.gen x)
      type del = M2.del
      let del c =
        if del_ign c then A M1.zero
        else C (M1.zero, M2.del c, M1.zero)
      type obs = M2.obs
      let obs = function
        | A x -> m1_to_m2 x |> M2.obs
        | C (l,c,r) -> M2.(op (m1_to_m2 l) (op c (m1_to_m2 r)) |> obs)
    end

Here, del_ign is the predicate for the ignored delimiters. This is the vertical monoid composition. When the monoids M1 and M2 are efficient -- their operations take constant time and space -- so is the ChunkMonoid.

ChunkMonoid produces an obs_del_monoid, which can be composed again, and again. The full Advent of Code 2022 problem is a reduction with the monoid

    struct
     include
      ChunkMonoid (struct
          module M1 = struct include Horner
                      type ix = char
                      let gen x = (Char.code x - Char.code '0') |> gen end
          let del_ign = Fun.const false
          module M2 = ChunkMonoid (struct
              module M1 = AsIs(Sum)
              let del_ign = ((=) ',')
              module M2 = struct
               include AsIs(Max)
               type del = char
               let del c = assert (c='|'); zero
               end
             end)
        end)
     let gen = function '0'..'9' as c -> gen c | c -> del c
    end

We have thus demonstrated nested grouping-aggregation as a monoid reduction with an efficient monoid. Not only can it be done in one pass, in constant space and linear time (in the size of input). It can also be done embarrassingly in parallel, with linear speed-up.

The problem presented here is reminiscent of the one described by Steele in his ICFP talk. He also used the idea of `chunks' and segments. However, our problem here is the triple nesting of the chunked processing -- three times more complex -- and deriving the general vertical monoid composition.

References

Nested Grouping-Aggregation
Introducing the problem and its sequential solution. Uncannily, it also used monoids and monoid composition, albeit not explicitly realized.

talk_FP_Madrid.ml [16K]
Complete OCaml code with various monoids and monoid compositions (including nested-aggregation composition)

par_nested_agg.ml [10K]
Complete multicore OCaml code for parallel nested-aggregation, using OCaml domains. For benchmark we take a 4.6GB file (of the same structure as the sample Advent of Code file). Each core gets a slice of the file to reduce; the main thread then collects the results. The parallel code needs no locking and is assuredly race-free.

This section casts the monoid and monoid reductions (specifically, the earlier Vertical monoid composition) into an alternative, imperative form -- suitable for generating high-performance, low-level first-order imperative code for shared-memory multiprocessing. (The earlier functional interface is more suitable for message-passing, such as MPI.)

For illustration, we generate OCaml code, using MetaOCaml. The generated OCaml code, however, is deliberately low level: no algebraic data types, no pattern-matching, and certainly no first-class functions. It may then be offshored to C.

As repeatedly emphasized, monoid reductions may be done in several ways. One is sequential: left-fold, or the accumulating for-loop. It needs only a subset of the observable monoid interface:

    module type state_machine = sig
      type t
      val zero : t
      type ix
      val opgen : t -> ix -> t
      type obs
      val obs : t -> obs
    end

This is the interface for a state machine: with the state t; initial state zero; input alphabet ix; state transition opgen computing a new state from the current state and current input; and the observation obs. This interface can be cast into an imperative form:

    module type upd_state_machine = sig
      type u
      type unt
      val newupd : (u -> unt) -> unt
      type ix
      val updgen  : u -> ix -> unt
      type obs
      val obs : u -> (obs -> unt) -> unt 
    end

The type u is the type of the updateable state (e.g., reference cell); newupd allocates/initializes the state, updgen updates it based on the current input, and obs observes it. The type unt is the type of `actions' so to speak; it is also a monoid, letting us compose actions. For example, unt can be unit (the type of imperative statements) or unit code (the code of imperative statements).

A monoid is more than just a state machine: monoid reductions may also be done slice-wise, so to speak. Therefore, monoid not only provides opgen : t -> ix -> t to update the state t based on the current input ix, but also op : t -> t -> t to concatenate, or merge the (partially accumulated) states. The imperative version of the observable monoid therefore provides operations to extract the current state from the updateable state u and update u with the imported state t:

    module type upd_monoid = sig
      type u
      type unt
      val newupd : (u -> unt) -> unt
      type ix
      val updgen  : u -> ix -> unt
      type obs
      val obs : u -> (obs -> unt) -> unt 
      
      type t                                (* exportable state *)
      val from_upd : u -> t                 (* export u         *)
      val upd   : u -> t -> unt             (* update u         *)
      val set   : u -> t -> unt
      val reset : u  -> unt
    
      type us                               (* plural of u *)
      type usix
      val newupds : usix -> (us -> unt) -> unt
      val us_set : us -> usix -> t -> unt
      val us_get : us -> usix -> t
    end

The slice-wise processing can be done independently, in parallel if possible. Therefore, the interface provides us for creating an independent updateable state for each of the workers (core). The type us and its operations hence support parallel-for (which we emulate in multicore OCaml using domains).

See below for the source code that implements several staged upd_monoids, including the vertical monoid composition. See also below for the generated OCaml code -- which is indeed first-order, low-level and imperative -- essentially, C in OCaml notation.

Overall, the code generation and update semantics consistently deliver 4x single-core speed-up, scaling to multiple cores.

Version

The current version is December 2025

References

spar_nested_agg.ml [25K]
The staged version of par_nested_agg.ml, also introducing the update-oriented monoids.

pa_na_red.ml [23K]
The generated code for parallel nested aggregation. It takes the file descriptor and the number of cores to use.

par_na_main.ml [2K]
The main function to run the generated code, also showing the benchmark results.

In his ICFP'09 keynote, Guy Steele exhorted us to use reduce (or, map_reduce) rather than fold, as far as possible. Unlike fold, reduce does not commit us to a particular evaluation strategy. It can be performed sequentially, embarrassingly parallel, or in a tree-like fashion. Like the Σ notation in Math, it specifies what should be summed up, but not how or in which sequence.

Guy Steele conclusions apply to the present article as well. Especially his final thoughts:

• Associative combining operators are a VERY BIG DEAL
• Inventing (and proving) new combining operators is a very, very big deal