Running latency percentiles

Developing enterprise software has taught me that "good enough" needs a definition. In B2B SaaS, the promises in an SLA eventually become engineering constraints. Before relying on a metric, I want to understand what it tells me and what it leaves out.

A running average needs only a count and a sum. Following a percentile also requires information about how request durations are ordered. P² (P-squared) keeps a small model of that order in five markers.

I would consider P² when I need a few predetermined percentiles, the per-stream state budget is tight, and a heuristic estimate is acceptable. Each tracker summarizes the observations seen since its latest initialization or reset.

For a concrete storage budget, consider 20,000 requests per second in total across two workers, with four 64-bit duration measurements per request. Retaining those values for one accumulation period costs about 549 MiB after 15 minutes or 2.15 GiB after an hour. The table compares that exact history with eight approximate, per-worker P95 trackers.

The fixed-state row uses one classical P95 tracker per worker and measurement. In this model, each tracker stores three five-element arrays (height, position, and fraction) with 64-bit entries, plus a 64-bit count: 3 * 5 * 8 + 8 = 128 bytes. Eight trackers total 8 * 128 = 1,024 bytes. Desired positions are recomputed, and the height array doubles as the startup buffer.

These are payload calculations for local estimates. Container, allocator, synchronization, and other runtime costs would add to both rows. Storage inputs.

Background

Consider two batches of 1,000 requests. In the first, every request takes 100 ms. In the second, 900 take 10 ms and the remaining 100 take 910 ms. Both average 100 ms, although the second batch has no request anywhere near that duration.

The average is correct. It just does not describe how the waiting time is spread across requests. To find P95, sort each batch from fastest to slowest and inspect its 950th request. The answers are 100 ms and 910 ms. This is the nearest-rank definition: for N > 0 observations and fraction 0 < p <= 1, take position ceil(p*N), counting from one.

That ordering information is what the percentile calculation needs.

Keeping all the values would solve that problem, but a service handling 20,000 requests per second accumulates 18 million observations in 15 minutes. One eight-byte duration per request takes roughly 137 MiB; four such measurements take about 549 MiB.

Instead, we can retain a small model of that ordering and update the model as each value arrives.

P²

Five markers for one requested percentile

For P95, p = 0.95 identifies the part of the ordered population we want to follow. The answer is a duration, perhaps 120 ms. The percentage specifies a target in the ranked population.

P² represents that ordering with a small set of markers. Each marker follows a fixed target fraction, written f, and maintains an estimated duration for that target. Those estimates are the state we carry forward.

The classical P² algorithm chooses five target fractions for a requested percentile 0 < p < 1: 1

0, p/2, p, (1+p)/2, 1

For a P95 tracker:

Here, p is the percentile requested from the tracker, while f is the target assigned to an individual marker. The lower support sits halfway between 0% and 95%; the upper support sits halfway between 95% and 100%. The midpoint calculations happen on the rank-fraction axis.

Each marker has a height, which is a duration estimate, and a current position, which is the rank attached to that estimate by the algorithm. A height of 120 ms is a value; a position of 94 is a place in the ordering. Positions count from one, while the marker indices in the table count from zero.

It also has a desired position: where its fixed target fraction says it should be after N observations:

desired_position = 1 + f * (N - 1)

This is the paper's and Boost's one-based position convention. Akinshin's zero-based version omits the leading 1. The two describe the same locations with a one-position offset.

For P95 at N = 100, the desired position is 95.05. A marker at current position 95 is 0.05 ranks behind that target. Desired positions move smoothly with the count; corrections change current positions in whole-rank steps.

The current position is maintained rank bookkeeping for an interpolated height, which can fall between observed durations. The exact nearest-rank definition gives us a reference against which to compare the estimate.

The third marker, at index 2, supplies P95. A height of 120 ms is the current P95 estimate. The other interior markers capture the local shape on either side so we can adjust it. Each marker's height and rank position evolve while its target fraction stays fixed.

