6.3. Parallel Reductions

rustyml::math::reduction does not ship a sum() or a mean() function. It ships 2 generic fold combinators, det_reduce and det_reduce_range, plus 1 constant, DET_REDUCE_BLOCK. The crate builds every parallel reduction in the library on top of these 2 functions. Examples include the sum of squared errors in linear regression, the global gradient norm for clip-by-global-norm, and the one-pass Welford moments in standardization. Other examples include the k-means inertia and the logistic log-loss. These functions exist to solve a problem that ordinary parallel summation cannot: a result that does not depend on the thread count.

6.3.1. What the module exposes

The public surface has 3 items. The module sits behind the math feature. Every other feature (machine_learning, neural_network, utils, metrics) pulls in the math feature. So these 3 items are always available when RustyML compiles (see 1.2. Installation and Feature Flags).

The full trait bounds appear below. The compiler error messages about them are hard to read without this context.

pub fn det_reduce<T, A, F, M>(slice: &[T], parallel: bool, fold_block: F, merge: M, identity: A) -> A
where
    T: Sync,
    A: Send,
    F: Fn(&[T]) -> A + Sync + Send,   // serial fold over 1 block
    M: Fn(A, A) -> A,                 // combines 2 partial results
{ /* ... */ }

fold_block reduces 1 block to a partial result of the accumulator type A. merge combines 2 partial results. identity is the value returned for an empty input, and it also seeds the final combine.

fold_block must be Sync + Send, because rayon can call it from any worker. This bound applies even when you pass parallel = false, since the bounds sit on the type, not on the flag. merge needs neither bound, because it always runs on 1 thread, in block order.

A can be anything Send: a scalar, a tuple such as (sum, sum_of_squares), a Welford triple, or an array of per-bucket sums.

det_reduce_range runs the same algorithm over an index range instead of a slice. Use it for reductions that read several arrays at once, or that index rows of a matrix. Its fold_block receives a Range<usize> instead of a &[T].

The module ships no sum wrapper on purpose. Below the parallel threshold, a 1-line slice.iter().sum() is already the right tool. Above the threshold, the caller almost always wants to fuse a map into the same pass. Examples include a square, an exp, or a distance function. Fusing beats building a separate intermediate array. The fold interface, not a fixed reduction, keeps that fusion at the call site.

6.3.2. Why naive parallel summation is non-deterministic

Floating-point addition is not associative. (a + b) + c and a + (b + c) can round to different f64 values. This is not a hardware bug. It is the definition of rounding to 53 bits after every operation. A sum that runs left to right on 1 thread has a fixed order, so the result is reproducible. Parallel execution removes that fixed order.

A bare slice.par_iter().sum::<f64>(), or fold().reduce(), splits the work adaptively. Rayon’s work-stealing scheduler decides which worker folds which sub-range. It also decides the order in which the partial sums combine.

A run on a machine with 4 idle cores produces 1 grouping. The same input with RAYON_NUM_THREADS=1 produces another grouping. 2 runs on a busy 16-core machine can disagree, because a thread got preempted at a different moment. Every one of these results is a correct sum of the same numbers. They only round differently, typically in the last few ULPs.

For a lot of numerical code, this jitter is harmless. For a machine-learning library, it is corrosive. A loss value that wobbles in its low bits can make an early-stopping check fire on a different iteration across runs. A gradient norm that depends on the thread count can make clip-by-global-norm clip each run slightly differently. It can also make 2 runs of the same fit produce 2 different models.

Reproducibility is a first-class promise in RustyML (see 7.1. Reproducibility and Random Seeds). A scheduler-dependent reduction breaks that promise, no matter how carefully you seed the RNG.

6.3.3. The blocked algorithm

The fix takes the grouping away from the scheduler and fixes it to a constant. det_reduce cuts the input into fixed DET_REDUCE_BLOCK-element chunks. It folds each chunk serially with your fold_block. It collects the per-block partial results in block order. Then it merges them left to right with your merge. The 2 paths, parallel and serial, differ only in how the blocks run:

if parallel {
    let parts: Vec<A> = slice.par_chunks(DET_REDUCE_BLOCK).map(fold_block).collect();
    parts.into_iter().fold(identity, merge)          // merge in block order
} else {
    slice.chunks(DET_REDUCE_BLOCK).map(fold_block).fold(identity, merge)
}

The key detail is that rayon’s par_chunks(...).collect::<Vec<_>>() is an indexed parallel iterator. No matter how work-stealing distributes the blocks across threads, the collected Vec comes back in the original block order. So the reduction tree is a pure function of the input length and DET_REDUCE_BLOCK. That tree defines which elements land in which block, and in what order the blocks combine. The tree does not depend on the thread count, on scheduling, or on the parallel flag. Both paths fold the same 16 384-element blocks in the same order, and they combine the blocks the same way.

