Implementing Planar Laplace Noise in Python

Sample the angle uniformly, sample the radius by inverting the Gamma(2, 1/ε) CDF through the $W_{-1}$ branch of the Lambert function, and add the resulting vector in a projected metric CRS. The whole mechanism is nine lines; the branch choice is the only place it goes wrong, and it goes wrong silently.

Core Calculation # permalink

The planar Laplace density for geo-indistinguishability is

Because it depends only on radial distance, sampling factorises into an angle and a radius. The angle is uniform on $[0, 2\pi)$. The radius has CDF

which is not invertible in elementary functions. Setting $F(r) = p$ and rearranging gives

so with $u = -(1 + \varepsilon r)$ we need $u\, e^{u} = (p-1)/e$, that is $u = W\!\left((p-1)/e\right)$, and therefore

The argument $(p-1)/e$ lies in $[-1/e, 0)$, where the Lambert function has two real branches. $W_0$ returns values in $[-1, 0)$, which makes $r$ negative or near zero. $W_{-1}$ returns values in $(-\infty, -1]$, which is the branch that yields $r \ge 0$. scipy.special.lambertw defaults to $W_0$, so an implementation that omits k=-1 produces displacements clustered at the origin while still looking like a working mechanism.

Worked numeric example # permalink

Take a 200 m privacy radius at a tolerable factor of $e$, so $\varepsilon = \ln e / 200 = 0.005\ \mathrm{m}^{-1}$.

The median displacement is 430 m and the mean is $2/\varepsilon = 400$ m — the distribution is right-skewed, so the median exceeds the mean here only because the mean is pulled by the definition rather than the tail; both are near twice the privacy radius, which is the sanity check to remember.

Python Implementation # permalink

from __future__ import annotations

import geopandas as gpd
import numpy as np
from scipy.special import lambertw

def planar_laplace_perturb(
    gdf: gpd.GeoDataFrame,
    privacy_radius_m: float,
    tolerable_factor: float = np.e,
    seed: int | None = None,
) -> gpd.GeoDataFrame:
    """Perturb point geometries under epsilon-geo-indistinguishability.

    The privacy rate is derived from the radius rather than set directly, so the
    reviewable parameter is a distance. Sampling happens in the local UTM zone:
    epsilon has units of inverse metres and is meaningless against degrees.

    Args:
        gdf: point GeoDataFrame with any CRS (reprojected internally).
        privacy_radius_m: distance within which two locations stay confusable.
        tolerable_factor: the posterior ratio granted at that distance.
        seed: omit in production — a reused seed lets two releases be subtracted.

    Returns:
        A copy in the input CRS, with perturbed geometries and the rate applied.
    """
    if privacy_radius_m <= 0:
        raise ValueError("privacy_radius_m must be positive")
    if gdf.crs is None:
        raise ValueError("input GeoDataFrame has no CRS")

    eps = float(np.log(tolerable_factor) / privacy_radius_m)   # per metre

    original_crs = gdf.crs
    metric = gdf.to_crs(gdf.estimate_utm_crs())
    if not metric.crs.is_projected:                            # guard, not decoration
        raise ValueError("failed to obtain a projected CRS for sampling")

    rng = np.random.default_rng(seed)
    n = len(metric)
    theta = rng.uniform(0.0, 2.0 * np.pi, n)
    p = rng.uniform(0.0, 1.0, n)

    # k=-1 selects the branch that yields non-negative radii. The scipy default
    # (k=0) returns radii collapsed toward zero and silently voids the guarantee.
    w = lambertw((p - 1.0) / np.e, k=-1).real
    r = -(w + 1.0) / eps

    out = metric.copy()
    out["geometry"] = gpd.points_from_xy(
        metric.geometry.x + r * np.cos(theta),
        metric.geometry.y + r * np.sin(theta),
        crs=metric.crs,
    )
    out["displacement_m"] = r
    out.attrs["epsilon_per_m"] = eps
    out.attrs["privacy_radius_m"] = privacy_radius_m
    return out.to_crs(original_crs)