Initialization

The classical initialization retains the first five values, sorts them, and uses them as the marker heights. Their current positions start at 1 through 5.

The first five observations form the startup phase. Until initialization is complete, an API can return an explicit pending status or an exact nearest-rank value from the buffered observations. An empty stream needs an explicit empty result. I would make that behavior clear to callers.

The initial heights need time to adapt. With the first five values 10, 10, 10, 910, 910, the middle marker starts at 10 ms, even though the exact nearest-rank P95 is 910 ms. 2

New value

An incoming duration gives the tracker two jobs: account for one more observation, then check whether its estimates still sit near their target fractions.

Imagine inserting a duration into sorted order. Later entries gain a rank while their durations stay the same. P² mirrors that change in its marker bookkeeping. After extending the minimum or maximum where needed, it finds the interval containing the new value and increments the positions of the markers after that interval.

The desired positions advance too, because the population has grown. A P95 target advances by 0.95 per observation; its current position advances by either zero or one during insertion. The difference accumulates as drift. Once its magnitude reaches one rank, a correction can move the marker a step toward its target.

For a concrete example, suppose N = 100, the P95 marker has height 100 ms and current position 95, and its upper support has height 150 ms and position 98. The desired P95 position is 95.05. Consider two alternative observations from this same starting state. Either one brings the count to 101 and the desired P95 position to 96:

• A new duration of 80 ms falls below P95. The P95 marker's current position advances from 95 to 96 because the new value belongs before it. That matches its new desired position, so the height stays at 100 ms.
• A new duration of 120 ms falls between P95 and its upper support. The upper support advances to position 99. P95 stays at position 95, one rank behind its new target. The four-rank gap to its upper neighbor leaves room to move to position 96. Its new height comes from interpolation between the neighboring markers, which the next section works through.

The position shift records where the sample landed. The height adjustment keeps the marker following its chosen percentile. A correction moves one rank at a time, provided the neighboring rank gap leaves room. Example state.

The following pseudocode recomputes desired positions from the count. It assumes initialized state and finite valid input. Check that one more observation fits the count and precision limits before entering this update:

record(x):
    height[0] = min(height[0], x)
    height[4] = max(height[4], x)
    k = largest index in 0..3 where height[k] <= x
    position[k+1 .. 4] += 1
    count += 1
 
    for i in 1, 2, 3:
        desired = 1 + fraction[i] * (count - 1)
        drift = desired - position[i]
        step = sign(drift)
        gap = position[i+1] - position[i] if step > 0
              else position[i] - position[i-1]
 
        if abs(drift) >= 1 and gap > 1:
            proposed = parabolic(i, step)
            if height[i-1] < proposed and proposed < height[i+1]:
                height[i] = proposed
            else:
                height[i] = linear(i, step)
            position[i] += step

Extending the minimum before the cell search ensures that a qualifying k exists. The minimum marker's position stays at 1; the maximum marker's position advances to the new count. The gap check keeps interior positions strictly ordered. Process the interior markers in ascending order, so each correction sees the updates already made in that pass. 3

Now we need a duration estimate for that new position.

The parabolic step

For the height calculation, take a smaller example. Suppose three neighboring markers, after the insertion bookkeeping, have these values:

rank position:       4      6      8
height in ms:       10     30     90

The middle marker needs to advance from position 6 to 7. A straight step toward its upper neighbor would move halfway from 30 to 90 ms, giving 60 ms.

The straight step uses only the upper neighbor. P² first fits a parabola through both neighbors and the current marker, accounting for the change from 10 to 30 ms per rank in the adjacent slopes. Evaluating it at position 7 gives 55 ms, inside the neighboring heights, so it can be accepted.

Writing h[i] for height, r[i] for position, and s for the direction, either -1 or +1, evaluate the following interpolation arithmetic in floating point. This includes the intermediates s, A and B, even though their values represent whole-number steps and rank distances:

