7.3. Performance Tuning and Parallelism

RustyML parallelizes with rayon, and only rayon. It does not link BLAS, and it does not use OpenMP. Every kernel in the crate runs on rayon’s global pool, except one.

The matrix-product backend, gemmkit, is the exception. gemmkit keeps a small number of persistent, private rayon pools, each sized to a fraction of the machine width. It installs the smallest pool that still holds every worker a product needs, so the fork-join step sees no idle slack.

These pools are still rayon pools. gemmkit builds each one once and reuses it warm. gemmkit skips its own pools when it already runs on another rayon pool. It falls back to the ambient pool when no tier is wide enough.

So the claim “rayon and only rayon” stays true. The claim “no thread pool of its own” does not. Every hot loop that can spread across cores does so, but none of them forks rayon without a check first.

Each parallel kernel sits behind a gate: a work estimate compared against a calibrated threshold. A small input stays on 1 core. The tuning module moves these gates at runtime. The shipped defaults may not fit your CPU. Recalibrate them without guesswork.

If you know scikit-learn’s n_jobs or Keras’ intra_op_parallelism_threads, note a difference. RustyML has no single knob. The crossover between “serial is faster” and “rayon is faster” differs by kernel and by element width.

7.3.1. The parallelism model: rayon behind a gate

Rayon is not free. Handing a task to the pool costs a fork, a join, and, for a reduction, some collecting. For a large matrix product or a million-element exp map, this overhead disappears into the work. For a 32x32 GEMM inside 1 LSTM timestep, or a ReLU over a few thousand activations, the fork/join step is the runtime. There the parallel version loses to a single core running at memory bandwidth.

So every gated kernel first estimates its own work. A convolution counts FLOPs. A map or a reduction counts elements. A tree walk counts node visits. Pooling counts window taps. Only when that estimate clears a threshold calibrated to the crossover point does the kernel use rayon.

RustyML does not gate the matrix product itself. Scheduling there belongs entirely to gemmkit. 6.2. Matrix Multiplication covers the backend and the small amount the crate adds to it. The shape of gemmkit’s decision differs from RustyML’s own gates.

First, a work gate, parallel_threshold, defaults to 48 * 48 * 256 = 589_824. It compares the m * n * k product, not a FLOP count. There is no factor of 2 to remember. Below the gate, the product runs on 1 thread, no matter how many workers the caller requested.

Above the gate, the worker count ramps with the total work, not with any single dimension. gemmkit assigns 1 worker per par_mnk_per_worker (2,000,000 by default) of m * n * k. The count is floored at 1 and capped by the core count and the job count. Exact crossover numbers are machine-specific. Inspect the current values through tuning::matmul::backend, or measure your own machine with the gemmkit-tune autotuner (see 7.3.5).

gemmkit then snaps the ramped width to one of the persistent pool tiers described above. Rayon’s fork-join tax scales with a pool’s slack, its width minus the workers doing real work. An 8-worker product in an exact 8-wide pool beats the same product in the 32-wide global pool by a wide margin.

A matvec shape (m == 1 or n == 1) leaves this path for a dedicated bandwidth-bound one. It stays serial below a byte floor derived from the machine’s L2 cache size, then splits across a worker cap sized to memory bandwidth.

RustyML’s own kernels all run on the global rayon pool, so they compose safely when nested. A parallel outer loop that calls into a gated convolution or reduction does not oversubscribe the machine. The inner work nests into the same pool instead of spawning a second wave of threads.

The matrix product behaves the same way from the outside. gemmkit checks whether it already runs on a rayon worker. If so, it stays in the caller’s pool instead of installing a private tier. This is also why the crate’s internal call sites can force a product serial (Parallelism::Serial) when they already sit inside a parallel region. That choice avoids forking twice. It has nothing to do with correctness.

7.3.2. What a gate changes, and what it never changes

The following contract from the tuning module docs makes tuning safe:

A gate only picks an execution strategy. It never changes what the code computes. The elementwise and reduction gates give the same result serial or parallel. The matrix-product scheduling lives in the gemmkit backend (see the matmul submodule). Its results reproduce on the same machine for a fixed configuration regardless of worker count. The matmul gates kept here only shape caller-side tiling. Retuning a gate does not change any result.

