LOESS Project

The fastest, most robust, and most feature-complete language-agnostic LOESS (Locally Estimated Scatterplot Smoothing) implementation for Rust, Python, R, Julia, JavaScript, C++, and WebAssembly.

[!IMPORTANT]

The loess-project contains a complete ecosystem for LOESS smoothing:

loess-rs - Core single-threaded Rust implementation with no_std support

fastLoess - Parallel Rust wrapper with ndarray integration

R bindings - extendr-based R binding

Python bindings - PyO3-based Python binding

Julia bindings - Native Julia binding with C FFI

JavaScript bindings - Node.js binding

WebAssembly bindings - WASM binding

C++ bindings - Native C++ binding with CMake integration

Installation

[!NOTE]

Currently available for R, Python, Rust, Julia, Node.js, WebAssembly, and C++. See the Installation Guide for detailed installation instructions.

Documentation

[!NOTE]

📚 View the full documentation

LOESS vs. LOWESS

Feature	LOESS (This Crate)	LOWESS
Polynomial Degree	Linear, Quadratic, Cubic, Quartic	Linear (Degree 1)
Dimensions	Multivariate (n-D support)	Univariate (1-D only)
Flexibility	High (Distance metrics)	Standard
Complexity	Higher (Matrix inversion)	Lower (Weighted average/slope)

[!TIP] Note: For a LOWESS implementation, use lowess-project.

Why this package?

Speed

The loess project beats the competition in terms of speed, whether in single-threaded or multi-threaded parallel execution. It is typically 5–20x faster than R's loess in serial mode, and up to 200x faster on large datasets with parallel execution.

For more details on the performance comparison, see the Benchmarks page.

Robustness

This implementation is more robust than R's loess due to two key design choices:

MAD-Based Scale Estimation:

For robustness weight calculations, this crate uses Median Absolute Deviation (MAD) for scale estimation:

s = median(|r_i - median(r)|)

In contrast, R's loess uses the median of absolute residuals (MAR):

s = median(|r_i|)

MAD is a breakdown-point-optimal estimator—it remains valid even when up to 50% of data are outliers.
The median-centering step removes asymmetric bias from residual distributions.
MAD provides consistent outlier detection regardless of whether residuals are centered around zero.

Boundary Padding:

This crate applies a range of different boundary policies at dataset edges:

Extend: Repeats edge values to maintain local neighborhood size.
Reflect: Mirrors data symmetrically around boundaries.
Zero: Pads with zeros (useful for signal processing).
NoBoundary: Original Cleveland behavior

R's loess does not apply boundary padding, which can lead to:

Biased estimates near boundaries due to asymmetric local neighborhoods.
Increased variance at the edges of the smoothed curve.

Features

A variety of features, supporting a range of use cases:

Feature	This package	R (stats)
Polynomial Degree	5 (0–4)	2 (1 or 2)
Kernel	7 options	only Tricube
Robustness Weighting	3 options	only Bisquare
Scale Estimation	3 options	only MAR
Distance Metric	6 options	normalized only
Boundary Padding	4 options	no padding
Zero Weight Fallback	3 options	no
Auto Convergence	yes	no
Online Mode	yes	no
Streaming Mode	yes	no
Confidence Intervals	yes	no
Prediction Intervals	yes	no
Diagnostics (RMSE, R², AIC)	yes	no
Cross-Validation	2 options	no
Parallel Execution	yes	no
`no-std` Support	yes	no

Validation

All implementations are numerical twins of R's loess:

Aspect	Status	Details
Accuracy	✅ EXACT MATCH	Max diff < 1e-12 across all scenarios
Consistency	✅ PERFECT	Multiple scenarios pass with strict tolerance
Robustness	✅ VERIFIED	Robust smoothing matches R exactly

API Reference

library(rfastloess)

