Planar Laplace Noise for Location Privacy

The planar Laplace distribution is the two-dimensional noise that a report-a-fuzzed-point mechanism draws from, and getting it right is almost entirely a sampling problem: you cannot add independent Laplace noise to latitude and longitude and call it done, because that produces a diamond-shaped, direction-leaking cloud rather than the isotropic one the privacy proof assumes. This page is the sampler deep-dive inside the parent guide on spatial noise mechanisms and the broader Differential Privacy for Geospatial Data architecture. It is deliberately narrow: the privacy semantics — why $\varepsilon$ -per-metre indistinguishability protects a location — live on the geo-indistinguishability pages, while here we cover only the density $D(\mathbf{x}) = \frac{\varepsilon^2}{2\pi} e^{-\varepsilon \lVert \mathbf{x} \rVert}$ , how to sample from it exactly through its polar decomposition, and how to project the resulting metric offset back to WGS84 degrees without introducing a fingerprint.

The trick that makes planar Laplace tractable is that in polar coordinates the two components are independent. The angle $\theta$ is uniform on $[0, 2\pi)$ — the mechanism is isotropic, so it leaks no bearing — and the radius $r$ has marginal density proportional to $r\,e^{-\varepsilon r}$ , a Gamma$(2, 1/\varepsilon)$ law whose mean is $2/\varepsilon$ metres. Sampling the radius is the only hard part, and it requires the Lambert $W$ function.

Parameter Configuration and Calibration

Planar Laplace has exactly one privacy knob, $\varepsilon$ , plus a handful of implementation constants that are correctness-critical rather than tunable. Because $\varepsilon$ here carries units of inverse metres, it behaves differently from the scalar $\varepsilon$ of the Laplace mechanism used for Gaussian and Laplace spatial aggregates: two locations $d$ metres apart are $\varepsilon d$ -indistinguishable, so the operationally meaningful quantity is not $\varepsilon$ itself but the privacy radius $\ell / \varepsilon$ at a chosen indistinguishability level $\ell$ .

epsilon — the per-metre privacy parameter. Pick it indirectly: decide a level $\ell$ (for example $\ell = \ln 2$ , “at most a 2× likelihood ratio”) and a radius $R$ at which that level should hold, then set $\varepsilon = \ell / R$ . With $\ell = \ln 2$ and $R = 200$ m, $\varepsilon \approx 3.47 \times 10^{-3}\,\text{m}^{-1}$ , giving a mean reported displacement of $2/\varepsilon \approx 577$ m. Smaller $\varepsilon$ means more noise and stronger privacy; the mean radius $2/\varepsilon$ is the number to sanity-check against your map.
Quantile clamp — numerical, not statistical. The inverse CDF is evaluated at a uniform quantile $p \in (0,1)$ . As $p \to 1$ the radius diverges (the density has an unbounded tail) and as $p \to 0$ the argument lands exactly on the Lambert $W$ branch point $-1/e$ where precision collapses. Clamp $p$ into $[10^{-12}, 1 - 10^{-12}]$ ; this bounds the tail without meaningfully perturbing the distribution.
Reference latitude for the metre→degree projection. The metric offset is generated in metres and must be converted to degrees at the local scale. One degree of latitude is ~111,320 m everywhere, but one degree of longitude shrinks by $\cos(\text{lat})$ . Use the true point’s latitude as the reference; skipping the cosine term inflates the eastward error by up to $2\times$ at 60° and is a silent correctness bug at high latitude.

The radius CDF you invert is $C(r) = 1 - (1 + \varepsilon r)\,e^{-\varepsilon r}$ . Setting $C(r) = p$ and solving gives the closed form used throughout:

r = -\frac{1}{\varepsilon}\left( W_{-1}\!\left( \frac{p - 1}{e} \right) + 1 \right)

The $-1$ branch of $W$ is mandatory. Because $1 + \varepsilon r \ge 1$ , the substitution forces $W(\cdot) \le -1$ , which is exactly the range of $W_{-1}$ ; the principal branch $W_0$ returns a value in $[-1, 0]$ and would silently yield the wrong (negative) root.

