Benchmarking pandas vs Polars for SPC Pipelines

This is the measured way to decide whether Polars is worth adopting for the limit-calculation stage of your pipeline, built as a fair, reproducible benchmark rather than a vibe. It extends the vectorized control limit calculation engines guide with a concrete head-to-head: the same X-bar and R limit calculation over the same synthetic dataset, timed on wall clock and peak memory, with an assertion that both engines produce identical limits within a tight tolerance. A benchmark that skips the parity check is worthless — a faster engine that computes a different number is not faster, it is wrong.

Prerequisites

  • Python 3.9+ with pandas, polars, and numpy installed (pip install pandas polars numpy).
  • The vectorized pandas engine from the parent guide, or an equivalent function that returns per-characteristic X-bar and R limits from a tidy frame.
  • A machine whose memory and CPU roughly match your production target; benchmark numbers do not transfer across very different hardware.
  • Familiarity with fixed subgroup size and the $A_2$/$D_3$/$D_4$ constants — see how to calculate control limits for X-bar R charts in Python.
  • Nothing else heavy running: close background jobs so the timer measures the engine, not contention.

Step 1 — Generate a fixed synthetic dataset

A fair benchmark feeds both engines the exact same rows. Generate a tidy long-format frame once, with a controllable number of characteristics and subgroups, so you can scale N and re-run without changing the shape of the work.

from __future__ import annotations

import numpy as np
import pandas as pd


def make_dataset(n_characteristics: int, n_subgroups: int, n: int = 5,
                 seed: int = 7) -> pd.DataFrame:
    """Build a tidy SPC dataset with a fixed subgroup size n.

    Returns one row per measurement with columns
    ``characteristic``, ``subgroup`` and ``value`` (float64).
    """
    if n < 2:
        raise ValueError("Subgroup size n must be at least 2.")
    rng = np.random.default_rng(seed)
    rows_per_char = n_subgroups * n
    frames = []
    for c in range(n_characteristics):
        center = 20.0 + c  # distinct center per characteristic
        values = rng.normal(center, 0.05, size=rows_per_char).astype("float64")
        frames.append(pd.DataFrame({
            "characteristic": f"char_{c:04d}",
            "subgroup": np.repeat(np.arange(n_subgroups), n),
            "value": values,
        }))
    return pd.concat(frames, ignore_index=True)


DATA = make_dataset(n_characteristics=500, n_subgroups=200, n=5)
print(f"rows: {len(DATA):,}")  # 500 * 200 * 5 = 500,000

Half a million rows is enough to separate the engines; scale n_characteristics up to stress memory and group cardinality, which is where Polars tends to pull ahead.

Step 2 — Define the pandas limit calculation

Keep the pandas path to a single grouped reduction. This mirrors the production engine: group by characteristic and subgroup, reduce to means and ranges, collapse to centerlines, and apply the constants.

_A2 = {2: 1.880, 3: 1.023, 4: 0.729, 5: 0.577, 6: 0.483,
       7: 0.419, 8: 0.373, 9: 0.337, 10: 0.308}
_D3 = {2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0,
       7: 0.076, 8: 0.136, 9: 0.184, 10: 0.223}
_D4 = {2: 3.267, 3: 2.574, 4: 2.282, 5: 2.114, 6: 2.004,
       7: 1.924, 8: 1.864, 9: 1.816, 10: 1.777}


def limits_pandas(df: pd.DataFrame, n: int = 5) -> pd.DataFrame:
    """Per-characteristic X-bar and R limits, computed with pandas."""
    sub = df.groupby(["characteristic", "subgroup"], sort=False)["value"].agg(
        mean="mean", _max="max", _min="min"
    )
    sub["range"] = sub["_max"] - sub["_min"]
    per = sub.groupby("characteristic").agg(
        x_double_bar=("mean", "mean"), r_bar=("range", "mean")
    ).reset_index()
    a2, d3, d4 = _A2[n], _D3[n], _D4[n]
    per["xbar_ucl"] = per["x_double_bar"] + a2 * per["r_bar"]
    per["xbar_lcl"] = per["x_double_bar"] - a2 * per["r_bar"]
    per["r_ucl"] = d4 * per["r_bar"]
    per["r_lcl"] = d3 * per["r_bar"]
    return per.sort_values("characteristic").reset_index(drop=True)

Step 3 — Define the equivalent Polars calculation

The Polars path expresses the same logic with group_by and agg. Use the lazy API so the query optimizer can fuse the two aggregations and prune columns; call .collect() once at the end.

import polars as pl


def limits_polars(df_pl: pl.DataFrame, n: int = 5) -> pl.DataFrame:
    """Per-characteristic X-bar and R limits, computed with Polars (lazy)."""
    a2, d3, d4 = _A2[n], _D3[n], _D4[n]
    lf = df_pl.lazy()
    sub = lf.group_by(["characteristic", "subgroup"]).agg(
        pl.col("value").mean().alias("mean"),
        (pl.col("value").max() - pl.col("value").min()).alias("range"),
    )
    per = sub.group_by("characteristic").agg(
        pl.col("mean").mean().alias("x_double_bar"),
        pl.col("range").mean().alias("r_bar"),
    ).with_columns(
        (pl.col("x_double_bar") + a2 * pl.col("r_bar")).alias("xbar_ucl"),
        (pl.col("x_double_bar") - a2 * pl.col("r_bar")).alias("xbar_lcl"),
        (d4 * pl.col("r_bar")).alias("r_ucl"),
        (d3 * pl.col("r_bar")).alias("r_lcl"),
    ).sort("characteristic")
    return per.collect()