That single paragraph holds 2 different guarantees. The elementwise maps (ReLU, sigmoid, scaling, normalization) are embarrassingly parallel. Each output element is independent. Serial and parallel give bit-identical results, so moving the gate cannot change a single bit.

The reductions (sums of squares, Welford moments, clip-by-global-norm) are trickier. Floating-point addition is not associative. A naive rayon sum would group its partial sums by work-stealing, and give a different rounding on every run. RustyML avoids this with the deterministic blocked fold in crate::math::reduction. The fold cuts the input into fixed-size blocks. The grouping depends only on a compile-time block size, never on the thread count or the gate.

So above a reduction gate, the parallel path still matches the serial result. Moving the gate never changes it.

The matrix product used to make a weaker promise. The old hand-rolled row split changed the block height with the thread count, and so changed the summation order. It no longer does. gemmkit’s blocking and job order do not depend on the worker count. For a fixed machine and configuration, the same product reproduces the same result bit-for-bit, no matter how many threads ran it. Re-running it is also deterministic.

The crate’s tests check this claim, at the following scope. For f64, a forced-serial product is compared on to_bits() against the same product forced onto 2, 4, 8, 16, and 32 workers. The check covers a square shape, a shallow-k shape, and a deep-k shape. For f32 the sweep is narrower: 2, 4, and 32 workers over 2 shapes.

The matvec check is narrower still. It compares forced-serial against automatic scheduling on 1 shape, with no worker sweep. A separate test compares a fused bias+ReLU epilogue bitwise against a plain product followed by the same scalar map.

Together these tests are good evidence for the guarantee. They are not an exhaustive proof of it, and the guarantee itself belongs to the backend, not to the tests. The matmul path is no longer the weak link. On 1 machine, it reproduces as well as the reductions do.

The caveat is the “fixed machine and configuration” clause. A different CPU picks a different SIMD width. A changed backend knob can change the blocking. Cross-machine bit-equality is still not a promise.

The deterministic reductions in 6.3 stay the stricter guarantee. They give the same answer by construction. They do not depend on the backend happening to block the same way on 2 machines.

The practical result stays the same: you can retune any gate freely, and your model’s outputs will not move. A gate is a performance knob, never a numerical one. For exact reproducibility across machines, see 7.1. Reproducibility and Random Seeds instead.

7.3.3. The tuning module: the runtime override surface

Every gate is a process-global AtomicUsize, initialized to its calibrated default. The tuning module is the facade that makes every gate discoverable. For each gate, it exposes a set_*(usize) function and a get_*() -> usize function, grouped into submodules by kernel family. Setters and getters are plain functions. Read the defaults your build shipped with like this:

use rustyml::tuning;

fn main() {
    // RustyML's own matmul knobs: the caller-side tiling policy, and nothing else.
    println!("chunk elems:   {}", tuning::matmul::get_chunk_elems());
    println!("cache bytes:   {}", tuning::matmul::get_cache_resident_max_bytes());

    // The product's scheduling knobs belong to the backend, reached through the alias.
    println!("mnk gate:      {}", tuning::matmul::backend::parallel_threshold());
    println!("mnk/worker:    {}", tuning::matmul::backend::par_mnk_per_worker());
    println!("pool tiers:    {}", tuning::matmul::backend::pool_classes());

    // Elementwise maps and deterministic reductions (element counts).
    println!("exp map f32:   {}", tuning::elementwise::get_exp_map_f32());
    println!("cheap map f64: {}", tuning::elementwise::get_cheap_map_f64());
    println!("sum f64:       {}", tuning::reduction::get_sum_f64());
    println!("exp reduce:    {}", tuning::reduction::get_exp_reduce());

    // Tree walks, conv/pool engines, normalization, metrics.
    println!("tree visits:   {}", tuning::tree::get_traversal_min_visits());
    println!("conv flops:    {}", tuning::conv::get_parallel_min_flops());
    println!("pool ops:      {}", tuning::pool::get_parallel_min_ops());
    println!("gn param grad: {}", tuning::norm::get_gn_param_grad());
    println!("silhouette:    {}", tuning::metrics::get_silhouette());
}