Design goal	Control	Concrete setting
$\ln 2$ privacy within 200 m	$\varepsilon = \ln 2 / 200$	mean radius $\approx 577$ m
CCPA “precise geolocation” blur (~564 m)	$\varepsilon = \ln 2 / 564$	mean radius $\approx 1.6$ km
Bounded tail for on-map reports	clamp $p \le 1 - 10^{-12}$	caps max radius, breaks pure DP — account for it
Correct eastward offset	divide `dx_m` by $111{,}320\cos(\text{lat})$	required above ~45° latitude

Reference Implementation

The function below is a single, focused planar-Laplace sampler. It uses scipy.special.lambertw on the $-1$ branch for the radius inverse-CDF — a hardened, well-tested evaluation that avoids hand-rolling a Halley iteration near the awkward branch point. Every line that touches the privacy guarantee is commented.

python

from __future__ import annotations

import math
import random
from typing import Tuple

from scipy.special import lambertw

# Metres per degree of latitude on the WGS84 ellipsoid (mean value).
_M_PER_DEG_LAT: float = 111_320.0


def sample_planar_laplace(
    lat: float,
    lon: float,
    epsilon: float,
    rng: random.Random | None = None,
) -> Tuple[float, float]:
    """Perturb one WGS84 point with 2D planar-Laplace noise.

    The reported point is drawn from the planar Laplace density
    D(x) = (epsilon**2 / (2*pi)) * exp(-epsilon * ||x||) centred on the true
    location. `epsilon` is a PER-METRE privacy parameter: two points d metres
    apart are (epsilon*d)-indistinguishable, so smaller epsilon means more
    noise and stronger privacy. The mean reported displacement is 2/epsilon
    metres -- the quantity to reason about operationally.

    Sampling uses the polar decomposition, in which the two coordinates are
    independent:
      * angle  theta ~ Uniform(0, 2*pi)                -- isotropic, leaks no bearing
      * radius r from density proportional to r*exp(-epsilon*r), obtained by
        inverting its CDF  C(r) = 1 - (1 + epsilon*r) * exp(-epsilon*r)  with
        the Lambert W function on its -1 branch.

    Returns the perturbed (lat, lon) in decimal degrees. Adding independent 1-D
    Laplace noise per axis would NOT reproduce this density and would leak the
    bearing of the true point -- do not substitute it.
    """
    if epsilon <= 0.0:
        raise ValueError("epsilon must be > 0 (per-metre privacy parameter)")
    if not (-90.0 <= lat <= 90.0) or not (-180.0 <= lon <= 180.0):
        raise ValueError(f"coordinate out of WGS84 range: ({lat}, {lon})")
    r_source = rng if rng is not None else random

    # 1. Isotropic bearing: no direction is preferred, so nothing about which
    #    way the true point lies is revealed.
    theta: float = r_source.uniform(0.0, 2.0 * math.pi)

    # 2. Uniform quantile, clamped away from the endpoints. p -> 1 sends the
    #    radius to infinity (heavy tail); p -> 0 lands on the Lambert W branch
    #    point -1/e where evaluation loses precision. Both guards are numerical,
    #    not statistical -- the clamp perturbs the distribution negligibly.
    p: float = min(max(r_source.random(), 1e-12), 1.0 - 1e-12)

    # 3. Invert the radial CDF. Solving 1 - (1 + eps*r)*exp(-eps*r) = p gives
    #    r = -(1/eps) * (W_{-1}((p-1)/e) + 1). The -1 branch is REQUIRED: since
    #    1 + eps*r >= 1 the root satisfies W <= -1, which only W_{-1} returns;
    #    the principal branch would hand back a wrong (negative-radius) root.
    w: float = lambertw((p - 1.0) / math.e, k=-1).real  # imag part is ~0, discard
    r: float = -(1.0 / epsilon) * (w + 1.0)             # radius in METRES
    r = max(r, 0.0)                                      # guard tiny negatives

    # 4. Metric offset (metres) -> degrees at the LOCAL scale. One degree of
    #    longitude shrinks by cos(lat); using the local factor keeps the
    #    eastward error correct away from the equator. The cos guard avoids a
    #    divide-by-zero at the poles (where longitude is degenerate anyway).
    dx_m: float = r * math.cos(theta)   # eastward component
    dy_m: float = r * math.sin(theta)   # northward component
    d_lat: float = dy_m / _M_PER_DEG_LAT
    cos_lat: float = math.cos(math.radians(lat))
    d_lon: float = dx_m / (_M_PER_DEG_LAT * max(cos_lat, 1e-12))

    return lat + d_lat, lon + d_lon

