Performance

numeris achieves competitive performance via SIMD intrinsics, register-blocked micro-kernels, and direct formulas for small sizes.

SIMD Architecture

SIMD is always-on for f32 and f64 — no feature flag, no runtime detection. Dispatch selects the widest available ISA at compile time via #[cfg(target_feature)]. Integer and complex types fall back to scalar loops via TypeId dispatch, with zero runtime overhead (dead-code eliminated at monomorphization).

Architecture	ISA	f64 tile (MR×NR)	f32 tile (MR×NR)
aarch64	NEON (128-bit)	8×4	8×4
x86_64	SSE2 (128-bit)	4×4	8×4
x86_64	AVX (256-bit)	8×4	16×4
x86_64	AVX-512 (512-bit)	16×4	32×4
other	scalar fallback	4×4	4×4

AVX and AVX-512 require compile-time opt-in:

RUSTFLAGS="-C target-cpu=native" cargo build --release
# or explicitly:
RUSTFLAGS="-C target-feature=+avx2,+avx512f" cargo build --release

SSE2 (x86_64) and NEON (aarch64) are always-on baselines.

Matrix Multiply Micro-Kernels

The matmul micro-kernels are inspired by nano-gemm and faer (Sarah Quinones). Each kernel:

1. Tiles the output matrix into MR×NR panels
2. Loads A and B into SIMD registers
3. Accumulates the full k-sum into MR×NR SIMD accumulators (never reading/writing C in the inner loop)
4. Writes the tile to C exactly once

This reduces memory traffic by O(n) compared to a naive C-read-per-k implementation. The k-loop is cache-blocked at KC=256 elements for L1 cache locality.

For each (block_i, block_j) tile of C:
  accum[0..MR*NR] = 0
  for k = 0..K:
    load A[block_i..block_i+MR, k] into MR SIMD vectors
    load B[k, block_j..block_j+NR] into NR scalars
    for row in 0..MR: accum[row*NR..] += A_row * B_col
  store accum → C[block_i..block_i+MR, block_j..block_j+NR]

Optimizations Applied

• Direct formulas for N ≤ 4: inverse() and det() use adjugate-based closed-form expressions, bypassing LU decomposition entirely.
• Unrolled LU/Cholesky for N ≤ 6: direct data[col][row] array indexing eliminates trait dispatch in inner loops.
• AXPY SIMD: LU elimination, QR Householder, SVD bidiagonalization, and Hessenberg reduction all use axpy_neg_dispatch on contiguous column slices.
• 4-accumulator dot product: reduces dependency chains in the dot product inner loop.

Benchmark Results

Platform: Apple Silicon M-series (aarch64 NEON), Rust stable.

Compared against nalgebra 0.34 and faer 0.24. All benchmarks run with cargo bench.

Benchmark	numeris	nalgebra	faer	Winner
matmul 4×4	4.9 ns	4.9 ns	58 ns	~tie
matmul 6×6	13.4 ns	20.0 ns	87 ns	numeris
matmul 50×50 (dyn)	5.76 µs	6.63 µs	6.3 µs	numeris
matmul 200×200 (dyn)	369 µs	361 µs	193 µs	faer
dot 100 (dyn)	11.6 ns	14.5 ns	—	numeris
LU 4×4	33.2 ns	28.2 ns	203 ns	nalgebra
LU 6×6	84.7 ns	82.1 ns	292 ns	nalgebra
LU 50×50 (dyn)	8.4 µs	7.5 µs	7.7 µs	nalgebra
Cholesky 4×4	25.2 ns	11.8 ns	139 ns	nalgebra
Cholesky 6×6	70.7 ns	39.6 ns	186 ns	nalgebra
QR 4×4	46.4 ns	90.6 ns	303 ns	numeris
QR 6×6	85.5 ns	207.9 ns	445 ns	numeris
SVD 4×4	299 ns	461 ns	1278 ns	numeris
SVD 6×6	1171 ns	925 ns	1858 ns	nalgebra
Inverse 4×4	27.6 ns	23.3 ns	—	nalgebra
Inverse 6×6	163 ns	127 ns	—	nalgebra
Eigen sym 4×4	165 ns	201 ns	578 ns	numeris
Eigen sym 6×6	287 ns	528 ns	1088 ns	numeris