The matmul submodule has a distinct shape. It owns 2 gates of its own, plus backend, a pub use gemmkit_ndarray::tuning re-export. Every GEMMKIT_* knob is reachable through a set_*/getter pair, with no direct gemmkit dependency in your Cargo.toml. The re-export forwards through the adapter on purpose: the knobs are process-global atomics, so a set_* call made on a separately resolved second gemmkit would write a copy that the adapter never reads. The getters under backend use bare names, for example parallel_threshold(), not the get_-prefixed style. That naming is gemmkit’s own, not RustyML’s facade.

Here is the full surface, with the shipped default and the unit each gate compares against:

The normalization layers add 7 more gates under tuning::norm::: set_/get_batch_norm, set_/get_bn_col_stats, set_/get_bn_plane_stats, set_/get_ln_row, set_/get_ln_col_stats, set_/get_gn_row, and set_/get_gn_param_grad. All 7 ship at 262_144. The clustering metric adds 1 more, tuning::metrics::set_/get_silhouette, also at 262_144.

A gate whose feature is not compiled in is simply not there. tuning::conv, tuning::pool, and tuning::norm need neural_network. tuning::tree needs machine_learning. tuning::metrics needs metrics. tuning::matmul and reduction::exp_reduce need math.

The elementwise gates split in 2: the f32 half needs neural_network, and the f64 half needs machine_learning or utils. Build with only the modules you use (see 1.2. Installation and Feature Flags and 7.4. Minimal Builds and Modular Integration). The irrelevant knobs then vanish.

This list is short by design. RustyML owns 2 matmul knobs and about 2 dozen kernel gates. Every one is a threshold on a work estimate.

There is no per-dtype matrix-product gate any more. f32 and f64 GEMM used to have separate crossovers here. They no longer do, because the backend gates on m * n * k, and the element width is the backend’s business, not yours. If a matrix product runs with the wrong parallelism on your machine, look in backend, not in this table.

7.3.4. Why the shipped defaults may be wrong for your machine

Most of these defaults were measured on the maintainer’s hardware. The tuning module docs record the machine: an AMD Ryzen 9 9950X with 16 cores, 32 threads, and 64 MiB L3. The few defaults that were not measured there are worse, not better. cache_resident_max_bytes is an educated guess. pool_parallel_min_ops is deliberately parked away from its measured bracket, for a reason given in 7.3.5.

A serial/parallel crossover is not a universal constant. It is a ratio: how fast 1 core runs the kernel, against how much fork/join overhead rayon adds. Both ends move with the machine. More cores lower the per-core share of a fixed problem, so the parallel side needs a bigger problem to break even. Faster single-thread SIMD raises the serial baseline. A larger L3 keeps more of a matrix resident, and this shifts where the tiled-GEMM-versus-GEMV-swarm decision flips.

The gate tied most obviously to one machine is cache_resident_max_bytes. The source documents it as: set this to the machine’s actual shared-L3 size. The default of 64 MiB is a guess at a typical L3. The band around it has not been calibrated.

The backend’s knobs carry the same caveat, one level down, and gemmkit is direct about it. par_mnk_per_worker defaults to 2,000,000, and pool_classes defaults to 2 tiers on x86. Both defaults come from that same Zen5 9950X, 32 hardware threads over 16 physical cores. There the 2 tiers are width/4 (8 workers) and width/2 (16 workers, the physical core count).

The aarch64 arm of pool_classes defaults to 1 tier instead. That tier was measured on an M4 Max (14 cores, 10 performance and 4 efficiency, no SMT). There the 1 tier is width/2 (7 workers).

Those are 2 real machines, not a model of yours. A knob whose crossover depends on the architecture ships with a cfg(target_arch)-split default. Each side is calibrated on its own reference machine. On any third architecture, the pool tiers default to off, pending on-device validation. Measure your own machine with the gemmkit-tune autotuner (see 7.3.5), and apply the result as a GEMMKIT_* profile.

Your CPU may differ a lot from the reference machines. Examples: an 8-core laptop against a 32-core workstation, half the L3, or an Apple-silicon NEON target against AVX-512. In that case, the defaults stay in the right neighborhood but miss the optimum.