model <- Loess(
    fraction = 0.67,
    iterations = 3L,
    weight_function = "tricube",
    robustness_method = "bisquare",
    zero_weight_fallback = "use_local_mean",
    boundary_policy = "extend",
    scaling_method = "mad",
    confidence_intervals = NULL,
    prediction_intervals = NULL,
    return_diagnostics = FALSE,
    return_residuals = FALSE,
    return_robustness_weights = FALSE,
    cv_fractions = NULL,
    cv_method = "kfold",
    cv_k = 5L,
    auto_converge = NULL,
    parallel = TRUE
)
custom_weights <- rep(1, length(x))
result <- model$fit(x, y, custom_weights = custom_weights)

# Result structure:
result$x,
result$y,
result$standard_errors,
result$confidence_lower,
result$confidence_upper,
result$prediction_lower,
result$prediction_upper,
result$residuals,
result$robustness_weights,
result$diagnostics,
result$iterations_used,
result$fraction_used,
result$cv_scores

Python:

from fastloess import Loess

model = Loess(
    fraction=0.67,
    iterations=3,
    weight_function="tricube",
    robustness_method="bisquare",
    zero_weight_fallback="use_local_mean",
    boundary_policy="extend",
    scaling_method="mad",
    confidence_intervals=None,
    prediction_intervals=None,
    return_diagnostics=False,
    return_residuals=False,
    return_robustness_weights=False,
    cv_fractions=None,
    cv_method="kfold",
    cv_k=5,
    auto_converge=None,
    parallel=True
)
custom_weights = [1.0] * len(x)
result = model.fit(x, y, custom_weights=custom_weights)

# Result structure:
result.x,
result.y,
result.standard_errors,
result.confidence_lower,
result.confidence_upper,
result.prediction_lower,
result.prediction_upper,
result.residuals,
result.robustness_weights,
result.diagnostics,
result.iterations_used,
result.fraction_used,
result.cv_scores

Rust:

use loess_rs::prelude::*;

let model = Loess::new()
    .fraction(0.67)
    .iterations(3)
    .weight_function("tricube")
    .robustness_method("bisquare")
    .zero_weight_fallback("use_local_mean")
    .boundary_policy("extend")
    .scaling_method("mad")
    .confidence_intervals(0.95)
    .prediction_intervals(0.95)
    .return_diagnostics()
    .return_residuals()
    .return_robustness_weights()
    .cv_fractions(vec![0.3, 0.5, 0.7])
    .cv_method("kfold")
    .cv_k(5)
    .cv_seed(123)
    .auto_converge(1e-4)
    .custom_weights(vec![1.0; x.len()])
    .parallel(true)             // fastLoess only
    .build()?;

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

// Result structure:
pub struct LoessResult<T> {
    pub x: Vec<T>,                           // Sorted x values
    pub y: Vec<T>,                           // Smoothed y values
    pub standard_errors: Option<Vec<T>>,
    pub confidence_lower: Option<Vec<T>>,
    pub confidence_upper: Option<Vec<T>>,
    pub prediction_lower: Option<Vec<T>>,
    pub prediction_upper: Option<Vec<T>>,
    pub residuals: Option<Vec<T>>,
    pub robustness_weights: Option<Vec<T>>,
    pub diagnostics: Option<Diagnostics<T>>,
    pub iterations_used: Option<usize>,
    pub fraction_used: T,
    pub cv_scores: Option<Vec<T>>,
}

Julia:

using FastLOESS

model = Loess(;
    fraction=0.67,
    iterations=3,
    weight_function="tricube",
    robustness_method="bisquare",
    zero_weight_fallback="use_local_mean",
    boundary_policy="extend",
    scaling_method="mad",
    confidence_intervals=NaN,
    prediction_intervals=NaN,
    return_diagnostics=false,
    return_residuals=false,
    return_robustness_weights=false,
    cv_fractions=Float64[], # e.g. [0.3, 0.5]
    cv_method="kfold",
    cv_k=5,
    auto_converge=NaN,
    parallel=true
)
custom_weights = ones(length(x))
result = fit(model, x, y; custom_weights=custom_weights)