The mechanism is memoryless per query: it does not compose across repeated reports of the same subject. If a device emits a location every few seconds, correlated planar-Laplace samples steadily erode the guarantee, which is why streaming pipelines route through w-event privacy for trajectories rather than calling this sampler per fix.

Validation Checkpoint

Planar-Laplace bugs — the wrong Lambert $W$ branch, a forgotten $\cos(\text{lat})$ , per-axis noise masquerading as planar — do not raise; they return plausible points with a subtly wrong distribution. The two invariants worth trusting are that the noise is unbiased (mean offset $\approx 0$ ) and that the empirical radius mean matches $2/\varepsilon$ . The harness below asserts both over a large sample and is safe to run in CI.

python

def _validate() -> None:
    import statistics

    lat0, lon0 = 48.8566, 2.3522          # Paris
    epsilon = math.log(2.0) / 200.0        # ln 2 privacy at 200 m; mean radius ~577 m
    n = 200_000
    rng = random.Random(20260715)          # deterministic for CI reproducibility

    lat_off: list[float] = []
    lon_off: list[float] = []
    radii: list[float] = []
    for _ in range(n):
        plat, plon = sample_planar_laplace(lat0, lon0, epsilon, rng=rng)
        d_lat = plat - lat0
        d_lon = plon - lon0
        lat_off.append(d_lat)
        lon_off.append(d_lon)
        # reconstruct the metric radius using the same local scale
        dy = d_lat * _M_PER_DEG_LAT
        dx = d_lon * _M_PER_DEG_LAT * math.cos(math.radians(lat0))
        radii.append(math.hypot(dx, dy))

    # 1. Radius law: empirical mean matches 2/epsilon within 2% (Gamma(2,1/eps)).
    mean_r = statistics.fmean(radii)
    expected_r = 2.0 / epsilon
    assert abs(mean_r - expected_r) / expected_r < 0.02, (mean_r, expected_r)

    # 2. Unbiasedness: each axis mean offset ~ 0. Per-axis std is sqrt(2)/eps,
    #    so the sample mean sits within ~5 sigma / sqrt(n) of zero -- convert
    #    that generous bound back to metres for the assertion.
    axis_tol_m = 5.0 * (math.sqrt(2.0) / epsilon) / math.sqrt(n)
    mean_dx = statistics.fmean(lon_off) * _M_PER_DEG_LAT * math.cos(math.radians(lat0))
    mean_dy = statistics.fmean(lat_off) * _M_PER_DEG_LAT
    assert abs(mean_dx) < axis_tol_m, mean_dx
    assert abs(mean_dy) < axis_tol_m, mean_dy

    print(f"radius mean {mean_r:.1f} m vs expected {expected_r:.1f} m; "
          f"bias dx={mean_dx:.2f} m dy={mean_dy:.2f} m -- all invariants hold")


if __name__ == "__main__":
    _validate()

A stronger test — worth adding when you touch the branch selection — bins the radii and checks the empirical CDF against $C(r) = 1 - (1 + \varepsilon r)e^{-\varepsilon r}$ with a Kolmogorov–Smirnov statistic; a wrong-branch bug passes neither the mean check nor the KS check.

Incident Response and Edge Cases