det_reduce_range does the identical thing over n.div_ceil(DET_REDUCE_BLOCK) index blocks. Block b covers the range b * BLOCK .. ((b + 1) * BLOCK).min(n).

This makes the parallel argument a pure performance hint. It never changes which numbers get added in which order. It only decides whether the blocks run on the rayon pool or in a plain sequential loop. The crate’s own tests check this with a bitwise == comparison, not an epsilon, across empty, sub-block, exactly-1-block, and ragged multi-block lengths. So on a given build, the 2 paths are bit-for-bit identical. Varying RAYON_NUM_THREADS cannot change the result.

The module docs still note that results are not always bit-for-bit reproducible. That caveat covers cross-machine and cross-build differences. Examples include a different libm sin, an FMA contraction toggled by a different target CPU, or a different rounding inside your own fold_block. It does not cover the thread count, which the blocking fixes completely.

DET_REDUCE_BLOCK is 16 384 because that size sits near the measured throughput plateau. On a 4.2 million element f64 sum of squares, the measured speedup rises from about 14x at a 1024-element block. It peaks at about 18x at a 32 768-element block. It then falls to about 15x at a 65 536-element block. It falls further, to about 11x, at a 262 144-element block. Too few blocks remain there to balance across cores.

The equivalent f32 benchmark, with an f64 accumulator, peaks higher, at about 21x, with a 65 536-element block. 16 384 sits close to both peaks: about 3 percent under the f64 peak, and about 8 percent under the f32 peak. It works well for both element types without a separate constant for each.

The constant counts elements, not bytes, and every element type shares it. A block of 16 384 f32 values is 64 KB. A block of 16 384 f64 values is 128 KB. Both sizes sit comfortably on the plateau for their element type.

The block size defines the grouping, so it is part of the reproducibility surface. Changing it changes the deterministic result in its low bits. For that reason, DET_REDUCE_BLOCK is a const, not a runtime knob.

6.3.4. Accuracy is a side effect, not the goal

Blocking was chosen for determinism. It also improves accuracy as a side effect. This applies to both paths. The serial path also chunks into blocks, so even parallel = false is not a naive whole-array left fold.

The worst-case rounding error of summing n floats left to right grows linearly in n. It follows roughly (n - 1) * eps * S, where eps is the machine epsilon and S is the sum of the absolute input values.

det_reduce uses a 2-level scheme. Each block folds b = 16 384 terms serially. Then the partial results from ceil(n / b) blocks fold serially. The error bound becomes roughly (b + n / b) * eps * S.

For a 4.2 million element sum, that bound is about (16 384 + 256) * eps. The naive bound for the same sum is about 4.2 million * eps. The blocked bound is about 250x tighter in the worst case. The improvement holds whether the blocks ran in parallel or in sequence.

This scheme is flat blocking plus a serial merge, not a full pairwise (log n) summation tree. In practice, the accumulator width matters more than the tree shape. global_grad_norm reduces f32 gradients into an f64 accumulator inside fold_block, so that squared-gradient sum stays in f64 end to end. This accumulator choice improves accuracy more than the blocking does. Blocking sets the deterministic floor. A wide accumulator is the right choice when accuracy itself is the concern.

6.3.5. When the parallel path turns on

det_reduce does not decide parallel. The caller passes that flag. Inside the crate, a calibrated size gate produces that boolean value. Below roughly 1 block, there is nothing to parallelize. An input shorter than 16 384 elements is a single block, so forking it onto rayon only adds join overhead.

The gates live in rustyml::tuning::reduction. Each gate is shared per cost class, rather than defined per call site:

Every call site follows the same pattern: it compares a work metric against gate() and passes the result as the flag. Most sites use slice.len() for that metric. The k-means centroid accumulation instead uses n_samples * n_features, since that product is the real element count its blocked fold walks. The gate moves the crossover point, but it never touches correctness, because the blocked fold gives the same answer on either side.

Each gate has a matching setter, such as set_sum_f64 or set_sq_sum_f32, to tune the crossover for different hardware. 7.3. Performance Tuning and Parallelism covers the mechanics and the calibration process. The exp-reduction gate sits lowest, at 32 768, because each element there pays for an exp and an ln. Parallelism amortizes sooner there than for a plain add.

6.3.6. Using them in your own code

The minimal call computes a fused sum of squares over a Vec<f64> and stays serial:

use rustyml::math::reduction::det_reduce;

fn main() {
    let data: Vec<f64> = (0..1_000).map(|i| (i as f64).sin()).collect();
    let sum_sq = det_reduce(
        &data,
        false, // performance hint: small input, stay serial
        |block| block.iter().map(|&x| x * x).sum::<f64>(),
        |a, b| a + b,
        0.0,
    );
    println!("sum of squares = {sum_sq}");
}

