# 6.3. Parallel Reductions
`rustyml::math::reduction` does not ship a `sum()` or a `mean()` function. It ships 2 generic
fold combinators, [`det_reduce`](https://docs.rs/rustyml/latest/rustyml/math/reduction/fn.det_reduce.html)
and [`det_reduce_range`](https://docs.rs/rustyml/latest/rustyml/math/reduction/fn.det_reduce_range.html),
plus 1 constant, [`DET_REDUCE_BLOCK`](https://docs.rs/rustyml/latest/rustyml/math/reduction/constant.DET_REDUCE_BLOCK.html).
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](../Chapter-01/1.2._Installation_and_Feature_Flags.md)).
| `DET_REDUCE_BLOCK` | `pub const DET_REDUCE_BLOCK: usize = 16_384` | Fixed block size, in elements, that sets the grouping |
| `det_reduce` | `fn det_reduce<T, A, F, M>(slice: &[T], parallel: bool, fold_block: F, merge: M, identity: A) -> A` | Folds a slice in fixed blocks |
| `det_reduce_range` | `fn det_reduce_range<A, F, M>(n: usize, parallel: bool, fold_block: F, merge: M, identity: A) -> A` | Folds the index range `0..n` in fixed blocks |
The full trait bounds appear below. The compiler error messages about them are hard to read
without this context.
```rust,ignore
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](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md)).
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:
```rust,ignore
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:
| `get_sum_f64` | 262 144 | `f64` sum-style reductions (SSE, Welford moments, k-means inertia) |
| `get_sq_sum_f32` | 65 536 | `f32` to `f64` square-sum for clip-by-global-norm |
| `get_scan_f64` | 262 144 | short `f64` per-row scans (arg-min, distance scans) |
| `get_exp_reduce` | 32 768 | the exp-heavy logistic log-loss reduction |
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](../Chapter-07/7.3._Performance_Tuning_and_Parallelism.md)
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:
```rust
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:
```rust
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`:
```rust
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:
```rust
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](../Chapter-01/1.3._Working_with_ndarray.md) for the interop details.
See [4.2. Standardization and Normalization](../Chapter-04/4.2._Standardization_and_Normalization.md)
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:
```rust
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:
```bash
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.1._Distance_Metrics.md), [6.2. Matrix Multiplication](./6.2._Matrix_Multiplication.md),
and the broader [6.0. Math Utilities](./6.0._Math_Utilities.md) overview.