<commit_before>\documentclass[11pt]{amsart}
\usepackage{geometry}                % See geometry.pdf to learn the layout options. There are lots.
\geometry{letterpaper}                   % ... or a4paper or a5paper or ... 
%\geometry{landscape}                % Activate for for rotated page geometry
%\usepackage[parfill]{parskip}    % Activate to begin paragraphs with an empty line rather than an indent
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{epstopdf}
\usepackage[linesnumbered,vlined,algoruled]{algorithm2e}
\usepackage{url}
\usepackage{mathtools}
\DeclarePairedDelimiter{\ceil}{\lceil}{\rceil}
\DeclarePairedDelimiter{\floor}{\lfloor}{\rfloor}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator{\argmin}{arg\,min}

\title{Computing Extremely Accurate Quantiles Using $t$-Digests}
\author{<NAME>}
\address{Ted Dunning \\ MapR Technologies, Inc \\ San Jose, CA}
\email{<EMAIL>}
\date{}                                           % Activate to display a given date or no date
\begin{document}
\begin{abstract}
An on-line algorithm for computing approximations of rank-based statistics is presented that allows controllable accuracy.  Moreover, this new algorithm can be used to compute hybrid statistics such as trimmed means in additional to computing arbitrary quantiles.  An unusual property of the method is that it allows a quantile $q$ to be computed with an accuracy relative to $q (1-q)$ rather than with an absolute accuracy as with most methods.  This new algorithm is robust with respect to highly skewed distributions or highly ordered datasets and allows separately computed summaries to be combined with no loss in accuracy.

An open-source implementation of this algorithm is available from the author.
\end{abstract}
\maketitle
\section{Introduction}
Given a sequence of numbers, it is often desirable to compute rank-based statistics such as the median, 95-th percentile or trimmed means in an on-line fashion, keeping only a small data structure in memory.  Traditionally, such statistics were often computed by sorting all of the data and then either finding the quantile of interest by interpolation or by processing all samples within particular quantiles.  This sorting approach is infeasible for very large datasets which has led to interest in on-line approximate algorithms.

One early algorithm for computing on-line quantiles is described in \cite{Chen2000}.  In that work specific quantiles were computed by incrementing or decrementing an estimate by a value proportional to the simultaneously estimated probability density at the desired quantile.  This method is plagued by a circularity in that estimating density is only possible by estimating yet more quantiles.  Moreover, this work did not allow the computation of hybrid quantities such as trimmed means.

\cite{Greenwald-space-efficient-online-quantiles}Munro and Paterson provided an alternative algorithm to get an accurate estimate of any particular quantile.  This is done by keeping $s$ samples from the $N$ samples seen so far where $s << N$ by the time the entire data set has been seen.  If the data are presented in random order and if $s = \theta(N^{1/2})$, then Munro and Paterson's algorithm has a high probability being able to retain a set of samples that contains the median.  This algorithm can be adapted to find a number of pre-specified quantiles at the same time at proportional cost in memory.  The memory consumption of Munro-Paterson algorithm is excessive if precise results are desired.  Approximate results can be had with less memory, however.  A more subtle problem is that the implementation in Sawzall\cite{sawzall} and the Datafu library\cite{datafu} uses a number of buckets computed from the GCD of the desired quantiles.  This means that if you want to compute the $99$-th, $99.9$-th and $99.99$-th percentiles, a thousand buckets are required.

An alternative approach is described in \cite{qdigest}.  In this work, incoming values are assumed to be integers of fixed size. Such integers can trivially be arranged in a perfectly balanced binary tree where the leaves correspond to the integers and the interior nodes correspond to bit-wise prefixes. This tree forms the basis of the data structure known as a Q-digest.  The idea behind a Q-digest is that in the uncompressed case, counts for various values are assigned to leaves of the tree.  To compress this tree, sub-trees are collapsed and counts from the leaves are aggregated into a single node representing the sub-tree such that the maximum count for any collapsed sub-tree is less than a threshold that is a small fraction of the total number of integers seen so far.  Any quantile can be computed by traversing the tree in left prefix order, adding up counts until the desired fraction of the total is reached.  At that point, count for the last sub-tree traversed can be used to interpolate to the desired quantile within a small and controllable error.  The error is bounded because the count for each collapsed sub-tree is bounded.