For most workloads the difference is small. The elementwise and reduction gates sit so far out that, at typical sizes, those kernels run serial regardless. The backend knobs matter only when matrix products dominate your runtime. Recalibrate after you measure that they do. Do not recalibrate on principle alone.

7.3.5. Recalibrating: the two tools, then the setters

Recalibration splits along the same line as the knobs, and each half uses a different tool. The matrix product belongs to gemmkit, so retune it with gemmkit’s autotuner. Everything else belongs to RustyML, so retune it with RustyML’s calibration bench.

To retune the product, install and run the backend’s sweeper on the target machine:

cargo install gemmkit-tune
gemmkit-tune

The sweeper runs on the machine it targets and emits a profile of GEMMKIT_* environment variables, ready to source. This is the deployment path, and it is the good one. Environment variables retune an already-built binary with no recompile, so 1 artifact can carry a different profile per host.

The programmatic equivalent is tuning::matmul::backend::set_parallel_threshold(..) and its siblings. Check the precedence rule before you use it. A knob resolves in this order: a per-call argument, then a programmatic set_* call, then a GEMMKIT_* variable, then the compiled default. A set_* call stores its value unconditionally. So once any code in the process calls a setter, the matching environment variable is dead for the rest of that process. RustyML never calls a setter on your behalf, for exactly this reason. Doing so would silently override a profile you had sourced.

A GEMMKIT_* value that fails to parse is not silently ignored either. The first access warns once on stderr and falls back to the default. It never panics.

For everything else, RustyML ships the calibration bench the maintainer used. It runs under cargo bench with harness = false, and prints straight to stdout:

# Elementwise / reduction / tree / conv / pool / normalization crossovers.
# Prints the tables and rewrites benches/calibrations/RESULTS.md.
cargo bench --bench parallel_gates

It needs the machine_learning and neural_network features. Its submodules cover the convolution engine, pooling, the elementwise and reduction kernels, the tree walks, and the normalization layers. For each kernel class, it forces the serial and the parallel implementation on either side of that class’s gate. It walks a ladder of shapes and reports the speedup at each rung.

Here is a real excerpt from the checked-in benches/calibrations/RESULTS.md. It was regenerated on 2026-07-26, on the 9950X at 32 rayon threads:

## conv engine FLOPs gate (CONV_PARALLEL_MIN_FLOPS), batch == 1

| shape | work (FLOPs) | serial (us) | parallel (us) | speedup |
|---|---:|---:|---:|---:|
| conv 3c->8f 16px k3   |     84672 |   12.1 |  25.2 | 0.48x |
| conv 8c->16f 32px k3  |   2073600 |   83.2 |  82.3 | 1.01x |
| conv 16c->32f 32px k3 |   8294400 |  109.2 |  85.0 | 1.29x |
| conv 32c->64f 64px k3 | 141705216 | 1437.6 | 298.2 | 4.82x |

**Takeaway:** crossover between 2073600 and 8294400 FLOPs.

Read the table the way its takeaway line does. The gate constant belongs at the crossover, with a small safety margin toward serial. A 2x penalty on a small tensor, to win 5% on a medium one, is a bad trade. The ladder makes this asymmetry visible. At 84,672 FLOPs, the parallel path runs at half the speed of serial. At 8.3M FLOPs, it is only 1.29x faster.

There is a trap here. The bench reports 1 work-estimate number per rung, but the serial cost of a kernel is not always a function of that number alone. POOL_PARALLEL_MIN_OPS stayed at 12,000 on purpose, even though the tool’s own bracket for it said 25K to 49K window taps.

Serial pooling speed depends strongly on the channel count. Window geometry gets amortized across channels. So the 25,088-tap 1x28x28x32 rung runs 14.4 us serial, while the smaller 16,384-tap 64x16x16x1 rung runs 76.7 us. Tap count alone does not predict serial cost.

Following the bracket would have pushed that 16K-tap shape back to serial. That would have given up a measured 45 us saving, only to avoid the 7 us loss the 12,288-tap 1x64x64x3 shape currently pays. Calibration output is evidence, not an instruction. Check which rungs of the ladder your workload actually sits on before you move a constant to match a bracket.