Verification # permalink

Three checks, run together, catch every implementation error this mechanism has.

def verify_planar_laplace(displacement_m: np.ndarray, eps: float) -> dict:
    """Distributional checks against the closed-form properties of Gamma(2, 1/eps)."""
    r = np.asarray(displacement_m, dtype=float)
    expected_mean = 2.0 / eps
    expected_median = 1.678 / eps          # numeric root of F(r) = 0.5
    return {
        "mean_m": float(r.mean()),
        "expected_mean_m": expected_mean,
        "mean_ratio": float(r.mean() / expected_mean),      # target ~1.0
        "median_m": float(np.median(r)),
        "expected_median_m": expected_median,
        "frac_below_10m": float((r < 10.0).mean()),         # target < 0.01
        "all_non_negative": bool((r >= 0).all()),
    }

• mean_ratio near 1.0. A ratio near zero is the $W_0$ branch. A ratio near $10^{-5}$ is degree-space sampling.
• frac_below_10m small. The planar Laplace density vanishes at the origin, so almost nothing should land on top of the true point. A spike near zero is the branch bug.
• all_non_negative true. A negative radius means the branch is wrong and the sign was patched with abs(), which produces a different distribution again.

Pair these with an adversarial check: fit a kernel density to the reported points, combine it with a population prior, and measure the attacker’s expected error. That is the number a re-identification risk assessment actually cares about.

Edge Cases & Adjustments # permalink

• Reports outside the service area. At a 200 m radius, roughly a quarter of reports land beyond 600 m. Resampling until the point is plausible is post-processing and preserves the guarantee, but it biases displacement inward — count the resampled fraction and record it.
• Very small privacy radii. Below about 25 m the mean displacement drops under 50 m, which is inside GPS error, and the mechanism stops being distinguishable from device noise. At that point it is providing a formal guarantee that no adversary was constrained by.
• Repeated reports from a stationary device. The noise is zero-mean and independent, so averaging $n$ reports shrinks the attacker’s error as $1/\sqrt{n}$. Memoize the report against the true location rather than re-drawing, as under local differential privacy for mobile clients.
• Seeding. The seed argument exists for tests. A seed that reaches production lets an adversary subtract two releases and recover the true movement exactly — the failure documented in Laplace and Gaussian noise for coordinate data.
• Multi-zone extents. estimate_utm_crs picks one zone from the whole frame. For a national dataset, group by zone and perturb each group in its own projection, or the effective radius drifts at the edges.

FAQ # permalink

Why not sample the radius from numpy.random.gamma(2, 1/eps) directly?

You can, and it is equivalent — Gamma(2, 1/ε) is exactly the radial distribution. The Lambert form is shown here because it is what the literature states and what most implementations copy, which makes the branch bug worth documenting. If you use the Gamma sampler, keep the same verification.

Does the mechanism need a truncation radius?

Not for the guarantee. Truncating is post-processing, so it is safe, but it changes the distribution and therefore the utility figures. If your application cannot tolerate a 1 km outlier, the honest conclusion is that it cannot tolerate a 200 m privacy radius.

Can I apply this to a trajectory rather than isolated points?

Not point by point. Independent perturbation of consecutive points is defeated by map matching and by averaging along the path — the mechanisms described in defending against map-matching attacks. Trajectories need sequence-aware controls.

What ε do I put in the release record?

Both numbers: the rate in inverse metres and the privacy radius it corresponds to. Recording only the rate leaves a reviewer to work out what it means, and recording only the radius loses the parameter the mechanism actually used.

Related # permalink

• Geo-Indistinguishability & the Planar Laplace Mechanism — the definition and where the parameter comes from
• Converting Epsilon per Metre to a Privacy Radius — the unit conversion, with a review-ready table
• Applying Laplace Noise to Latitude/Longitude Pairs — the central-model counterpart and its CRS traps
• Local Differential Privacy for Mobile Clients — the discrete on-device alternative

← Back to Geo-Indistinguishability & the Planar Laplace Mechanism