Fast implementation of the moving quantile based on the partitioning heaps, Part 3
Previous posts explored fast implementations of moving quantiles based on partitioning heaps:
Fast implementation of the moving quantile based on the partitioning heaps, Part 1
By Andrey Akinshin
·
2020-12-29
The Hardle-Steiger method to estimate the moving median and its generalization for the moving quantilesPart 1,
Fast implementation of the moving quantile based on the partitioning heaps, Part 2
By Andrey Akinshin
·
2021-01-19
Improvements of the Hardle-Steiger method that allows estimating moving quantiles using linear interpolationPart 2.
This post presents a cleaner implementation that is simpler to use.
Compared to previous attempts, this version makes several important changes:
- Quantile-focused API: The approach estimates only the target quantile, dropping the $k^\textrm{th}$ order statistics approach
- Consistent quantile evaluation: The algorithm maintains quantile estimates throughout initialization, even when observed elements are less than the window size
- Linear interpolation only: All quantile estimates use linear interpolation rather than single order statistics
- Hyndman-Fan Type 7 standard: The implementation uses the Type 7 quantile estimator from Hyndman–Fan TaxonomyHyndman–Fan Taxonomy for consistency with statistical software
- Invariant maintenance: Elements distribute between lower and upper heaps so linear interpolation occurs between the root and the first upper heap element, preserving this relationship when the window fills
The following sections provide detailed algorithm descriptions and a reference implementation in C#.
Algorithm Description
Consider the problem of tracking quantiles in streaming data. When new values arrive continuously, applications need the 75th percentile, median, or any other quantile of the most recent 1000 values without storing and sorting the entire window every time. Quantiles work especially well for noisy data because they resist outliers — a single extreme value cannot distort the entire estimate, unlike moving averages. The partitioning heaps algorithm solves this problem efficiently.
The core insight comes from a simple observation: to calculate any quantile, the algorithm needs at most two specific order statistics from the data. For the median of 9 values, the algorithm needs the 5th smallest element. For the 75th percentile, the calculation may need interpolation between the 7th and 8th smallest elements. Instead of sorting the entire window repeatedly, the partitioning heaps algorithm maintains these crucial elements in their correct relative positions while allowing efficient updates.
This extends the original Hardle-Steiger method, which was designed only for medians. By adjusting the heap sizes dynamically, the same data structure can track any quantile efficiently, extending the algorithm far beyond its original median-smoothing purpose.
The Fundamental Idea
Consider a row of elements arranged by value, where only one specific position matters for the target quantile. The algorithm creates a specialized data structure where this target position becomes the “root” of two connected heaps. Elements to the left (smaller values) form a max-heap, and elements to the right (larger values) form a min-heap.
This arrangement has a useful property: the root element and its immediate right neighbor in the upper heap are always the two order statistics needed for linear interpolation. The algorithm preserves these relationships regardless of how data shifts and changes.
Core Data Structure
The algorithm organizes all elements in a single array that logically functions as two connected heaps:
Lower heap (max-heap) | Root | Upper heap (min-heap)
[0] ... [rootIndex-1] | [rootIndex] | [rootIndex+1] ... [windowSize-1]
The lower heap is a max-heap containing elements smaller than the root.
It occupies positions from 0 to rootHeapIndex - 1.
The parent is always greater than or equal to its children, so larger values bubble up toward the root boundary.
The root element serves as the pivot that separates the two heaps.
It sits at position rootHeapIndex, calculated using:
rootHeapIndex = floor((windowSize - 1) * probability)
This calculation ensures the correct proportion of elements falls on each side.
The upper heap is a min-heap containing elements larger than the root.
It spans positions from rootHeapIndex + 1 to windowSize - 1.
Here, the parent is always less than or equal to its children, so smaller values bubble up toward the root boundary.
This design ensures that elements naturally segregate around the target quantile position, with heap properties maintaining the necessary ordering without full sorting.
The heap indexing uses a specific pattern for efficient parent-child navigation.
For any element at heap position i, the parent sits at position rootHeapIndex + (i - rootHeapIndex) / 2.
Children of element i are located at positions rootHeapIndex + (i - rootHeapIndex) * 2 and rootHeapIndex + (i - rootHeapIndex) * 2 + sign(i - rootHeapIndex).
This indexing scheme lets the sift operations traverse both heap types with the same logic.
Why This Structure Works
This approach works by using heap properties to maintain quantile relationships. Traditional heaps efficiently find the minimum or maximum, but not arbitrary order statistics. By creating two heaps that meet at the target quantile position, the algorithm gets direct access to the elements needed for quantile calculation.
The partitioning property is essential: all elements in the lower heap are less than or equal to the root, which is less than or equal to all elements in the upper heap. This creates a natural ordering where the root represents the target quantile position, and the first element of the upper heap gives the next order statistic for interpolation.
The heap properties automatically maintain these relationships when elements are added or removed. A newly inserted value sifts to its appropriate position, potentially displacing other elements, but the fundamental partitioning around the target quantile remains intact.
Element Addition and Window Management
When a new element arrives, the algorithm handles two scenarios. For unfilled windows, it places the element in the appropriate heap while building toward the target proportion. For full windows, it replaces the oldest element and restores the heap properties.
For unfilled windows, the algorithm distributes elements between heaps intelligently. The algorithm calculates the desired lower heap size:
floor((totalElements - 1) * probability)
New elements are added to whichever heap needs to grow. This ensures that even during the initialization phase, the quantile estimates remain meaningful and approach the target quantile as more elements arrive.
For full windows, the replacement process requires careful management. The algorithm maintains two mapping arrays: one mapping heap positions to original elements, and another mapping original elements to heap positions. When replacing the oldest element, the algorithm knows exactly where that element lives in the heap structure and can overwrite it directly, then restore heap properties.
This mapping system makes the algorithm practical. Without these mappings, finding and removing an arbitrary element from a heap would require expensive search operations. With mappings, element replacement becomes a constant-time operation followed by logarithmic heap maintenance.
The removal process follows what the original paper calls “hole propagation.” When an old element is removed, it creates a “hole” in the heap structure. The algorithm fills this hole by moving the replacement element into position, then sifting to restore heap properties. This approach avoids the expensive search operations that would otherwise be required to locate and remove arbitrary elements from heap structures.
The Sift Operation: Restoring Order
After inserting a new value, heap properties may be violated. The sift operation restores these properties. This operation is more complex than traditional heap operations because it handles two different heap types connected at a root.
The sift process starts at the insertion point and moves in two phases. First, the operation sifts upward: if the newly inserted value violates the heap property with respect to its parent, the algorithm swaps their positions and continues up the tree. This handles cases where a large value is inserted into the lower heap or a small value into the upper heap.
If upward sifting fails to resolve all violations, the algorithm sifts downward. The algorithm compares the current element with its children and swaps with the child that best satisfies the heap property. For the lower heap (max-heap), this means swapping with the larger child. For the upper heap (min-heap), it means swapping with the smaller child.
Complexity arises at the root position, which connects both heaps. The root can be displaced by elements from either side, requiring careful handling of these boundary conditions. The implementation uses a nested helper function to manage different cases of child relationships and heap types.
The sift operation is convergent. The sift process eventually finds the correct element position regardless of insertion location or initial heap property violations. This robustness keeps quantile relationships intact as the window slides.
Quantile Calculation Through Linear Interpolation
With the heap structure set up, calculating the quantile is straightforward. The algorithm uses the Hyndman-Fan Type 7 method, the standard in most statistical software including R, Python’s numpy, and Excel.
The root element represents one boundary of the quantile estimate.
If elements exist in the upper heap, the first element at position rootHeapIndex + 1 provides the other boundary.
These two elements provide the order statistics needed for linear interpolation.
Here’s how the interpolation formula works.
For a window of size n and probability p, the theoretical position is:
h = (n - 1) * p
If this position falls between two elements at positions floor(h) and ceil(h), the quantile becomes:
quantile = element[floor(h)] + (h - floor(h)) * (element[ceil(h)] - element[floor(h)])
This interpolation gives smooth quantile estimates instead of the discrete jumps you get from simple order statistics. When the theoretical position falls exactly on an integer, the formula returns that specific order statistic. When it falls between positions, the formula gives a weighted average that reflects how close the desired quantile is to each boundary element.
Handling Edge Cases and Small Windows
The algorithm handles several edge cases that can complicate quantile estimation. When the window contains only one element, the quantile equals that element. When the upper heap is empty (which occurs with extreme quantiles or small windows), the quantile equals the root element.
For small windows during initialization, the algorithm maintains meaningful estimates throughout the process. Instead of waiting until the window fills, the algorithm calculates quantiles based on elements seen so far. This helps applications that need immediate quantile estimates, even with limited data.
The mapping between probability and heap sizes adjusts automatically as the window fills. For a 0.5 quantile (median) in a window of size 21, the root sits at position 10, with 10 elements in each heap. For a 0.9 quantile in the same window, the root moves to position 18, with 18 elements in the lower heap and only 2 in the upper heap. This flexibility lets the same algorithm handle any quantile efficiently.
Performance Characteristics and Trade-offs
The algorithm achieves O(log windowSize) complexity for element addition, which is optimal for heap-based approaches. This logarithmic scaling means doubling the window size adds only one more heap level, keeping updates efficient even for large windows.
Memory usage is O(windowSize), storing exactly the elements needed plus mapping structures. This eliminates the overhead of sorting or maintaining multiple data structures. This makes the algorithm suitable for memory-constrained environments.
Quantile retrieval runs in O(1) because necessary elements remain available at known positions. This constant-time access matters for high-frequency applications where quantile values are queried much more often than new data arrives.
The algorithm works best for moderate window sizes (roughly 10 to 10,000 elements). For very small windows, simpler approaches may suffice. For very large windows, more sophisticated techniques like approximate quantile sketches may be more appropriate. Within this range, the partitioning heaps approach balances accuracy, efficiency, and simplicity well.
This method has clear advantages over alternatives. Naive recomputation by sorting the entire window requires O(n log n) time per update. Maintaining a sorted list with insertions and deletions requires O(n) time for pointer maintenance despite O(log n) binary search. The partitioning heaps approach achieves the theoretically optimal O(log n) complexity while keeping exact quantile values instead of approximations.
Reference Implementation
/// <summary>
/// A moving selector based on a partitioning heaps.
/// Memory: O(windowSize).
/// Add complexity: O(log(windowSize)).
/// GetValue complexity: O(1).
///
/// <remarks>
/// Based on the following paper:
/// Hardle, W., and William Steiger. "Algorithm AS 296: Optimal median smoothing." Journal of the Royal Statistical Society.
/// Series C (Applied Statistics) 44, no. 2 (1995): 258-264.
/// <br /><br />
/// (c) 2020–2025 Andrey Akinshin, MIT License
/// </remarks>
/// </summary>
public class PartitioningHeapsMovingQuantileEstimator
{
private readonly int windowSize;
private readonly double probability;
private readonly double[] heap;
private readonly int[] heapToElementIndex;
private readonly int[] elementToHeapIndex;
private readonly int rootHeapIndex;
private int lowerHeapSize, upperHeapSize, totalElementCount;
public PartitioningHeapsMovingQuantileEstimator(int windowSize, double probability)
{
if (windowSize <= 0)
throw new ArgumentOutOfRangeException(nameof(windowSize), windowSize,
$"{nameof(windowSize)} must be positive");
if (probability < 0 || probability > 1 || !double.IsFinite(probability))
throw new ArgumentOutOfRangeException(nameof(probability), probability,
$"{nameof(probability)} must be between 0 and 1");
this.windowSize = windowSize;
this.probability = probability;
heap = new double[windowSize];
heapToElementIndex = new int[windowSize];
elementToHeapIndex = new int[windowSize];
rootHeapIndex = (int)Math.Floor((windowSize - 1) * probability);
}
private void Swap(int heapIndex1, int heapIndex2)
{
int elementIndex1 = heapToElementIndex[heapIndex1];
int elementIndex2 = heapToElementIndex[heapIndex2];
double value1 = heap[heapIndex1];
double value2 = heap[heapIndex2];
heap[heapIndex1] = value2;
heap[heapIndex2] = value1;
heapToElementIndex[heapIndex1] = elementIndex2;
heapToElementIndex[heapIndex2] = elementIndex1;
elementToHeapIndex[elementIndex1] = heapIndex2;
elementToHeapIndex[elementIndex2] = heapIndex1;
}
private void Sift(int heapIndex)
{
int SwapWithChildren(int heapCurrentIndex, int heapChildIndex1, int heapChildIndex2, bool isUpperHeap)
{
bool hasChild1 = rootHeapIndex - lowerHeapSize <= heapChildIndex1 &&
heapChildIndex1 <= rootHeapIndex + upperHeapSize;
bool hasChild2 = rootHeapIndex - lowerHeapSize <= heapChildIndex2 &&
heapChildIndex2 <= rootHeapIndex + upperHeapSize;
if (!hasChild1 && !hasChild2)
return heapCurrentIndex;
if (hasChild1 && !hasChild2)
{
if (heap[heapIndex] < heap[heapChildIndex1] && !isUpperHeap ||
heap[heapIndex] > heap[heapChildIndex1] && isUpperHeap)
{
Swap(heapIndex, heapChildIndex1);
return heapChildIndex1;
}
return heapCurrentIndex;
}
if (hasChild1 && hasChild2)
{
if ((heap[heapIndex] < heap[heapChildIndex1] || heap[heapIndex] < heap[heapChildIndex2]) &&
!isUpperHeap ||
(heap[heapIndex] > heap[heapChildIndex1] || heap[heapIndex] > heap[heapChildIndex2]) && isUpperHeap)
{
int heapChildIndex0 =
heap[heapChildIndex1] > heap[heapChildIndex2] && !isUpperHeap ||
heap[heapChildIndex1] < heap[heapChildIndex2] && isUpperHeap
? heapChildIndex1
: heapChildIndex2;
Swap(heapIndex, heapChildIndex0);
return heapChildIndex0;
}
return heapCurrentIndex;
}
throw new InvalidOperationException();
}
while (true)
{
if (heapIndex != rootHeapIndex)
{
bool isUpHeap = heapIndex > rootHeapIndex;
int heapParentIndex = rootHeapIndex + (heapIndex - rootHeapIndex) / 2;
if (heap[heapParentIndex] < heap[heapIndex] && !isUpHeap ||
heap[heapParentIndex] > heap[heapIndex] && isUpHeap)
{
Swap(heapIndex, heapParentIndex);
heapIndex = heapParentIndex;
continue;
}
else
{
int heapChildIndex1 = rootHeapIndex + (heapIndex - rootHeapIndex) * 2;
int heapChildIndex2 = rootHeapIndex + (heapIndex - rootHeapIndex) * 2 +
Math.Sign(heapIndex - rootHeapIndex);
int newHeapIndex = SwapWithChildren(heapIndex, heapChildIndex1, heapChildIndex2, isUpHeap);
if (newHeapIndex != heapIndex)
{
heapIndex = newHeapIndex;
continue;
}
}
}
else // heapIndex == rootHeapIndex
{
if (lowerHeapSize > 0)
{
int newHeapIndex = SwapWithChildren(heapIndex, heapIndex - 1, -1, false);
if (newHeapIndex != heapIndex)
{
heapIndex = newHeapIndex;
continue;
}
}
if (upperHeapSize > 0)
{
int newHeapIndex = SwapWithChildren(heapIndex, heapIndex + 1, -1, true);
if (newHeapIndex != heapIndex)
{
heapIndex = newHeapIndex;
continue;
}
}
}
break;
}
}
public void Add(double value)
{
int elementIndex = totalElementCount % windowSize;
int Insert(int heapIndex)
{
heap[heapIndex] = value;
heapToElementIndex[heapIndex] = elementIndex;
elementToHeapIndex[elementIndex] = heapIndex;
return heapIndex;
}
if (totalElementCount++ < windowSize) // Heap is not full
{
if (totalElementCount == 1) // First element
{
Insert(rootHeapIndex);
}
else
{
int desiredLowerHeapSize = (int)Math.Floor((totalElementCount - 1) * probability);
if (lowerHeapSize < desiredLowerHeapSize)
{
lowerHeapSize++;
int heapIndex = Insert(rootHeapIndex - lowerHeapSize);
Sift(heapIndex);
}
else
{
upperHeapSize++;
int heapIndex = Insert(rootHeapIndex + upperHeapSize);
Sift(heapIndex);
}
}
}
else
{
// Replace old element
int heapIndex = elementToHeapIndex[elementIndex];
Insert(heapIndex);
Sift(heapIndex);
}
}
public double Quantile()
{
if (totalElementCount == 0)
throw new InvalidOperationException();
if (totalElementCount == 1 || upperHeapSize == 0)
return heap[rootHeapIndex];
double a = heap[rootHeapIndex];
double b = heap[rootHeapIndex + 1];
int n = Math.Min(totalElementCount, windowSize);
double h = (n - 1) * probability;
return a + (h - Math.Floor(h)) * (b - a);
}
}
using Xunit.Abstractions;
public class PartitioningHeapsMovingQuantileEstimatorTests(ITestOutputHelper output)
{
private const int MaxN = 1000;
[Fact]
public void MainTest()
{
var random = new Random(1729);
for (int windowSize = 1; windowSize < 20; windowSize++)
{
output.WriteLine($"WindowSize = {windowSize}");
int m = Math.Max(3, windowSize + 3 * (windowSize - 1));
double[] probabilities = new double[m];
var estimators = new PartitioningHeapsMovingQuantileEstimator[m];
for (int k = 0; k < m; k++)
{
probabilities[k] = 1.0 * k / (m - 1);
estimators[k] = new PartitioningHeapsMovingQuantileEstimator(windowSize, probabilities[k]);
}
double[] sample = new double[MaxN];
for (int n = 0; n < MaxN; n++)
{
output.WriteLine($" n = {n + 1}");
sample[n] = random.NextDouble();
double[] window = sample.Take(n + 1).TakeLast(windowSize).ToArray();
Array.Sort(window);
for (int k = 0; k < m; k++)
{
var estimator = estimators[k];
double p = probabilities[k];
output.WriteLine($" p = {p}");
estimator.Add(sample[n]);
double actual = estimator.Quantile();
double expected = GetTrueQuantile(window, p);
Assert.Equal(expected, actual, 9);
}
}
}
}
// Assuming sample is sorted
private static double GetTrueQuantile(double[] sample, double p)
{
if (sample.Length == 1)
return sample[0];
double h = (sample.Length - 1) * p;
int hFloor = (int)Math.Floor(h);
int hCeil = (int)Math.Ceiling(h);
return sample[hFloor] + (h - hFloor) * (sample[hCeil] - sample[hFloor]);
}
}
References
Algorithm AS 296: Optimal Median Smoothing · 1995
· W. Hardle
et al.
Running Median Algorithm and Implementation for Integer Streaming Applications · 2019-06-01
· Oswaldo Cadenas
et al.
An Optimal Algorithm for Sliding Window Order Statistics · 2023-01-01
· Pavel Raykov