The two problems with the Q-digest are that it depends on the tree structure being known ahead of time and that the error bounds do not necessarily apply if the algorithm is used in an on-line fashion.  Adapting the Q-digest to use an balanced tree over arbitrary elements is difficult.  This difficulty arises because rebalancing the tree involves sub-tree rotations and these rotations may require reapportionment of previously collapsed counts in complex ways.  This reapportionment could have substantial effects on the accuracy of the algorithm and in any case make the implementation much more complex because the concerns of counting cannot be separated from the concerns of maintaining a balanced tree.  Another problem with Q-digests is that if they are subjected to compression during building, it isn't entirely clear how to handle compressed counts that move high above the leaves, but which eventually have non-trivial counts at a lower level in the tree.  The proof of correctness for Q-digests ignores this complication by only considering the case where counts are compressed after being collected on the leaves.  It would be desirable to have error bounds that apply to a completely on-line data structure.  These limitations do not apply in the original formulation of the Q-digest as a compression algorithm for quantile information, but current trends towards the use of on-line algorithms make these limitations awkward.

The work described here shows how the fundamental idea of a Q-digest can be easily extended to values in $\mathbb R$ without the complexity of apportioning counts during tree rebalancing.  Indeed, this new data structure eliminates the idea of a tree for storing the original samples, maintaining only the concept of collapsing groups of observations in a way that preserves accuracy.  This new algorithm, known as $t$-digest, has well-defined and easily proven error bounds and allows parallel on-line operation.  A particularly novel aspect of the variant of $t$-digest described here is that accuracy for estimating the $q$ quantile is relative to $q(1-q)$.  This is in contrast to earlier algorithms which had errors independent of $q$.  The relative error bound of the $t$-digest is convenient when computing quantiles for very small $q$ or for $q$ near $1$.  As with the Q-digest algorithm, the accuracy/size trade-off can be controlled by setting a single compression parameter.

The accuracy bounds are tight regardless of data ordering for the non-parallel on-line case.  In the case of parallel execution, the error bound is somewhere between the constant error bound and the relative error bound except in the case of highly ordered input data.  For randomly ordered data, parallel execution has expected error bounded by the same bound as for sequential execution.

\section{The Algorithm}
The basic outline of the algorithm for constructing a $t$-digest is quite simple.  An initially empty ordered list of centroids, $C = [ c_1 \ldots c_m ]$ is kept.  Each centroid consists of a mean and a count.  To add a new value $x_n$ with a weight $w_n$, the set of centroids is found that have minimum distance to $x_n$.  This set is reduced by retaining only centroids with a count less than $4\delta q (1-q) n $ where $\delta$ controls the accuracy of quantile estimates and $q$ is the estimated quantile for the mean of the centroid.  If more than one centroid remains, one is selected at random.  If a centroid is found, then $(x_n,w_n)$ is added to that centroid.  It may happen that the weight $w_n$ is larger than can be added to the selected centroid.  If so, as much weight as possible is allocated to the selected centroid and the selection is repeated with the remaining weight.  If no satisfactory centroid is found or if there is additional weight to be added after all centroids with minimum distance are considered, then $x_n$ is used to form a new centroid and the next point is considered.  This procedure is shown more formally as Algorithm \ref{alg:full}.
 \begin{algorithm}[tb]
\SetKw{KwTo}{in}\SetKwFor{For}{for}{\string:}{}
\SetKwIF{If}{ElseIf}{Else}{if}{:}{elif}{else:}{}
\SetKwFor{While}{while}{:}{fintq}
\SetKwFor{For}{for}{\string:}{}
 \label{alg:full}
\SetNoFillComment
\KwIn{Ordered set of weighted centroids $C = \lbrace \rbrace$, sequence of real-valued, weighted points $X = \lbrace (x_1, w_1),\ldots (x_N, w_N)\rbrace$, and accuracy tolerance $\delta$}
\KwOut{final set $C=[c_1 \ldots c_m]$ of weighted centroids} 
\For{$(x_n, w_n) \in X$} {
  $z = \min | c_i.\mathtt{mean} - x |$\;
  $S = \lbrace c_i  :  |c_i.\mathtt{mean} - x| = z \wedge c_i.\mathtt{count} + 1 \le 4n\delta q(c_i) (1-q(c_i)) \rbrace $\;
   \While {$S \ne \lbrace \rbrace\wedge w_n > 0$} {
       Sample $c_j \sim \mathrm{Uniform}( S)$\;
       $\Delta w = \min(4n\delta q(c_i) (1-q(c_i))-c_j.\mathtt{count}, w_n)$\;
       $c_j.\mathtt{count} \gets c_j.\mathtt{count} + \Delta w$\;
       $c_j.\mathtt{mean} \gets c_j.\mathtt{mean} + \Delta w (x_n - c_j.\mathtt{mean})/c_j.\mathtt{count}$\;
       $w_n \gets w_n - \Delta w$\;
       $S \gets S - c_j$
     } 
     \If {$w_n > 0$} {
        $C \gets C + (x_n, w_n)$\;
      }
      \If {$|C| > K/\delta$} {
         $C \gets \mathtt{cluster}(\mathtt{permute}( C )) $\;
       }
} 
$C \gets \mathtt{cluster}(\mathtt{permute}( C )) $\;
\Return $ C $\\
\caption{Construction of a $t$-Digest}
\end{algorithm}