2 details need attention here. First, det_reduce takes a &[T], so the data must be a contiguous slice. An ndarray array yields a contiguous slice only through as_slice(), and that method returns None for a view that is not in standard layout. The crate’s own idiomatic pattern reduces through det_reduce on the contiguous fast path, and it falls back to ndarray’s serial kernel otherwise. This pattern gates the flag on the size class:

use ndarray::Array1;
use rustyml::math::reduction::det_reduce;
use rustyml::tuning::reduction::get_sum_f64;

fn main() {
    let v: Array1<f64> = (0..10_000).map(|i| i as f64).collect();
    let sum = match v.as_slice() {
        Some(slice) => det_reduce(
            slice,
            slice.len() >= get_sum_f64(),
            |block| block.iter().sum::<f64>(),
            |a, b| a + b,
            0.0,
        ),
        None => v.sum(), // non-contiguous: ndarray's serial fold
    };
    println!("sum = {sum}");
}

Second, the accumulator does not have to be a scalar. That is the reason the module exposes the fold instead of a fixed reduction. A single pass can return both the sum and the sum of squares, enough for a mean and a variance. Use a tuple accumulator with a tuple merge:

use rustyml::math::reduction::det_reduce;

fn main() {
    let data: Vec<f64> = (0..10_000).map(|i| (i as f64).sin()).collect();
    let (sum, sum_sq) = det_reduce(
        &data,
        false,
        |block| block.iter().fold((0.0f64, 0.0f64), |(s, sq), &x| (s + x, sq + x * x)),
        |(sa, sqa), (sb, sqb)| (sa + sb, sqa + sqb),
        (0.0, 0.0),
    );
    let n = data.len() as f64;
    let mean = sum / n;
    let variance = sum_sq / n - mean * mean;
    println!("mean = {mean}, variance = {variance}");
}

Some reductions need to read more than 1 array at once, for example a dot product, a distance accumulation, or a per-row pick. For these cases, use det_reduce_range and index inside the block:

use rustyml::math::reduction::det_reduce_range;

fn main() {
    let xs: Vec<f64> = (0..5_000).map(|i| i as f64).collect();
    let ys: Vec<f64> = (0..5_000).map(|i| (i as f64).cos()).collect();
    let dot = det_reduce_range(
        xs.len(),
        false,
        |range| range.map(|i| xs[i] * ys[i]).sum::<f64>(),
        |a, b| a + b,
        0.0,
    );
    println!("dot = {dot}");
}

det_reduce versus ndarray’s .sum()

ndarray’s .sum(), .dot(), and .mean() are serial and single-threaded. Their internal grouping is their own, so it will not generally match det_reduce’s blocking bit-for-bit. For small arrays, ndarray is the right choice. It is shorter to write, it needs no closures, and below the gate det_reduce would run serially anyway, with more setup code.

Use det_reduce when 3 conditions hold. The buffer is large and contiguous. The reduction should run in parallel. The parallel result must be reproducible. ndarray does not offer that combination. Even with ndarray’s rayon feature, a bare parallel sum is scheduler-dependent.

det_reduce is also the right tool when you want to fuse a map into the reduction, or to accumulate something richer than a scalar. Use ndarray by default, for convenience and small data. Switch to det_reduce at the point where you were about to write par_iter().sum() and you need the answer to stay stable. See 1.3. Working with ndarray for the interop details. See 4.2. Standardization and Normalization for a real Welford reduction built on this fold.

6.3.7. Verifying determinism across thread counts

You can test this claim from outside the crate. The following program reduces 2 million values on the rayon path and prints the sum at full precision:

use rustyml::math::reduction::det_reduce;

fn main() {
    let data: Vec<f64> = (0..2_000_000).map(|i| (i as f64 * 0.001).sin()).collect();
    let sum = det_reduce(
        &data,
        true, // force the parallel path
        |block| block.iter().sum::<f64>(),
        |a, b| a + b,
        0.0,
    );
    // Full-precision print so a low-bit difference would show
    println!("{sum:.17e}");
}

Build the program once. Then run it under different thread counts by setting RAYON_NUM_THREADS, which caps rayon’s global pool:

RAYON_NUM_THREADS=1 ./target/release/demo
RAYON_NUM_THREADS=2 ./target/release/demo
RAYON_NUM_THREADS=8 ./target/release/demo

All 3 runs print the identical 17-digit mantissa. Every run folds the same 16 384-element blocks and merges them in the same block order, regardless of how many workers carried them. Setting the flag to false does not change the output either. A version that swaps the body for data.par_iter().sum::<f64>() behaves differently. Under different RAYON_NUM_THREADS values, and on a large enough input, that version prints sums that differ in their last digits. This is the exact failure mode that det_reduce closes.

The related numeric primitives share the same determinism discipline. See 6.1. Distance Metrics, 6.2. Matrix Multiplication, and the broader 6.0. Math Utilities overview.