# Result structure:
result.x,
result.y,
result.standard_errors,
result.confidence_lower,
result.confidence_upper,
result.prediction_lower,
result.prediction_upper,
result.residuals,
result.robustness_weights,
result.diagnostics,
result.iterations_used,
result.fraction_used,
result.cv_scores

Node.js:

import { Loess } from "fastloess"

const model = new Loess({
    fraction: 0.67,
    iterations: 3,
    weight_function: "tricube",
    robustness_method: "bisquare",
    zero_weight_fallback: "use_local_mean",
    boundary_policy: "extend",
    scaling_method: "mad",
    confidence_intervals: 0.95,
    prediction_intervals: 0.95,
    return_diagnostics: true,
    return_residuals: true,
    return_robustness_weights: true,
    cv_fractions: [0.3, 0.5, 0.7],
    cv_method: "kfold",
    cv_k: 5,
    auto_converge: 1e-4,
    parallel: true
})
const custom_weights = Array(x.length).fill(1.0)
const result = model.fit(x, y, custom_weights)

// Result structure:
result.x,
result.y,
result.standard_errors,
result.confidence_lower,
result.confidence_upper,
result.prediction_lower,
result.prediction_upper,
result.residuals,
result.robustness_weights,
result.diagnostics,
result.iterations_used,
result.fraction_used,
result.cv_scores

WebAssembly:

import { Loess } from "fastloess-wasm"

const model = new Loess({
    fraction: 0.67,
    iterations: 3,
    weight_function: "tricube",
    robustness_method: "bisquare",
    zero_weight_fallback: "use_local_mean",
    boundary_policy: "extend",
    scaling_method: "mad",
    confidence_intervals: 0.95,
    prediction_intervals: 0.95,
    return_diagnostics: true,
    return_residuals: true,
    return_robustness_weights: true,
    cv_fractions: [0.3, 0.5, 0.7],
    cv_method: "kfold",
    cv_k: 5,
    auto_converge: 1e-4,
    parallel: true
})
const custom_weights = new Float64Array(x.length).fill(1)
const result = model.fit(x, y, custom_weights)

// Result structure:
result.x,
result.y,
result.standard_errors,
result.confidence_lower,
result.confidence_upper,
result.prediction_lower,
result.prediction_upper,
result.residuals,
result.robustness_weights,
result.diagnostics,
result.iterations_used,
result.fraction_used,
result.cv_scores

C++:

fastloess::LoessOptions options;
options.fraction = 0.67;
options.iterations = 3;
options.weight_function = "tricube";
options.robustness_method = "bisquare";
options.zero_weight_fallback = "use_local_mean";
options.boundary_policy = "extend";
options.scaling_method = "mad";
options.confidence_intervals = 0.95;
options.prediction_intervals = 0.95;
options.return_diagnostics = true;
options.return_residuals = true;
options.return_robustness_weights = true;
options.cv_fractions = {0.3, 0.5, 0.7};
options.cv_method = "kfold";
options.cv_k = 5;
options.auto_converge = 1e-4;
options.parallel = true;

fastloess::Loess model(options);
std::vector<double> custom_weights(x.size(), 1.0);
const auto result = model.fit(x, y, custom_weights).value();

// Result structure:
result.x_vector(),
result.y_vector(),
result.standard_errors(),
result.confidence_lower(),
result.confidence_upper(),
result.prediction_lower(),
result.prediction_upper(),
result.residuals(),
result.robustness_weights(),
result.diagnostics(),
result.iterations_used(),
result.fraction_used(),
result.cv_scores()

Contributing

Contributions are welcome! Please see the Contributing Guide for more information.

Changelog

See the Changelog for a history of changes.

License

Licensed under MIT or Apache-2.0.

Citation

If you use this software in your research, please cite it using the CITATION.cff file or the BibTeX entry below:

@software{loess_project,
  author = {Valizadeh, Amir},
  title = {LOESS Project: High-Performance Locally Estimated Scatterplot Smoothing},
  year = {2026},
  url = {https://github.com/thisisamirv/loess-project},
  license = {MIT OR Apache-2.0}
}

loess-rs 0.9.0