Summary

• numeris wins: matmul 6×6 (1.5×), QR (2×), SVD 4×4 (1.5×), symmetric eigendecomposition (1.2–1.8×), dot product, matmul 50×50
• nalgebra wins: Cholesky, LU, inverse, SVD 6×6 — Cholesky gap is a measurement artifact (see below)
• faer wins: large dynamic matmul (200×200) — A+B packing, cache-aware blocking
• faer has high overhead at small sizes due to dynamic dispatch / runtime machinery
• matmul 4×4: dead heat with nalgebra (4.9 ns each)

Cholesky benchmark artifact

The ~2× Cholesky gap is a measurement artifact, not a real performance difference. Micro-benchmarking shows raw computation is within 4% of nalgebra. The gap comes from Criterion's black_box reading Result<CholeskyDecomposition, LinalgError> byte-by-byte (48 ldrb instructions) vs nalgebra's Option<Cholesky> using word-sized ldr (17 instructions).

Remaining Opportunities

• Large matmul: A-panel packing + larger tile sizes could close the remaining ~2× gap with faer
• SVD 6×6: 27% behind nalgebra — likely dominated by Givens rotations in bidiagonal QR
• LU: small margin behind nalgebra — possibly similar Result vs Option artifact

FFT

The fft module's DynFft tier deinterleaves into structure-of-arrays re/im buffers and runs radix-2 butterflies through the same per-ISA SIMD kernel macros as the element-wise ops; the no-alloc fixed tier stays scalar (its audience is embedded, small N). Bluestein (arbitrary / prime lengths) runs its inner power-of-two transforms on the same SIMD core. Real transforms are half-size in both directions; 1D and 2D FFT convolution pad to a power of two and use the real plans; the 2D row pass runs on a cache-blocked transposed copy so both passes are contiguous.

Apple Silicon (aarch64, NEON), f64, cargo bench -p numeris-bench --bench fft, single thread. Measured on a busy machine, so treat the absolute times as approximate — the ratios are what the optimizations bought relative to the first landing (gather/scatter 2D, scalar Bluestein, full-length real inverse, Bluestein-at-exact-length convolution):

Transform	Size	Time	vs. first landing
DynFft forward, power of two	1024	~4 µs	1.2×
DynFft forward, power of two	65536	~0.9 ms	1.2×
DynFft forward, f32	4096	~12 µs	1.7×
DynFft inverse (conj folded into copies)	1024	~4 µs	1.6×
DynFft Bluestein (prime)	1009	~26 µs	2.4×
DynFft Bluestein (prime)	4093	~130 µs	2.3×
DynRealFft forward	65536	~0.6 ms	1.3×
DynRealFft inverse (half-size)	65536	~0.6 ms	2.0×
fft_convolve	4096 ∗ 512	~0.17 ms	11×
fft_convolve	10000 ∗ 1000	~0.35 ms	12×
DynFft2 forward (sequential)	512²	~3.3 ms	1.4×
DynFft2 forward (sequential)	1024²	~16 ms	1.6×
DynRealFft2 forward (sequential)	1024²	~9 ms	—

A prime length costs roughly 6× the next power of two (two length-≥ 2N transforms plus the chirp multiplies), so pad to a power of two when you control the length — fft_convolve and fft_convolve2d do this for you.

No-std Performance

On embedded targets with no hardware FPU, float operations fall back to the libm software implementation. Performance in this mode is entirely determined by the target's ALU throughput — SIMD code paths are not compiled for targets without SIMD (the TypeId dispatch compiles down to scalar loops).

Parallelism (rayon)

SIMD parallelizes within a core (vector lanes); the optional rayon feature parallelizes across cores (threads). The two are orthogonal and compose — each worker thread still runs the SIMD kernels. Parallelism is opt-in because rayon requires std and a thread pool, which the no-std / embedded baseline cannot assume:

numeris = { version = "0.5", features = ["rayon"] }   # implies std

MSRV