Wrong Lambert $W$ branch. Calling lambertw(z) without k=-1 uses the principal branch $W_0$ , which on the interval $(-1/e, 0)$ returns values in $[-1, 0]$ and yields a near-zero or negative radius. Symptom: reported points cluster on top of the true location and the radius mean falls far below $2/\varepsilon$ . Remediation: pin k=-1, and keep the radius-mean assertion in CI so a regression fails the build rather than shipping a mechanism with almost no noise.
Heavy-tail blow-up. Because the radial density has an unbounded tail, an unclamped quantile near $1$ can place a report thousands of kilometres away — across an ocean, or as inf/NaN if $W$ overflows. Symptom: occasional absurd coordinates or non-finite output. Remediation: clamp $p \le 1 - 10^{-12}$ as shown; if you additionally truncate the radius at a hard maximum to keep points on-map, record that you have done so, because truncation biases the distribution and breaks pure geo-indistinguishability — treat it as an approximate mechanism and account for the relaxation.
Longitude scaling at high latitude. Converting the eastward offset with the latitude factor ( $111{,}320$ m/deg) instead of $111{,}320\cos(\text{lat})$ makes real displacements smaller than intended east–west, weakening privacy in that axis. Symptom: the empirical radius mean drifts below $2/\varepsilon$ and grows with latitude. Remediation: always divide dx_m by $\cos(\text{lat})$ ; for work above ~60° or spanning wide areas, project to a local metric CRS (UTM), sample the offset there, and reproject.
Antimeridian and pole wraparound. A perturbation near $\pm180°$ longitude or $\pm90°$ latitude can push the reported point off the valid range. Symptom: longitudes like $180.4$ or latitudes above $90$ leak into downstream joins. Remediation: wrap longitude modulo $360$ back into $[-180, 180]$ ; do not clamp latitude at the pole, since hard clamping concentrates probability mass there and leaks that the true point was polar — reflect across the pole or handle high-latitude data in a projected CRS instead.

Frequently Asked Questions

Why must I use the -1 branch of the Lambert W function?

The radius satisfies 1 + epsilon*r >= 1, so when you rearrange the CDF equation the Lambert W argument must map to a value less than or equal to -1. Only the -1 branch (k=-1) produces values in that range; the principal branch W_0 returns values in [-1, 0] on the same input interval and therefore hands back a wrong root, typically a near-zero or negative radius. The bug is silent because the output still looks like a coordinate, which is why the radius-mean check belongs in CI.

What does epsilon mean as a per-metre parameter, and how do I pick it?

Unlike the scalar epsilon of the standard Laplace mechanism, here epsilon has units of inverse metres and sets how fast indistinguishability decays with distance: two points d metres apart are (epsilon*d)-indistinguishable. Pick it indirectly by choosing an indistinguishability level l (for example ln 2) and the radius R at which it should hold, then set epsilon = l / R. The mean reported displacement is 2/epsilon metres, so always sanity-check that value against your map before deploying.

Why sample the radius from r*exp(-epsilon*r) instead of a plain exponential?

The planar Laplace density is exp(-epsilon*r) in the plane, but converting to polar coordinates multiplies by the Jacobian r (the circumference factor). The marginal density of the radius is therefore proportional to r*exp(-epsilon*r), a Gamma(2, 1/epsilon) distribution with mean 2/epsilon, not an exponential with mean 1/epsilon. Sampling a plain exponential radius concentrates points too close to the centre and produces the wrong displacement distribution.

Can I just add independent Laplace noise to latitude and longitude separately?

No. Independent per-axis Laplace noise produces a diamond-shaped, anisotropic density whose contours leak the bearing of the true point and does not satisfy geo-indistinguishability. The planar Laplace is radially symmetric by construction; you get that symmetry only by sampling a uniform angle and a radius from the r*exp(-epsilon*r) law, then projecting the offset onto the two axes.

Spatial Noise Mechanisms — the parent guide covering the full family of location and aggregate noise mechanisms this sampler belongs to.
Geo-Indistinguishability for Location Privacy — the privacy semantics and epsilon-per-metre definition that this distribution enforces.
Implementing Geo-Indistinguishability in Python — wiring this sampler into an end-to-end obfuscation service with budgeting.
Calibrating Gaussian Noise for Spatial Aggregates — the sibling mechanism for counts and heatmaps rather than single points.
w-Event Privacy for Streaming Trajectories — how to compose this per-point noise safely across a continuous stream.

Up: Spatial Noise Mechanisms · Differential Privacy for Geospatial Data