The lazy engine wins most clearly when the query has structure to optimize: here it fuses the reshape and the two aggregations, and it keeps only the columns each stage needs, so peak memory stays below what the eager pandas path holds.

Step 4 — Measure wall time and peak memory

Time both engines under identical conditions. Use time.perf_counter for wall time and tracemalloc for peak Python allocation, warm up once to exclude import and first-touch costs, and take the best of several repeats to reduce noise.

import time
import tracemalloc
from typing import Callable


def benchmark(fn: Callable[[], object], repeats: int = 5) -> dict[str, float]:
    """Return best wall time (seconds) and peak memory (MiB) over repeats."""
    fn()  # warm up: exclude import / first-touch allocation
    best_time = float("inf")
    peak_mem = 0.0
    for _ in range(repeats):
        tracemalloc.start()
        t0 = time.perf_counter()
        fn()
        elapsed = time.perf_counter() - t0
        _, peak = tracemalloc.get_traced_memory()
        tracemalloc.stop()
        best_time = min(best_time, elapsed)
        peak_mem = max(peak_mem, peak / 1024 / 1024)
    return {"seconds": best_time, "peak_mib": peak_mem}


DATA_PL = pl.from_pandas(DATA)
pd_result = benchmark(lambda: limits_pandas(DATA))
pl_result = benchmark(lambda: limits_polars(DATA_PL))

rows = len(DATA)
print(f"pandas : {pd_result['seconds']:.3f}s  "
      f"{rows / pd_result['seconds']:,.0f} rows/s  {pd_result['peak_mib']:.1f} MiB")
print(f"polars : {pl_result['seconds']:.3f}s  "
      f"{rows / pl_result['seconds']:,.0f} rows/s  {pl_result['peak_mib']:.1f} MiB")

Report throughput as rows per second alongside raw seconds, so the number stays meaningful when you change N. Peak memory matters as much as speed for a pipeline that runs many characteristics concurrently — an engine that is 20% faster but doubles peak memory may not fit the box.

Step 5 — Keep the two paths numerically identical

Exclude the pl.from_pandas conversion from the timed region — it is a one-time cost your production code pays once, not per calculation. Convert once, before benchmarking, so the timer measures the aggregation and not the marshaling.

def timed_conversion_note() -> None:
    """Conversion is a fixed setup cost, kept out of the timed loop."""
    # DATA_PL was built once in Step 4, before either benchmark() call.
    # In production, ingest directly into the engine's native frame to skip it.
    ...

If your ingestion already lands data in Arrow or Parquet, read it straight into Polars and skip the pandas hop entirely; forcing a pandas -> polars conversion on every batch would erase the advantage the benchmark reveals.

Verification

The benchmark only means something if both engines agree on the limits. Align both results by characteristic and assert every column matches within a tight float64 tolerance.

def assert_parity(a: pd.DataFrame, b_pl: pl.DataFrame, atol: float = 1e-9) -> None:
    """Assert pandas and Polars limits agree within tolerance."""
    b = b_pl.to_pandas().sort_values("characteristic").reset_index(drop=True)
    a = a.sort_values("characteristic").reset_index(drop=True)
    assert list(a["characteristic"]) == list(b["characteristic"]), "key mismatch"
    for col in ["x_double_bar", "r_bar", "xbar_ucl", "xbar_lcl", "r_ucl", "r_lcl"]:
        max_diff = float(np.abs(a[col].to_numpy() - b[col].to_numpy()).max())
        assert max_diff <= atol, f"{col} differs by {max_diff:.2e} (> {atol:.0e})"
    print("parity OK: pandas and Polars agree within tolerance")


assert_parity(limits_pandas(DATA), limits_polars(DATA_PL))

Expect the maximum difference to be at the level of floating-point rounding (below 1e-12 in practice), because both engines apply the same arithmetic to the same float64 inputs. If parity fails, stop and diagnose before trusting any timing — a difference here means one path is computing the wrong statistic, not that one is faster. Both engines feed the same downstream out-of-control rule detection, so a parity gap would silently change which points alarm.

Root-Cause Table

Symptom Cause Fix
Polars looks slower than pandas Timing includes pl.from_pandas conversion each run, or the eager API is used Convert once outside the timed loop; use the lazy API and a single .collect()
Parity assertion fails by a large margin A different statistic in one path — for example range as std, or a wrong constant Align the exact formula; range is max minus min, and the same $A_2$/$D_3$/$D_4$ table feeds both
Parity fails by a tiny margin near the tolerance Legitimate float ordering differences in summation Loosen atol to a still-tight value such as 1e-9; do not paper over large gaps
Neither engine scales as N grows A hidden per-group Python loop, or object-dtype keys forcing slow hashing Keep the reduction fully vectorized; cast key columns to categorical or string dtype
Memory numbers are noisy or wrong tracemalloc misses native allocations outside Python Corroborate with an OS-level resident-memory sample; treat tracemalloc peak as a lower bound
Results flip run to run Background load or too few repeats Close other jobs, warm up once, take the best of several repeats

This how-to belongs to vectorized control limit calculation engines, within automated control chart generation and calculation.