The rayon feature pulls in the rayon crate (Rust ≥ 1.80) — though the crate's MSRV is 1.80 regardless, so the feature does not raise it.

The feature is purely additive: builds without it are byte-for-byte unchanged, and the signatures are unconstrained (the Send + Sync element requirement is carried by a marker bound that is empty unless rayon is enabled, and is satisfied automatically by f32 / f64).

What is parallelized

Only heap-backed, runtime-sized paths with independent, disjoint output columns — never small fixed-size Matrix ops (thread dispatch would dwarf the work) and never order-sensitive reductions (which would sacrifice reproducibility):

Area	Routines	Axis
optim	finite_difference_jacobian_dyn_par, finite_difference_gradient_dyn_par	columns = independent function evaluations
imageproc separable convolution	convolve2d_separable / _into, gaussian_blur / _into, box_blur, unsharp_mask, laplacian_of_gaussian, canny, Harris / Shi-Tomasi corners, DoG, Gaussian pyramid	bands of output columns (each band's vertical + horizontal pass on a per-worker scratch slab)
imageproc rank/median	rank_filter, percentile_filter, median_filter	output columns (quickselect per pixel)
imageproc geometric/stats	resize_bilinear, local_mean / local_variance / local_stddev, adaptive_threshold	output columns
imageproc morphology	dilate / erode, opening / closing, max/min_filter, gradient, top-hat, black-hat	output columns (horizontal pass via transpose)
fft 2D	DynFft2 / DynRealFft2 (hence fft_convolve2d / fft_correlate2d)	columns of each pass — independent 1D transforms sharing one read-only plan, one DynFftScratch per worker; row pass via a blocked transpose (also banded in parallel). A single 1D FFT is never multithreaded.

The optim parallel routines are separate _par functions (e.g. finite_difference_jacobian_dyn_par) requiring Fn + Sync + Send, rather than a feature-gated change to the sequential finite_difference_jacobian_dyn (which keeps its FnMut signature). This keeps the feature purely additive — enabling rayon anywhere in a dependency tree never alters an existing signature. The imageproc routines take no user closure, so they parallelize in place (the element-type Send + Sync requirement is invisible for f32/f64).

Determinism and gating

• Deterministic results. Each worker writes a disjoint slice of the output, so the result is identical regardless of thread count. (Parallel reductions, where summation order would change the floating-point answer, are intentionally not used.)
• Work-aware gating. Each routine decides whether to fan out based on total work (e.g. nrows · ncols · kernel_size), not raw column count — so a small image, or a cheap operation on a medium image, stays sequential and never pays thread-dispatch overhead. The crossover that matters is the cost of the work, not its shape.
• Core-count aware. The work threshold scales with rayon::current_num_threads(): more cores need more total work to amortize the extra coordination (and more chunks to feed), so a 2-core laptop parallelizes smaller inputs than a 64-core server. The budget constants are tuned on an 8-core reference and normalized to it, so this is a coarse machine-independent guard — not a per-machine optimum (the cost band around the crossover is wide). For the finite-difference Jacobian, where the dominant cost is your own f rather than the machine, the tuning signal is simply which function you call — the opt-in _par variant.

Measured speedups

Apple M3 — 8 cores (4 performance + 4 efficiency), aarch64, f64. Rayon uses its default thread pool (all 8 logical cores), so the heterogeneous P+E cores cap practical scaling below 8×. Parallel vs. the same build with rayon off:

Workload	Size	Speedup
Separable Gaussian blur	512²	~2.6×
Median filter (5×5)	256²	~3.4×
Median filter (3×3)	256²	~3.7×
Morphological dilation (r=3)	1024²	~4.1×
2D FFT (DynFft2, complex)	512²	~2.6×
2D FFT (DynFft2, complex)	1024²	~2.8×
2D real FFT (DynRealFft2, forward / inverse)	1024²	~2.9× / ~2.7×
Finite-difference Jacobian	n=256, expensive f	~4×

Wins scale with total work; below each routine's gate the parallel build matches the sequential one (no regression). See bench/ for the harnesses — each is feature-toggled (--no-default-features flips rayon off) so Criterion can diff the two baselines.