Efficient Quantile Tracking Using an Oracle

For incremental quantile estimators the step size and possibly other tuning parameters must be carefully set. However, little attention has been given on how to set these values in an online manner. In this article we suggest two novel procedures that address this issue. The core part of the procedures is to estimate the current tracking mean squared error (MSE). The MSE is decomposed in tracking variance and bias and novel and efficient procedures to estimate these quantities are presented. It is shown that estimation bias can be tracked by associating it with the portion of observations below the quantile estimates. The first procedure runs an ensemble of $L$ quantile estimators for wide range of values of the tuning parameters and typically around $L = 100$. In each iteration an oracle selects the best estimate by the guidance of the estimated MSEs. The second method only runs an ensemble of $L = 3$ estimators and thus the values of the tuning parameters need from time to time to be adjusted for the running estimators. The procedures have a low memory foot print of $8L$ and a computational complexity of $8L$ per iteration. The experiments show that the procedures are highly efficient and track quantiles with an error close to the theoretical optimum. The Oracle approach performs best, but comes with higher computational cost. The procedures were further applied to a massive real-life data stream of tweets and proofed real world applicability of them.


Introduction
The volumes of automatically generated data are constantly increasing (Ramírez-Gallego et al.; with more urgent demand for being analyzed in real-time (Kejariwal et al.;2015). Conventional statistical and data mining techniques are constructed for offline situations and therefore are often not applicable for such real-time analysis (Krempl et al.;. Thus a wide range of streaming algorithms are continuously being developed by the data mining community addressing a range of real-time tasks such as clustering, filtering, cardinality estimation, estimation of moments or quantiles, predictions, dimensionality reduction and anomaly detection (Kejariwal et al.;2015).
An inherent part of almost any real-time method is how rapid the history of the data stream should be forgotten. The forgetting mechanism should be adjusted depending on the properties of the data stream which again can vary over time. A prominent example to address this issue is the adaptive windowing approach (ADWIN) (Bifet and Gavalda; where a historic and a current window of observations are compared to look for changes in the properties of the data stream. The sizes of the windows are further adapted, and properties of the data stream can be computed from the adaptive windows. Costa and Vázquez-Abad (2007) suggests methods to adjust the step sizes of stochastic approximation algorithms related to the frequency of regime switches. There is also a wide range of event detection based methods that rapidly adjust the current estimate when an event is detected, see e.g. Ross (2014); Hammer and Yazidi (2018).
In this paper we address the issue when the objective is to track quantiles of data streams. To the best of our knowledge this is the first article addressing this. Quantiles are useful to characterize data stream distribution in a flexible and non-parametric way (Luo et al.; and have been used for a wide range of applications such as portfolio risk measurement in the stock market (Gilli et al.;2006;Abbasi and Guillen;, fraud detection (Zhang and Guan;2008), signal processing and filtering (Stahl et al.;2000), climate change monitoring 2011), SLA violation monitoring , 2010, network monitoring (Choi et al.;Liu et al.;, Monte Carlo simulation , structural health monitoring (Gregory et al.;, non-parametric statistical testing (Lall;2015), concept drift detection 2019) and Tukey depth estimation 2020).
We focus on a family of lightweight and efficient methods called incremental quantile estimators, see Yazidi and Hammer (2017); Hammer et al. (2019) for recent reviews. The methods are based on performing small updates of the quantile estimate every time a new sample is received from the data stream and the methods can document state-of-the-art performance on quantile tracking. We suggest two new procedures that select and tune the update size, respectively, by the guidance of an estimate of the current tracking mean squared error (MSE). Naturally the MSE cannot be computed directly since the true quantile is unknown. The MSE is estimated by decomposing it in estimation variance and bias and novel and efficient procedures to estimate these quantities are presented. In particular we show that estimation bias can be associated with the portion of observations below the quantile estimates. The procedures are highly efficient with O(1) memory requirement and O(1) computational complexity per iteration.
The paper is organized as follows. Section 2 gives a short presentation of incremental quantile estimators, Section 3 presents how to estimate the current tracking MSE and how to use it to adapt the quantile estimates. Sections 4 and 5 evaluate the suggested procedure using synthetic and real-life data streams while Section 6 gives some closing remarks.

Incremental Quantile Estimators
Let X t ∼ f t (x) represent possible outcomes from a data stream at time t, x t a random sample and Q t,q the quantile associated with probability q, i.e P (X t ≤ Q t,q ) = F t (Q t,q ) = q.
Incremental quantile estimator update a quantile estimate every time a new observation is received. The algorithms are initiated with an estimate Q 0,q (λ) and further recursively updated where D 1 and D 2 are positive and can be deterministic or random. The estimator is intuitive in the sense that if the received observation is above (below) the current estimate, the estimate is increased (decreased). The functions are typically further constructed to ensure that the estimator converges to the underlying true quantile (Tierney;1983). A prominent example is the deterministic based multiplicative incremental quantile estimator (DUMIQE)  where D 1 (q, Q t,q (λ)) = q Q t,q (λ) an D 2 (q, Q t,q (λ)) = (1 − q) Q t,q (λ). Another example is the Frugal estimator (Ma et al.; where D 1 (q, Q t,q (λ)) = I(q < U ) and D 2 (q, Q t,q (λ)) = I(1 − q < U ) and U uniformly distributed on [0, 1] and I(·) the indicator function.
The tuning parameter λ determines update size. If the data stream distribution changes rapidly (slowly) with time, a high (small) should be used. Further, the step size should be adjusted to the scale of the data.

Adaptive Updating of Tuning Parameters
In this section we explain how to estimate the current quantile tracking mean squared error (MSE) and how to use this to efficiently track the true quantile. The MSE estimator uses the well-known bias and variance decomposition where µ(·) and σ 2 (·) refer to expectation and variance, respectively. For notational convenience, references to λ are avoided in this part.
The variance in Eq.
(2) is tracked using exponential smoothing, but with some careful considerations. For dynamically changing data streams, σ 2 ( Q t,q ) typically is very small relative to the dynamic changes of the data stream 6 making this a hard tracking problem. A common approach is to rely on the result that σ 2 ( Q t,q ) = µ( Q 2 t,q ) − µ( Q t,q ) 2 and separately track the two expectations, but this approach is prone to loss of precision and cannot guarantee that the estimate will be positive. Instead we suggest the following estimator which is a generalization of Knuth (2014) that used β = α. When α = β the estimator will be biased, but our experiments showed that the bias is small even for large difference between α and β. If a too small value of α is used in Eq. (3), the estimator will not be able the track the expectation resulting in serious overestimation of the variance in Eq. (4). By rather using a large value, say α = 0.5, combined with a small value of β to efficiently smooth the volatile local variance estimates worked well in our experiments. Naturally the Bias 2 term in Eq.
(2) cannot be computed since Q t,q is unknown. However, the following computation shows that it can be associated with the expected portion of observations below the quantile estimate which is possible to estimate: where Eq. (7) is derived using a first order Taylor approximation. H q (µ ( q t )) = (µ ( q t ) − q) 2 in Eq. (9) is tracked in a two stage procedure similar to the variance procedure above Finally G t,q is estimated by tracking an auxiliary quantile for a second probability q close to q using the same value of λ An estimate of the MSE is now obtained by substituting the estimates into (2) A useful property of the MSE estimator is that it does not use the values of a data stream directly making it robust to outliers if the quantile estimator is robust to outliers which is the case for most incremental estimators. We suggest two strategies to efficiently track the true quantile using Eq. (1). Oracle: Let λ 1 < λ 2 < · · · < λ L span all reasonable values for the tuning parameter and track Figure 1: An overview of the Oracle approach. In the bottom row, the gray dots show observations from the data stream. The black curves show tracking of the q = 0.7 quantile for different values of λ (λ 1 , λ 2 and λ L ). The information in the panels are used to estimate MSE for each value of λ using the procedure in Section 3. The MSE estimates are sent to the Oracle which selects the current quantile estimate.
Q t,q (λ l ), l = 1, . . . , L using Equations (1). For each Q t,q (λ l ) estimate the associated MSEs, MSE t (λ l ), l = 1, . . . , L, using (13). Let λ t = arg min λ l ,l∈1,...,L MSE t (λ l ) and let the current quantile estimate be given by Q t,q ( λ t ). Naturally, we may get large fluctuations in the values of λ used and one may limit updates only to neighbouring values (friction). We denote this approach the Oracle approach, since we can imagine an oracle that monitors the individual quantile tracking procedures and use the estimated MSEs to select the best current quantile estimate, without disturbing the quantile or MSE tracking procedures. The procedure is illustrated in Figure 1 where the objective is to track the q = 0.7 quantile of the data stream (gray dots). In the bottom left panel the quantile is tracked with a small update size, λ 1 , resulting in a high estimation bias according to Eq. (9) (the portion of observations below the quantile estimates is far below the target q = 0.7). Further in the bottom right panel the quantile is tracked with a large update size, λ L , resulting in high tracking variance according to Eq. (4). The estimator in the middle panel seem to have both fairly small bias and variance resulting in a small MSE. A potential challenge with the approach, is that with no knowledge about the data stream it may be difficult to select the range of λ's. A potential solution is that if the oracle selects the quantile estimate for a λ close to λ 1 (or λ L ) additional estimators with values of λ below λ 1 (or above λ L ) are included. The approach needs to store a total of 8L values and perform a total of 8L operations for each received sample from the data stream where L = 100 is a representative value. Thus for extremely massive data streams, this may pose computational challenge. An example is when incremental quantile estimators are used to track depth contours requiring to track thousands or even millions of quantile estimates 2020). HIL: A computationally less demanding approach is to only track the quantile for a high, intermediate and low (HIL) value of λ. Thus, the procedure only requires to store a total of 8 · 3 values and perform a total of 8 · 3 operations for each received data stream sample.
• Run three quantile tracking algorithms in parallel using tuning parameters λ 1ow = λ/a, λ intermediate = λ and λ high = aλ, a > 1 and track the MSE for each of them.
• Every M iterations, update λ: -If MSE is smallest for λ 1ow (or λ high ), reduce (or increase) the value of the tuning parameter for the three estimators by setting λ ← λ/a (or λ ← aλ). Restart the three quantile estimators initialized with the currently best quantile estimate, i.e. for λ 1ow (or λ high ).
-If MSE is smallest using λ intermediate , no updates are done.
In practice changes in data stream dynamics happen slowly or rarely (else it would be part of the original dynamics) and it makes sense to use a high value of M , say 1, 000. Further the values of the tuning parameters β, γ, κ and η should be chosen such that when an update of λ is performed, the estimate of the MSE has converged and at the same time most of the information since last update of λ is used. A simple rule of thumb is to set β = γ = κ = η = 1− M √ 0.01 which means that the weight of the M th term of the exponentially weighted sum is 0.01. This worked well in our experiments. Often data streams follow different periodic patterns, thus randomly selecting when to update λ can sometimes be useful.
The main challenge with the HIL approach is that the quantile estimates and associated MSEs are restarted after an update and must have converged before a new update is performed. This limits how rapidly and smoothly the procedure can adapt to changes in the data stream dynamics. Further since λ is updated rarely, a fairly large value of a, say 2, must be used limiting any fine tuning of λ.

Synthetic Experiments
DUMIQE was used for all the experiments in this section. Consider a normally distributed data stream where the expectation changes between slow and rapid dynamics with τ 1 = 500, τ 2 = 10 4 , µ = 8, b = 2, σ = 1, T = 10 4 . Figure 2 shows tracking of Bias 2 , Variance, MSE and values of λ using the Oracle approach with β = γ = κ = η = 1 − 1000 √ 0.01 = 0.005. Further, the following values of λ were used λ 1 = exp(−7), λ 2 = exp(−6.95), . . . , exp(0). Figure 3 shows the resulting tracking. Top left panel of Figure 2 shows that using a high value of λ, i.e. 2λopt, (solid gray curve) results in small bias since the quantile algorithm is able to efficiently track the dynamic changes of the data stream. On the other, as shown in the in the upper right panel, the estimation variance will be high since the step size is high. The MSE in the bottom left panel is further the sum of these to quantities and finding the best λ can be interpreted as a continuous battle between tracking bias and variance. The bottom right panel shows that the procedure quite rapidly finds values centered around the theoretically optimal λ. When the dynamics of the data stream changes after 10, 000 iterations, the procedure rapidly selects an estimator with a smaller value of λ however slightly higher then the theoretical optimal λ. It is important to recall that Eq. (13) is only an approximation of the true MSE that is further estimated. We see that both bias and variance are smaller after 10, 000 iterations which is as expected since slow dynamics are easier to track.
In Figure 3 we see that for the first iterations, the tracking is erratic until a reliable estimates of MSE is obtained.
We now analyze the performance of the suggested procedure in more detail. Consider again the data stream in (14) as well as a χ 2 distributed stream f t (x) = χ 2 ν + b sin 2π τ (n) n τ (n) = τ 1 I(mod n < T ) + τ 2 I(T ≤ mod n < 2T )   (14) and the black line is tracking of the q = 0.7 quantile with the Oracle approach.
using the same values τ 1 , τ 2 , T and b as above and ν = 6. χ 2 (ν) represents the χ 2 distribution with ν degrees of freedom. The χ 2 distributed stream is challenging since both the expectation and variance change with time. For the HIL approach, a = 1.5 and M = 1000 + U was used, where U was uniformly distributed on the interval [0, 1000], i.e. λ was updated on average every 1500 iteration. We tracked the q = 0.5, 0.7 and 0.9 quantiles and used q = 0.6, 0.6 and 0.8, respectively. To remove any Monte Carlo error, data streams were ran for a total of N = 10 7 iterations and observed tracking MSE were computed The results are shown in Figure 4. Let MSE theo τ 1 and MSE theo τ 2 represent the theoretically minimum tracking MSE for a data stream with constant fast dynamics (τ 1 ) and constant slow dynamics (τ 2 ), respectively. The data streams above, consist of an equal amount of fast and slow dynamics and thus the theoretically minimum tracking error will be 0.5 MSE theo τ 1 + 0.5 MSE theo τ 2 and is shown as gray solid lines in Figure 4. We further computed the theoretically minimum tracking MSE for the streams above using a constant value of λ and is shown as gray dashed lines. The results reveal that the Oracle approach is close to the theoretically optimal tracking error for all the cases. The HIL approach also performs well, but not as good as Oracle approach, which is as expected given the discussions at the end of Section 3. Optimal values of the tuning parameters in the procedures (β, γ, κ, η) are close to the suggested rule of thumb 1 − M √ 0.01 = 0.0031.

Real-Life Data Example -Tracking of Twitter Traffic
In this section we present a real-world example for benchmarking the procedures. Analysis of Twitter steaming data is an important application in research but also in practice. Efficient analysis of Twitter streams allows to for example perform detection of natural disasters or is used to predict election outcomes. For our experiment, let T t represent the time point when tweet number n is received from a stream of tweets. In this section we will track quantiles of the quantity R t = (T t − T t−1 ) −1 which can be interpreted as the frequency of received tweets. The quantity can for instance be used to real world event detection in the sense that if the number of received tweets is increasing, R t will increase (Atefeh and Khreich;2015;Hasan et al.;.
We consider a dataset consisting of the time points of every posted tweet by Norwegian users before and after the horrific terrorist attack July 22 2011 (Sollid et al.;2012). The time stamp for each tweet were given as whole second. To reconstruct a representation of the true time stamps a uniformly drawn value between zero and one second were added to each time stamp. The data stream is massive and we assumed no knowledge about suitable values for λ and thus used the HIL approach. We tracked quantiles of R t using the Frugal estimator (Ma et al.;. Since the update size of the Frugal estimator is independent of the scale of the data, the suggested procedure must adjust λ if the scale of the R t changes as well as if the dynamics changes. The value of λ was adjusted after every M = 1000 tweet received, which on average was about every 18 th minute in real-time. The value of λ was scaled by a = 1.5 and the tuning parameters were β, γ, κ and η was set to 0.005 in accordance with the rule of thumb in Section 3.
The results are shown in Figure 5. The attack was initiated by a bomb going of in Oslo July 22 at 3:25 p.m. local time and is marked on the x axis. The black curve in the upper panel shows the tracking of the q = 0.7 quantile of R t . The gray dots show the observed R t . The bottom panel shows the value of λ in every iteration. The Frugal estimator was initiated with λ = 1 and the procedure rapidly adjusted it to a more suitable value. During the terrorist attack at some point a bomb placed by the terrorist exploded which had an impact on the tweets.
For the period until the bomb went of it is a trend that λ is adjusted to lower values during night time which makes sense since both scale and dynamics of R t are small. When the bomb went off, R t changed rapidly. The value λ was rapidly increased to be able to efficiently track under the increased scale and dynamics. However, when the data stream stabilizes on a higher scale a few hours of the bomb went of the value of λ is slightly decreased. The value of λ is again increased when R t rapidly decreases (rapid dynamics) during the evening of July 22. More or less the same trend in λ is observed during July 23, a high value of λ is used in the morning to track the increased scale and dynamics, reduction of λ during daytime since the dynamics stabilized and again an increase in the value of λ in the evening of July 23. Overall the experiment demonstrates that the suggested procedure is able to adjust λ to obtain efficient quantile tracking in all parts of the complex data stream.

Closing Remarks
Surprisingly little attention has been given to automatic adjustment of the values of tuning parameters of incremental quantile estimators. In this paper we develop two simple procedures to achieve this. Both procedures are based on estimating the current tracking MSE and use this to efficiently track the true quantiles. The Oracle approach tracks the quantile and associated MSE for a wide range of values of the tuning parameter and in each iteration select the quantile estimate with minimal estimated MSE. The second approach only tracks the quantile for three values of the tuning parameter and repeatedly forget the estimate with the highest estimated MSE and add a quantile estimator for another value of the tuning parameter. Both methods are computationally and memory efficient since only a limited set of quantities needs to be computed and stored in each iteration.
The results show that the methods are highly efficient in adjusting the value of λ to achieve efficient tracking. The synthetic experiments showed that the resulting tracking error is close to the theoretical minimum. The Oracle performs the best, but at a higher computational cost. The real-life data example demonstrated that the procedure is able to handle complex and massive real-life data streams.