One LOESS to Rule Them All: Rebuilding a Classic Smoother in Rust

LOESS (locally estimated scatterplot smoothing) is one of those algorithms almost everyone has used without thinking about it. It draws the trend line through a noisy scatterplot, it normalizes microarray data, and it is the default smoother behind a lot of plots in R. It is also considered finished. The reference implementations descend from Fortran written decades ago (Cleveland's LOWESS in 1979, LOESS in 1988), and for most people "LOESS" means "whatever stats::loess does".

I wanted to see what happens if you treat a finished algorithm as a starting point instead. The result is loess-project: a Rust implementation with a modular core, a much larger feature set, and bindings for nine ecosystems.

What LOESS does, in one paragraph

For every point where you want a smoothed value, LOESS takes the nearest fraction of the data, weights those neighbours by distance, fits a small polynomial by weighted least squares, and reads off the fitted value. Repeat for every point and you get a smooth curve that follows the data without assuming a global shape. Because it is local and non-parametric, it is useful for trend estimation in noisy time series, normalization in genomics, and exploratory analysis of scatterplots.

Step 1: match the reference exactly

Before adding anything, the implementation has to agree with the thing people already trust. The project validates against R's loess, and the README reports a maximum difference below 1e-12 across the test scenarios. That gave me a safe baseline: every later feature or optimization is checked against it.

Step 2: build on it in a modular way

Once the baseline matched, the design goal was to make each decision in the algorithm a swappable part instead of a hard-coded constant. R's loess is fixed in several places; the Rust core exposes them:

	This project	R's loess
Polynomial degree	0 to 4	1 or 2
Kernel	7 options	tricube only
Robustness weighting	3 options	bisquare only
Scale estimation	3 options (including MAD)	median absolute residual only
Distance metric	6 options	normalized only
Boundary padding	4 options	none

It also adds things the original never had: streaming and online modes for data that arrives continuously, confidence and prediction intervals, diagnostics (RMSE, R², AIC), cross-validation, and no_std support.

The crates reflect this layering:

• loess-rs is the core: single-threaded, no_std-capable, and dependency-light.
• fastLoess wraps it with parallel execution and ndarray integration.

Two choices that matter statistically: robustness weights use the median absolute deviation (MAD) for scale, which stays valid with up to 50% outliers, and boundary padding policies reduce the asymmetric-neighbourhood bias that plain LOESS has at the edges of the data.

Inside the core: six layers

Inside loess-rs, the modules form six numbered layers, from the bottom up: primitives (data structures), math (distances, kernels, linear algebra), algorithms (the core LOESS fitting, including specialized 1D to 3D paths), evaluation (cross-validation, intervals, diagnostics), engine (orchestration), and adapters (the batch, streaming and online modes users pick from). A public API sits on top, then fastLoess, then the language bindings.

Two rules hold throughout: a layer depends only on the layers below it, never above, and modules within the same layer don't depend on each other. For an algorithm with this many moving parts, that has practical benefits:

• Isolated testing. A kernel or distance function can be tested, and its SIMD path compared against the scalar one, without building a whole smoother.
• Contained change. A new kernel touches the math layer, a new robustness method the algorithms layer, and a new way of feeding data in is just an adapter.
• No cycles, easier reasoning. Data and control flow only one way, so you can reason about each layer on its own.
• Reuse. Batch, streaming and online share the same engine and algorithms instead of reimplementing them.

Step 3: make it fast

Rust gives two straightforward levers here. The first is parallelism across the points being smoothed (via the fastLoess wrapper). The second is SIMD: the specialized 1D, 2D and 3D regression accumulators, the distance kernels and the batch residual operations use the wide crate rather than relying on auto-vectorization.

On the project's benchmark harness against R's stats::loess, the README reports roughly 47x in serial mode and 236x in parallel mode. Those numbers vary by workload, and parallel is not always faster: for datasets of about a thousand points or fewer, serial is the better choice.

Mean runtime of R's stats::loess versus the serial and parallel R binding (rfastloess) across scalability, financial, scientific, genomic and pathological workloads. CPU only, and results depend on hardware and thread count. Reproduce with make compare in benchmarks/.

Step 4: bring the same engine to everyone else

The part I find most interesting is not the algorithm but the distribution. A statistician in R, a data scientist in Python and a front-end developer in the browser should all get the same numerical engine in a form that feels native to them. Rust's binding ecosystem makes that practical:

Community	Binding approach
R	extendr
Python	PyO3 with maturin
Julia	C ABI called with ccall
Node.js	napi-rs
WebAssembly	wasm-bindgen
C++	C ABI generated by cbindgen, plus CMake packaging
Go	C ABI from cbindgen through cgo
Java	JNI

Every binding wraps the same core, so a fix or a new smoothing option lands everywhere at once, and the Python and R bindings are checked numerically against R's loess in CI.

Here is the same robust fit in three languages. Six points, one obvious outlier (50.0), identical parameters, and a check for which points the robustness iterations downweighted.

Rust

use fastLoess::prelude::*;

fn main() -> Result<(), LoessError> {
    let x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
    let y = vec![2.0, 4.0, 6.0, 50.0, 10.0, 12.0];

    let model = Loess::new()
        .fraction(0.7)
        .iterations(5)
        .robustness_method("bisquare")
        .outputs(["weights"])
        .build()?;

    let result = model.fit(&x, &y)?;

    if let Some(weights) = &result.robustness_weights {
        for (i, w) in weights.iter().enumerate() {
            if *w < 0.5 {
                println!("Point {} is likely an outlier (weight: {:.3})", i, w);
            }
        }
    }
    Ok(())
}

Python

import fastloess as fl
import numpy as np

x = np.linspace(1, 6, 6)
y = np.array([2.0, 4.0, 6.0, 50.0, 10.0, 12.0])

model = fl.Loess(
    fraction=0.7,
    iterations=5,
    robustness_method="bisquare",
    outputs=["weights"],
)
result = model.fit(x, y)

for i, w in enumerate(result.robustness_weights):
    if w < 0.5:
        print(f"Point {i} is likely an outlier (weight: {w:.3f})")

R

library(rfastloess)

x <- 1:6
y <- c(2.0, 4.0, 6.0, 50.0, 10.0, 12.0)

model <- Loess(
    fraction = 0.7,
    iterations = 5,
    robustness_method = "bisquare",
    outputs = "weights"
)
result <- fit(model, x, y)

for (i in seq_along(result$robustness_weights)) {
    if (result$robustness_weights[i] < 0.5) {
        cat(sprintf("Point %d is likely an outlier (weight: %.3f)\n",
                    i, result$robustness_weights[i]))
    }
}

Each one flags the 50.0 value (index 3 in Rust and Python, 4 in R because R counts from 1). The option names are the same in every language, so code and documentation carry over between communities with very little translation.

Step 5: keep it all in one repo

Nine bindings could easily become nine repositories drifting apart. Instead everything lives in one monorepo, organized by role:

• crates/ holds the Rust libraries.
• bindings/ holds one directory per language, each with its own build, docs and tests.
• validation/ holds the cross-implementation checks.
• benchmarks/ and dev/ hold the comparison harness and release tooling.

Each binding can be built and tested on its own, but they all consume the same core. A root Makefile delegates to per-component Makefiles, so make python-dev, make julia-dev or make go-dev runs the full format, lint, build and test cycle for just that binding, and make all-dev runs everything.

What it takes to keep nine bindings in step

Automation per language and per registry. Every binding has its own CI workflow, and every registry gets its own release workflow (crates.io, PyPI, npm, conda-forge, Maven Central, the Julia registry, the Go module proxy, and a GitHub release for the C++ artifacts). A tagged release fans out from one place, but each package ecosystem still has to be driven in its own way.

One command for versions. A version bump has to land in every crate and manifest at once, so a script (dev/bump_version.py) does it in one step. Without that, version drift between bindings is the first thing that breaks. A weekly workflow also checks the hardcoded tool and library version pins against their latest releases and opens an issue when any are outdated. The bump script has to know ecosystem quirks too: a Go module tagged v2 or later must carry a /vN suffix in its module path, or the Go toolchain and pkg.go.dev silently ignore the tag.

Docs are tested, not just written. Documentation lives beside each crate and binding, in whatever tool that community expects (rustdoc, Sphinx, pkgdown, Documenter.jl, Doxygen, and others). To keep it honest, dev/verify_snippets.py extracts the fenced code blocks from the docs and runs them, and dev/check_links.py validates the cross-references. make docs-test runs the snippet check across all languages, which means an API change that breaks an example fails CI instead of surviving as a stale doc.

Packaging has its own surprises. A few that cost real time:

• With the napi-rs v3 CLI, a plain napi build emits only a generic .node file. The platform-suffixed binary that the generated loader actually requires appears only if you pass --platform.
• A Go module that wraps a Rust library through cgo has to ship its native artifacts inside the module release, which needed its own packaging script (dev/package_go_module.sh).
• Registries have their own rules. A Spack review wanted the recipe's url to stay on a stable earlier release tag rather than the newest one, which meant fixing the release workflow that updates the recipe. A vcpkg review pushed back on a port built from prebuilt binaries because of link compatibility with other ports.

What I would tell someone starting this. Automate the boring parts (versions, link checks, snippet checks) before adding the next language, because each new binding multiplies the amount of drift you have to catch. Treat each registry as a separate product with its own reviewers and conventions. And keep the numerical validation against the reference implementation in CI (here make validate runs the Python and R checks against R's loess) so that "same engine everywhere" stays a tested claim rather than an assumption.

Takeaways

• A "finished" algorithm still has room: matching the reference first gives you permission to extend it.
• Modular design turns hard-coded constants into options, which is where most of the new functionality came from.
• I think this approach could extend well beyond LOESS. Much of scientific computing rests on mature algorithms with decades-old cores that are wrapped separately for each language. A single memory-safe, fast Rust core with thin bindings for each community could modernize that kind of package without asking anyone to leave the language they already use.
• One well-organized monorepo can serve many language communities without the bindings drifting apart.

The code is on GitHub at thisisamirv/loess-project under MIT or Apache-2.0, and contributions are welcome.

Which other "finished" algorithm do you think deserves this treatment: a trusted classic whose reference implementation is decades old, and that every language community wraps in its own way? I would like to hear which ones you would pick, and what would make you trust a rewrite enough to switch.