Apply the numbers at program start, before any real work touches a gate. The atomics are process-global, and every kernel reads them live:

use ndarray::{Array1, Array2};
use rustyml::machine_learning::LinearRegression;
use rustyml::tuning;

fn main() {
    // RustyML's own gates, from a `parallel_gates` run on this machine.
    tuning::matmul::set_cache_resident_max_bytes(32 * 1024 * 1024); // this CPU's real L3
    tuning::conv::set_parallel_min_flops(6_000_000);
    tuning::pool::set_parallel_min_ops(30_000);
    tuning::reduction::set_scan_f64(131_072);

    // The backend's scheduling knobs, only if you are not shipping a GEMMKIT_* profile:
    // a setter shadows the matching env var for the rest of the process.
    tuning::matmul::backend::set_parallel_threshold(2_000_000);
    tuning::matmul::backend::set_par_mnk_per_worker(4_000_000);

    // Each store is one relaxed atomic write. Read it back to confirm it took.
    assert_eq!(tuning::conv::get_parallel_min_flops(), 6_000_000);
    assert_eq!(tuning::matmul::backend::parallel_threshold(), 2_000_000);

    // Same API, same numbers. Only the serial/parallel strategy shifted.
    let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
    let y = Array1::from_vec(vec![3.0, 5.0, 7.0, 9.0, 11.0]);
    let mut model = LinearRegression::new(true);
    model.fit(&x, &y).unwrap();
    let preds = model.predict(&x).unwrap();
    println!("prediction: {:.3}", preds[0]);
}

RustyML’s own gates have no config file and no environment variable. The override API is the whole mechanism there. GEMMKIT_* variables reach only the backend. Setters are global and last for the life of the process. Call them once, in main or in a OnceLock-guarded init, rather than per model.

To see where threading actually pays off on your machine before you touch anything, run cargo bench --bench matmul_kernels.

7.3.6. Controlling rayon itself

The gates decide whether to go parallel. Rayon decides how wide. RustyML’s own kernels use rayon’s global pool and never build one of their own, so the standard rayon controls apply to them unchanged. The exception is the backend’s private tier pools, covered 2 paragraphs down.

The simplest control is the RAYON_NUM_THREADS environment variable, read when the pool is first touched:

RAYON_NUM_THREADS=8 ./my_program

For programmatic control, build the global pool once, before any rustyml call reaches it. rustyml does not re-export rayon, so add rayon as your own dependency (rayon = "1") to do this:

fn main() {
    rayon::ThreadPoolBuilder::new()
        .num_threads(8)
        .build_global()
        .unwrap(); // build_global fails if the pool was already initialized

    // ... rustyml calls now run on an 8-thread pool ...
}

One interaction here catches people by surprise. Shrinking the ambient pool does not rescale the matrix product. gemmkit sizes its worker count from the work, m * n * k divided by par_mnk_per_worker, and caps it by the machine’s core count. It reads that core count once, from std::thread::available_parallelism(), and caches it for the process. It never consults rayon::current_num_threads().

gemmkit’s pool tiers derive from that same cached machine width. So RAYON_NUM_THREADS=8 on a 32-thread box does not narrow gemmkit. gemmkit still ramps toward 32 workers. It still snaps to tiers of 8 and 16, in its own private pools, with no regard for your 8-wide global pool.

Your pool reasserts itself only at the top of the range. Once a product is large enough to want the full 32 workers, no tier is wide enough to hold it. The work then falls back to the ambient pool, for example 32 jobs over your 8 threads. To make the product narrower, raise backend::set_par_mnk_per_worker to demand more work per worker, or source a GEMMKIT_* profile. Do not resize the rayon pool for this.

RustyML’s own gates have the mirror problem. They are fixed numbers, calibrated at 1 thread count (32, on the 9950X, per RESULTS.md). They read nothing at all about the pool. They just compare a work estimate to a constant.

If you halve the pool, the crossover where parallel starts to win moves up in practice. Each remaining core now carries a bigger share, so the fork/join overhead pays off later. Yet conv::parallel_min_flops still holds its original value. The gate is not wrong, only no longer optimal.