In this algorithm, a centroid object contains a mean and and a count.  Updating such an object with a new data-point $(x,w)$ is done using Welford's method \cite{wiki:welford, knuth2welford, welford62}.

The number of points assigned to each centroid is limited to $\max(1, \floor{4N\delta q(1-q})$ where $q$ is the quantile for the approximate mean of the centroid and $\delta$ is a parameter chosen to limit the number of points that can be assigned to a centroid.  Typically, this compression factor is described in terms of its inverse, $1/\delta$ in order to stay compatible with the conventions used in the Q-Digest.  The algorithm approximates $q$ for centroid $c_i$ by summing the weights for all of the centroids ordered before $c_i$:

\[q(c_i)=\frac{c_i.\mathtt{count}/2 + \sum_{j<i}c_j.\mathtt{count}} {\sum_j c_j.\mathtt{count}}\]

In order to compute this sum quickly, the centroids can be stored in a data structure such as a balanced binary tree that keeps sums of each sub-tree.  For centroids with identical means, order of creation is used as a tie-breaker to allow an unambiguous ordering.  Figure \ref{fig:gamma-sizes} shows actual centroid weights from multiple runs of this algorithm and for multiple distributions plotted against the ideal bound.
<commit_msg>update t-digest algo in tex<commit_after>\documentclass[11pt]{amsart}
\usepackage{geometry}                % See geometry.pdf to learn the layout options. There are lots.
\geometry{letterpaper}                   % ... or a4paper or a5paper or ... 
%\geometry{landscape}                % Activate for for rotated page geometry
%\usepackage[parfill]{parskip}    % Activate to begin paragraphs with an empty line rather than an indent
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{epstopdf}
\usepackage[linesnumbered,vlined,algoruled]{algorithm2e}
\usepackage{url}
\usepackage{mathtools}
\DeclarePairedDelimiter{\ceil}{\lceil}{\rceil}
\DeclarePairedDelimiter{\floor}{\lfloor}{\rfloor}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator{\argmin}{arg\,min}

\title{Computing Extremely Accurate Quantiles Using $t$-Digests}
\author{<NAME>}
\address{<NAME> \\ MapR Technologies, Inc \\ San Jose, CA}
\email{<EMAIL>}
\date{}                                           % Activate to display a given date or no date
\begin{document}
\begin{abstract}
An on-line algorithm for computing approximations of rank-based statistics is presented that allows controllable accuracy.  Moreover, this new algorithm can be used to compute hybrid statistics such as trimmed means in additional to computing arbitrary quantiles.  An unusual property of the method is that it allows a quantile $q$ to be computed with an accuracy relative to $q (1-q)$ rather than with an absolute accuracy as with most methods.  This new algorithm is robust with respect to highly skewed distributions or highly ordered datasets and allows separately computed summaries to be combined with no loss in accuracy.

An open-source implementation of this algorithm is available from the author.
\end{abstract}
\maketitle
\section{Introduction}
Given a sequence of numbers, it is often desirable to compute rank-based statistics such as the median, 95-th percentile or trimmed means in an on-line fashion, keeping only a small data structure in memory.  Traditionally, such statistics were often computed by sorting all of the data and then either finding the quantile of interest by interpolation or by processing all samples within particular quantiles.  This sorting approach is infeasible for very large datasets which has led to interest in on-line approximate algorithms.

One early algorithm for computing on-line quantiles is described in \cite{Chen2000}.  In that work specific quantiles were computed by incrementing or decrementing an estimate by a value proportional to the simultaneously estimated probability density at the desired quantile.  This method is plagued by a circularity in that estimating density is only possible by estimating yet more quantiles.  Moreover, this work did not allow the computation of hybrid quantities such as trimmed means.

\cite{Greenwald-space-efficient-online-quantiles}Munro and Paterson provided an alternative algorithm to get an accurate estimate of any particular quantile.  This is done by keeping $s$ samples from the $N$ samples seen so far where $s << N$ by the time the entire data set has been seen.  If the data are presented in random order and if $s = \theta(N^{1/2})$, then Munro and Paterson's algorithm has a high probability being able to retain a set of samples that contains the median.  This algorithm can be adapted to find a number of pre-specified quantiles at the same time at proportional cost in memory.  The memory consumption of Munro-Paterson algorithm is excessive if precise results are desired.  Approximate results can be had with less memory, however.  A more subtle problem is that the implementation in Sawzall\cite{sawzall} and the Datafu library\cite{datafu} uses a number of buckets computed from the GCD of the desired quantiles.  This means that if you want to compute the $99$-th, $99.9$-th and $99.99$-th percentiles, a thousand buckets are required.

An alternative approach is described in \cite{qdigest}.  In this work, incoming values are assumed to be integers of fixed size. Such integers can trivially be arranged in a perfectly balanced binary tree where the leaves correspond to the integers and the interior nodes correspond to bit-wise prefixes. This tree forms the basis of the data structure known as a Q-digest.  The idea behind a Q-digest is that in the uncompressed case, counts for various values are assigned to leaves of the tree.  To compress this tree, sub-trees are collapsed and counts from the leaves are aggregated into a single node representing the sub-tree such that the maximum count for any collapsed sub-tree is less than a threshold that is a small fraction of the total number of integers seen so far.  Any quantile can be computed by traversing the tree in left prefix order, adding up counts until the desired fraction of the total is reached.  At that point, count for the last sub-tree traversed can be used to interpolate to the desired quantile within a small and controllable error.  The error is bounded because the count for each collapsed sub-tree is bounded.

The two problems with the Q-digest are that it depends on the tree structure being known ahead of time and that the error bounds do not necessarily apply if the algorithm is used in an on-line fashion.  Adapting the Q-digest to use an balanced tree over arbitrary elements is difficult.  This difficulty arises because rebalancing the tree involves sub-tree rotations and these rotations may require reapportionment of previously collapsed counts in complex ways.  This reapportionment could have substantial effects on the accuracy of the algorithm and in any case make the implementation much more complex because the concerns of counting cannot be separated from the concerns of maintaining a balanced tree.  Another problem with Q-digests is that if they are subjected to compression during building, it isn't entirely clear how to handle compressed counts that move high above the leaves, but which eventually have non-trivial counts at a lower level in the tree.  The proof of correctness for Q-digests ignores this complication by only considering the case where counts are compressed after being collected on the leaves.  It would be desirable to have error bounds that apply to a completely on-line data structure.  These limitations do not apply in the original formulation of the Q-digest as a compression algorithm for quantile information, but current trends towards the use of on-line algorithms make these limitations awkward.

The work described here shows how the fundamental idea of a Q-digest can be easily extended to values in $\mathbb R$ without the complexity of apportioning counts during tree rebalancing.  Indeed, this new data structure eliminates the idea of a tree for storing the original samples, maintaining only the concept of collapsing groups of observations in a way that preserves accuracy.  This new algorithm, known as $t$-digest, has well-defined and easily proven error bounds and allows parallel on-line operation.  A particularly novel aspect of the variant of $t$-digest described here is that accuracy for estimating the $q$ quantile is relative to $q(1-q)$.  This is in contrast to earlier algorithms which had errors independent of $q$.  The relative error bound of the $t$-digest is convenient when computing quantiles for very small $q$ or for $q$ near $1$.  As with the Q-digest algorithm, the accuracy/size trade-off can be controlled by setting a single compression parameter.

The accuracy bounds are tight regardless of data ordering for the non-parallel on-line case.  In the case of parallel execution, the error bound is somewhere between the constant error bound and the relative error bound except in the case of highly ordered input data.  For randomly ordered data, parallel execution has expected error bounded by the same bound as for sequential execution.

\section{The Algorithm}
The basic outline of the algorithm for constructing a $t$-digest is quite simple.  An initially empty ordered list of centroids, $C = [ c_1 \ldots c_m ]$ is kept.  Each centroid consists of a mean and a count.  To add a new value $x_n$ with a weight $w_n$, the set of centroids is found that have minimum distance to $x_n$.  This set is reduced by retaining only centroids with a count less than $4\delta q (1-q) n $ where $\delta$ controls the accuracy of quantile estimates and $q$ is the estimated quantile for the mean of the centroid.  If more than one centroid remains, one is selected at random.  If a centroid is found, then $(x_n,w_n)$ is added to that centroid.  It may happen that the weight $w_n$ is larger than can be added to the selected centroid.  If so, as much weight as possible is allocated to the selected centroid and the selection is repeated with the remaining weight.  If no satisfactory centroid is found or if there is additional weight to be added after all centroids with minimum distance are considered, then $x_n$ is used to form a new centroid and the next point is considered.  This procedure is shown more formally as Algorithm \ref{alg:full}.
 \begin{algorithm}[tb]
\SetKw{KwTo}{in}\SetKwFor{For}{for}{\string:}{}
\SetKwIF{If}{ElseIf}{Else}{if}{:}{elif}{else:}{}
\SetKwFor{While}{while}{:}{fintq}
\SetKwFor{For}{for}{\string:}{}
 \label{alg:full}
\SetNoFillComment
\KwIn{Ordered set of weighted centroids $C = \lbrace \rbrace$, sequence of real-valued, weighted points $X = \lbrace (x_1, w_1),\ldots (x_N, w_N)\rbrace$, and accuracy tolerance $\delta$}
\KwOut{final set $C=[c_1 \ldots c_m]$ of weighted centroids} 
\For{$(x_n, w_n) \in X$} {
  $z = \min | c_i.\mathtt{mean} - x |$\;
  $S = \lbrace c_i  :  |c_i.\mathtt{mean} - x| = z \wedge c_i.\mathtt{count} + 1 \le 4n\delta q(c_i) (1-q(c_i)) \rbrace $\;
   \While {$S \ne \lbrace \rbrace\wedge w_n > 0$} {
       Sample $c_j \sim \mathrm{Uniform}( S)$\;
       $\Delta w = \min(4n\delta q(c_j) (1-q(c_j))-c_j.\mathtt{count}, w_n)$\;
       $c_j.\mathtt{count} \gets c_j.\mathtt{count} + \Delta w$\;
       $c_j.\mathtt{mean} \gets c_j.\mathtt{mean} + \Delta w (x_n - c_j.\mathtt{mean})/c_j.\mathtt{count}$\;
       $w_n \gets w_n - \Delta w$\;
       $S \gets S - c_j$
     } 
     \If {$w_n > 0$} {
        $C \gets C + (x_n, w_n)$\;
      }
      \If {$|C| > K/\delta$} {
         $C \gets \mathtt{cluster}(\mathtt{permute}( C )) $\;
       }
} 
$C \gets \mathtt{cluster}(\mathtt{permute}( C )) $\;
\Return $ C $\\
\caption{Construction of a $t$-Digest}
\end{algorithm}

In this algorithm, a centroid object contains a mean and and a count.  Updating such an object with a new data-point $(x,w)$ is done using Welford's method \cite{wiki:welford, knuth2welford, welford62}.

The number of points assigned to each centroid is limited to $\max(1, \floor{4N\delta q(1-q})$ where $q$ is the quantile for the approximate mean of the centroid and $\delta$ is a parameter chosen to limit the number of points that can be assigned to a centroid.  Typically, this compression factor is described in terms of its inverse, $1/\delta$ in order to stay compatible with the conventions used in the Q-Digest.  The algorithm approximates $q$ for centroid $c_i$ by summing the weights for all of the centroids ordered before $c_i$:

\[q(c_i)=\frac{c_i.\mathtt{count}/2 + \sum_{j<i}c_j.\mathtt{count}} {\sum_j c_j.\mathtt{count}}\]

In order to compute this sum quickly, the centroids can be stored in a data structure such as a balanced binary tree that keeps sums of each sub-tree.  For centroids with identical means, order of creation is used as a tie-breaker to allow an unambiguous ordering.  Figure \ref{fig:gamma-sizes} shows actual centroid weights from multiple runs of this algorithm and for multiple distributions plotted against the ideal bound.

{'fraction_non_alphanumeric': 0.07059965396463687, 'fraction_numerical': 0.004135544583702578, 'mean_word_length': 4.168593238822247, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5280123', 'n_tokens_mistral': 6395, 'n_tokens_neox': 6073, 'n_words': 3335}