A = r[i] - r[i-1]
B = r[i+1] - r[i]
 
parabolic(i, s) = h[i] + s/(A+B) * (
    (A+s) * (h[i+1]-h[i])/B
  + (B-s) * (h[i]-h[i-1])/A
)

This includes s/(A+B): integer division would evaluate 1/(2+2) as zero and produce 30 ms instead of 55 ms. The division's operands must already be floating point; assigning its result to a floating-point height happens too late.

The curve can also overshoot. Keep the same positions, but change the heights to 10, 30, 31. The parabolic proposal at position 7 becomes 32.875 ms, above the upper neighbor's 31 ms.

The implementation rejects it and uses the linear step:

linear(i, +1) = h[i] + (h[i+1]-h[i]) / (r[i+1]-r[i])
linear(i, -1) = h[i] - (h[i]-h[i-1]) / (r[i]-r[i-1])

Both divisions are floating point, and both denominators are positive rank distances. The upward fallback gives 30.5 ms, preserving the ordering. For a downward example at the same ranks, heights 10, 11, 31 give a parabolic proposal of 8.125 ms at rank 5, below the lower neighbor; the fallback gives 10.5 ms. The fallback keeps the heights ordered through the same one-rank move. 3

Reading the estimate

Once the tracker is initialized, reading its configured percentile is simply reading:

estimate = height[2]

Reading is the cheap part: return that marker's height. The work happens as observations arrive.

Tracking more than one percentile

For P90, P95 and P99, three independent classical trackers are the simplest extension: each receives the same duration, but uses its own target fractions. That is 15 marker heights plus their bookkeeping, still independent of stream length.

Extended P² can share markers across requested percentiles. For three targets, one established layout uses nine markers:

0%, 45%, 90%, 92.5%, 95%, 97%, 99%, 99.5%, 100%

The requested values then sit at indices 2, 4 and 6. The grid follows the requested targets and their midpoints. Desired positions use those fractions, and reads use the matching target indices. Shared markers stay ordered by construction; independent trackers can produce crossing estimates. 4

Reliability

A p95 field is the visible result. Behind it, the tracker maintains a small amount of related state.

Desired positions can be stored and incremented, or recomputed from the count as above. Repeated floating-point additions can accumulate rounding drift. Recomputing avoids that accumulation; the approximation error still depends on the data and its arrival order. 5

A classical tracker can reuse its five-height array during startup, as in the storage model above. An implementation that retains desired positions or a separate startup buffer needs to include those fields in its byte budget.

The invariants are small enough to write down and enforce:

heights are nondecreasing
rank positions are strictly increasing
the endpoints track the minimum and maximum
the observation count advances once per accepted value

Repeated durations must be allowed: equal heights are valid.

Finally, an update changes several fields together. Give each tracker a single owner, or synchronize reads and updates as a unit. An atomic count alone leaves readers exposed to a mix of old and new marker state.

Conclusion

For a fixed set of targets, P² keeps storage and per-observation work bounded. Each value updates the ordering information, corrects drifting markers, and can then be discarded. Reading the estimate becomes a marker lookup.

How useful that estimate is depends on the population, arrival order, and sample count. I would check its error on representative workloads before relying on it. Resetting starts a fresh estimate with its own warm-up. For rolling-window queries or a service-wide view across workers, I would use the metric system's histogram or mergeable-sketch support. 6

References

1. Raj Jain and Imrich Chlamtac, The P² Algorithm for Dynamic Calculation of Quantiles and Histograms Without Storing Observations, 1985; Andrey Akinshin's implementation-focused introduction.
2. Andrey Akinshin, P² initialization strategy.
3. Boost.Accumulators, classical P² implementation.
4. Boost.Accumulators, extended P² implementation.
5. Andrey Akinshin, P² rounding issue.
6. Prometheus, Histograms and summaries, especially aggregation and time-window trade-offs.