Recalibrate with parallel_gates under the pool size you plan to deploy, or accept that the defaults assume the calibration machine’s thread count. The rule from both halves: the pool size is a resource limit, not a tuning parameter. Nothing rescales itself when you change it.

Oversubscription is the failure mode to watch for. RustyML’s kernels nest into 1 global pool. gemmkit deliberately stays in the caller’s pool instead of installing a tier, whenever it finds itself already on a rayon worker. So rustyml calls inside a rayon region stay safe.

The danger is building 2 sources of parallelism outside that pool. One example is an OS thread pool. Another is a std::thread::spawn fan-out. A third is a second rayon pool, built with ThreadPoolBuilder::build() instead of build_global. In each case, every worker independently calls into rustyml. Now each worker’s gated kernels each try to fill the global pool, and you get threads * threads contention.

If you already parallelize at the application level, give rustyml a smaller pool. Or keep your outer fan-out serial per item, and let rustyml’s own gates spread each item’s work. Do not stack pools.

7.3.7. Why the gate reads are free: relaxed atomics

Every gate read on a hot path is a single relaxed atomic load. Every setter is a single relaxed store. There is no lock, no fence, no contention. This design is deliberate. It makes runtime-tunable gates cost nothing more than the compile-time constants they replaced.

The reasoning comes from the macro that generates them. The gate only selects a strategy and never changes a result, so it needs no stronger memory ordering. A relaxed load has no ordering obligation to satisfy. On every mainstream architecture, it compiles to an ordinary load with no barrier. So a kernel that checks flops >= conv_parallel_min_flops() before every convolution pays almost nothing for the indirection.

The flip side of Relaxed is that a set_* call does not synchronize with in-flight kernels. Say you change a gate from 1 thread while another thread is mid-computation. Then no guarantee exists about which reads see the old value, and which see the new. This is harmless, because the gate only picks serial versus parallel, and both paths compute the same result. A mid-run flip at worst makes 1 kernel pick a slightly suboptimal strategy. It never gives a wrong answer.

Still, set gates at startup rather than while under load. This keeps behavior predictable, and every kernel sees 1 consistent policy.

7.3.8. What not to tune, and a sane workflow

Reach for the gates last, not first. Before you touch any of them, measure end to end. Time the actual fit, predict, or training loop you care about, and find where the time goes.

9 times out of 10, the answer is not a mis-set gate. More often it is a feature build with more modules than you need, dragging in code you never call. Or it is an f64 pipeline where f32 would halve the memory traffic. Or a model gets refit in a loop when it could fit once. Or the problem is simply small enough that it runs serial by design, and no gate change touches it.

The elementwise and reduction gates sit especially far out. At ordinary preprocessing and layer sizes, those kernels run serial no matter what you set. Moving cheap_map_f64 does nothing for a standardization over 10,000 rows and 20 features. That is 200,000 elements, still an order of magnitude below its 4,000,000 gate.

When profiling does point at parallelism, tune in this order.

First, confirm the pool is the size you expect, with rayon::current_num_threads(). A wrong RAYON_NUM_THREADS, or an accidental second pool, dwarfs the effect of any gate.

Second, if matrix products dominate, run cargo install gemmkit-tune on the deploy machine, and source the GEMMKIT_* profile it emits. For most classical-ML and dense-network work, this step gives the biggest gain, and it costs no recompile. Set the values as an environment profile, not as backend::set_* calls, unless you have a reason to hard-code them. A setter permanently shadows the environment variable, and takes the deployment knob away from whoever runs the binary.

Third, if you lean on the tiled-product paths (KNN, t-SNE, MeanShift), set cache_resident_max_bytes to your CPU’s real shared-L3 size.

Take those 3 steps first. Only then, check whether parallel_gates shows your kernels crossing over at a size that differs from the shipped default. If it does, adjust the elementwise, reduction, tree, conv/pool, or normalization gates. Change 1 gate at a time. Re-measure the end-to-end number, and keep the change only if it helped. Retuning never changes your results, so the only cost of a wrong guess is the time spent tuning a gate that was never the bottleneck.