Fast post-hoc method for updating moments of large datasets

TL;DR

This paper introduces a fast, memory-efficient post-hoc method for updating statistical moments of large datasets by only using previous moments and new data, significantly reducing storage and computation time.

Contribution

The proposed method allows updating dataset moments without recalculating from scratch, improving efficiency for large datasets and low-order moments.

Findings

Reduces data storage needs for moment updates

Speeds up computation for large datasets

Effective for low-order moments (n ≤ 10)

Abstract

Moments of large datasets utilise the mean of the dataset; consequently, updating the dataset traditionally requires one to update the mean, which then requires one to recalculate the moment. This means that metrics such as the standard deviation, $R^2$ correlation, and other statistics have to be `refreshed' for dataset updates, requiring large data storage and taking long times to process. Here, a method is shown for updating moments that only requires the previous moments (which are computationally cheaper to store), and the new data to be appended. This leads to a dramatic decrease in data storage requirements, and significant computational speed-up for large datasets or low-order moments (n $\lesssim$ 10).

Figures

Equations

$$M_n[\{x_i\}]\equiv \frac{1}{Z}\sum _{i=1}^{N}f(x_i)(x_i-\bar{x}{)}^n$$

$$Z\equiv \sum _{i=1}^{N}f(x_i)$$

$$\bar{x}\equiv \frac{1}{Z}\sum _{i=1}^{N}f(x_i)x_i$$

$$\frac{T}{u}=\mathcal{O}(N)$$

$$S_{\mathbf{X}}=\sum _{i=1}^N[\text{size}(x_i)+\text{size}(f(x_i))]$$

$$S_n=\text{size}(M_n)$$

$$\begin{array}[]{r c l}p\circ q\in\mathbb{G}&;&p\times q\in\mathbb{G}\end{array}$$

$$p\circ q\equiv (p+q)\text{mod}10$$

$$p\circ q\equiv [\Re (p)+\Re (q)]+i[\Im (p)+\Im (q)]$$

$$\mathbf{A}\circ \mathbf{B}\equiv \{a_i\circ b_i\}$$

$$\mathbf{A}\times \mathbf{B}\equiv \{a_i\times b_i\}$$

$$\mathbf{A}^{n}\equiv\big{\{}(a_{i})^{n}\big{\}}$$

$$\begin{array}[]{r l}\circ&:\mathbb{D}\times\mathbb{D}\to\mathbb{D}\\ \times&:\mathbb{G}\times\mathbb{D}\to\mathbb{D}\end{array}$$

$$\mathbf{X}\subset {\mathbf{X}}^{\prime }$$

$$\Delta \equiv N^{\prime }-N$$

$$\sum _{i=N+1}^{N^{\prime }}\equiv \sum ^{\Delta }$$

$$M_n^{\prime }=\frac{Z}{Z^{\prime }}\left[(\bar{x}-{\bar{x}}^{\prime }{)}^n+\sum _{k=0}^{n-2}(\genfrac{}{}{0ex}{}{n}{k})M_{n-k}(\bar{x}-{\bar{x}}^{\prime }{)}^k\right]+\frac{1}{Z^{\prime }}\sum ^{\Delta }\left[f(x_i)(x_i-{\bar{x}}^{\prime }{)}^n\right]$$

$$Z^{\prime }=Z+\sum ^{\Delta }f(x_i)$$

$${\bar{x}}^{\prime }=\frac{Z}{Z^{\prime }}\bar{x}+\frac{1}{Z^{\prime }}\sum ^{\Delta }f(x_i)x_i$$

$$S_{\Delta }=\text{size}(M_n)\cdot (n-1)+\sum ^{\Delta }[\text{size}(f(x_i))+\text{size}(x_i)]$$

$$\frac{T}{u}=\mathcal{O}[(n-1)\cdot \Delta ]$$

$$n-1\lesssim \frac{N^{\prime }}{\Delta }$$

$$n_{\text{thr.}}\sim \frac{N^{\prime }}{\Delta }+1$$

$$\begin{array}[]{r l}M_{n}^{\prime}&=\dfrac{Z}{Z^{\prime}}\bigg{[}(\bar{x}-\bar{x}^{\prime})^{n}+\displaystyle\sum\limits_{k=0}^{n-2}\binom{n}{k}M_{n-k}(\bar{x}-\bar{x}^{\prime})^{k}\\ &\hskip 50.0pt+\displaystyle\sum\limits_{k=n+1}^{\infty}\binom{n}{k}M_{n-k}(\bar{x}-\bar{x}^{\prime})^{k}\bigg{]}+\dfrac{1}{Z^{\prime}}\displaystyle\sum\limits^{\Delta}\left[f(x_{i})(x_{i}-\bar{x}^{\prime})^{n}\right]\end{array}$$

$$|A|=\left\{\sqrt{a_i\times a_i}\right\}$$

$$\begin{array}[]{r l}M_{n}^{\prime}&\approx\dfrac{Z}{Z^{\prime}}\bigg{[}(\bar{x}-\bar{x}^{\prime})^{n}+\displaystyle\sum\limits_{k=0}^{n-2}\binom{n}{k}M_{n-k}(\bar{x}-\bar{x}^{\prime})^{k}\\ &\hskip 50.0pt+\displaystyle\sum\limits_{k=n+1}^{n^{*}}\binom{n}{k}M_{n-k}(\bar{x}-\bar{x}^{\prime})^{k}\bigg{]}+\dfrac{1}{Z^{\prime}}\displaystyle\sum\limits^{\Delta}\left[f(x_{i})(x_{i}-\bar{x}^{\prime})^{n}\right]\end{array}$$

$$S_{\Delta }=\text{size}(M_n)\cdot (n^{\ast }-1)+\sum ^{\Delta }[\text{size}(f(x_i))+\text{size}(x_i)]$$

$$\frac{T}{u}=\mathcal{O}[(n^{\ast }-1)\cdot \Delta ]$$

$$W=\frac{1}{Z}\sum _{i=0}^{N}f(x_i)g(x_i)$$

$$W=\sum _{n=0}^{\infty }{c}_n{M}_n$$

Full text

1 Introduction

This year, the Oak Ridge National Laboratory machine SUMMIT[1] exceeded 1.0 exaops ($10^{18}$ operations per second) computational speed, bringing computing into the exascale era. In this era, scientific simulations will produce larger datasets than before, requiring large amounts of memory and data storage. In the plasma physics community, there has been much recent chatter about the strength of exascale computing and the associated data complications [2, 3]. Furthermore, the recent increase in popularity of machine learning and deep learning across a variety of different scientific sectors [4, 5, 6] has lead to a growing requirement to produce reduced metrics from datasets in an increasingly efficient fashion.

Computing moments with an updated dataset proves problematic for two reasons: the computational time taken to compute moments can scale with the size of the total data, and one has to retain a large amount of data in storage to compute new moments.

Here, we give an algorithm for computing moments of datasets which is instead computed in time proportional to the size of the data to be appended, and requires only a finite and typically small number of moments to be placed in data storage.

1.1 Definitions

We start by stating that we have a dataset $\mathbf{X}$ containing data objects $\{x_{i}\}$. Let us define the $n^{\textrm{th}}$ moment $M_{n}[\mathbf{X}]$ as the following:

[TABLE]

where $\mathbf{X}=\{x_{i}\}$, and $Z$ is the normalising function, defined here as the sum over $i$ of the weight function $f(x_{i})$:

[TABLE]

Accordingly, $M_{0}$ is equal to unity. $\bar{x}$ is the weighted mean of the dataset $\mathbf{X}$, such that $M_{1}$ is identically zero:

[TABLE]

1.2 Difficulties with calculating moments

Firstly, if one computes the original moment (1), this is completed in $\mathcal{O}(N)$ time. This is problematic, because if the dataset is updated $u$ times, then $\mathcal{O}(N)$ time is required on each update:

[TABLE]

where $T$ is the total time taken. The time taken scales linearly with the number of updates, but also with the full dataset size. With data warehouses, this can lead to significant latency issues; in the limit that the increase in dataset size tends to zero, $T/u$ is still finite and non-zero.

This means that each update requires a full ‘refresh’ which discourages small updates, but also each update requires more computational time than any update before it.

Secondly, if one has a dataset which has moments calculated beforehand, the original information is lost; the mapping $\mathbf{X}\to M_{n}$ is many-to-one, and therefore the mapping is not invertible. This means that one requires an amount of available data storage $S_{\mathbf{X}}$ given by:

[TABLE]

where $\textrm{size}(x_{i})$ is the size of the data $x_{i}$ and $\textrm{size}(f_{i})$ is the size of the data $f(x_{i})$. $S_{\mathbf{X}}$ increases with each update, and means that we are required to use larger and large amounts of data to calculate these moments.

However, in principle, the moments themselves take up a much smaller amount of data storage $S_{n}$:

[TABLE]

where $\textrm{size}(M_{n})$ is the size of the data required to store each moment. For comparison, if $x_{i}$, $f(x_{i})$ and $M_{n}$ are datatypes of the same size, then $S_{\mathbf{X}}=2NS_{n}$.

1.3 Dataset operations

We begin by defining operations on datasets. To do so, we will take the concept of datasets as a ‘vector’ and generalise the concept of vector space over a field (the field being the data objects, and the vector being the dataset); this will allow us to perform algebra using entire data structures, with a great deal of freedom in the available operations between data.

First, we define a dataset $\mathbf{A}$ as a set of data objects $\{a_{i}\}$. The dataset is categorized by the type of data objects that form the set, and the size of the dataset. The size of the dataset (cardinality) is referred to as $N$.

It becomes useful to consider the data objects themselves as belonging to a field under addition and multiplication, which we will call the ‘datatype field’. That is to say, for data objects $p$ and $q$ belonging to a set $\mathbb{G}$ corresponding to a given ‘datatype’:

[TABLE]

The addition and multiplication operators ($\circ$ and $\times$) can be defined differently for different datatypes. For example, if we desire that the datatype is a set of cyclic numbers no bigger than 10, one might find:

[TABLE]

whereas if the datatype were a set of a complex numbers:

[TABLE]

Next, we define the data addition operator between two datasets $\mathbf{A}=\{a_{i}\}$ and $\mathbf{B}=\{b_{i}\}$ by the following:

[TABLE]

Formally, the set of datasets $\mathbb{D}\equiv\{\mathbf{A},\mathbf{B},\dots\}$ and the data addition operator ‘$+$’ form an Abelian group with identity element $\{0\}$. There is no formal requirement for each member of $\mathbf{A}$ to have the same datatype.

The data multiplication operator between two datasets is defined by the following:

[TABLE]

The set of datasets $\mathbb{D}$ and the data multiplication operator ‘$\times$’ form an Abelian group with identity element $\{1\}$. As both groups are Abelian, one can now define a mathematical field formed by $\mathbb{D}$, data addition, and data multiplication. Exponentiation is simply an extension of multiplication:

[TABLE]

It becomes useful to define a vector space over the field $\mathbb{G}$ which the data objects belong to, using the field $\mathbb{D}$ as elements of the vector space. Accordingly, for addition and multiplication:

[TABLE]

For example, if $\mathbf{X}$ is a set of 10 floating point numbers (floats), then the corresponding data structure is all possible sets of 10 floats. Operations between $\mathbf{X}$ and another set of 10 floats are now defined, and each member of $\mathbf{X}$ is a float.

Using a more complicated example, one can examine $\mathbf{X}$ as a set of 10 lists, each containing 3 floats and 1 modulo number. In this case, $\mathbb{D}$ is all possible combinations of 10 lists each with 3 floats and 1 modulo number, and each member of $\mathbf{X}$ is a list of 3 floats and 1 modulo.

2 Primed moments

Suppose that we define the ‘updated’ dataset containing new values as $\mathbf{X}^{\prime}=\{x_{i}^{\prime}\}$ such that:

[TABLE]

$\mathbf{X}^{\prime}$ is now said to have a cardinality of $N^{\prime}$. The difference in cardinality is worth noting, and is herein referred to as $\Delta$:

[TABLE]

We will also introduce the following shorthand:

[TABLE]

2.1 Integer $n$

For integer $n$, the primed moment can be represented in the following form (see Appendix for proof):

[TABLE]

where $\binom{n}{k}$ are binomial coefficients. $Z^{\prime}$ can be determined using only a simple update:

[TABLE]

which is computed in $\mathcal{O}(\Delta)$ time. Similarly, $\bar{x}^{\prime}$ can be computed in $\mathcal{O}(\Delta)$ time:

[TABLE]

The only points from $\mathbf{X}^{\prime}$ that are required are the points not in $\mathbf{X}$; viz. the only data required to update the moment is the new data. If all moments use the same amount of memory, then the data storage required is only:

[TABLE]

This reduced data requirement leads to a computational speed up as well. The computational time taken per update is given by:

[TABLE]

One quickly show that by comparing to (4) for speed up, one requires at worst:

[TABLE]

This yields a threshold moment value:

[TABLE]

For this value of $n$ or greater, (14) produces the primed moment slower than one could achieve using the original dataset; the speed-up is greatest for small updates to the dataset. However, in all cases, this method uses significantly less data storage.

2.2 Non-integer $n$

For non-integer values of $n$, one can use the fractional continuation of the binomial expansion:

[TABLE]

Although it appears that the infinite sum is singular for $\bar{x}^{\prime}=\bar{x}$, we have already assumed a priori that $\bar{x}^{\prime}\neq\bar{x}$. As $\mathbf{X}$ exists in a vector space, it has a well defined norm. We expect that the infinite sum is Abel convergent for $|x_{i}-\bar{x}|>1$, where $|A|$ denotes the Euclidean norm of the dataset $A$ on the vector space established in Section 1.3:

[TABLE]

This condition is therefore related to the standard deviation of the original dataset. In short, we posit that for a set of data with a low standard deviation ($\sqrt{M_{2}}$), the infinite sum in (21) should converge. Therefore, one can truncate the sum so as to keep enough terms to yield the answer to a desired level of numerical accuracy:

[TABLE]

where $n^{*}$ is the cutoff moment. Therefore, to numerical accuracy, one requires for data storage:

[TABLE]

and the computational time taken per update is given by:

[TABLE]

3 Weighted metrics

Suppose that we define the metric $W$ as:

[TABLE]

where $g(x_{i})$ is an arbitrary function. Then, one can Taylor expand the function $g(x_{i})$ of the dataset around the mean value:

[TABLE]

where $c_{n}$ are Taylor coefficients.

3.1 Updating metrics

If one updates the function $g\to g^{\prime}$, and the dataset $\mathbf{X}\to\mathbf{X}^{\prime}$ then:

[TABLE]

where $c_{n}^{\prime}$ denotes the new Taylor coefficients. This leads to a very powerful tool; provided that the moments of the original dataset are known it is possible for us to generate any arbitrary metric.

It is worth noting that if one rearranges for the change between $W^{\prime}$ and $W$:

[TABLE]

In principle, if one had $g^{\prime}(x_{i})\equiv g(x_{i})$, then this method is not particularly useful; we are essentially Taylor expanding at two different points (because the mean shifts), making our life much harder. This is not immediately intuitive; our expansion for $g$ employed here is always around the mean point of the dataset, rather than a free choice.

However, if $W$ is unknown and $\{M_{n}\}$ are known, then (27) is useful.

One can swap the order of summation on the first term. If one does so:

[TABLE]

But as the binomial coefficients are monotonically increasing (at roughly $k^{n}$), the inside sum is not necessarily convergent. To guarantee Abel convergence, as $M_{n}$ should be monotonically increasing also, $c_{k+n}$ must satisfy:

[TABLE]

This requirement means that we can only update $W$ via this method if the function $g(x)$ has a convergent Taylor series. In such a case, one can use a cutoff for the first sum over $n$ in (27). Then:

[TABLE]

The first term is iteratively calculated in $\mathcal{O}(n^{*}!)$ time (such that one reaches numerical accuracy), while the second term is calculated in $\mathcal{O}(\Delta)$ time.

4 Demonstration

Here, moments of datasets are compared using Python. Python forms a nice framework for adding data structures as one can freely redefine the addition and multiplication operators on objects of an arbitrary class. This makes the formalism from Section 1.3 easy to code, allowing us to simply add and multiple entire datasets as if they were just scalars.

4.1 Single precision floating point

For single precision floating points with boundaries at 0 and 1:

[TABLE]

Two sets of floats, a 16 float set and a 256 float set were used to compute the original integer moments. Then, extra points were added to each set, and the new integer moments were calculated. For each speed test, a new, random set of floats were initialised using the numpy.random subpackage.

The scalar function $f(x_{i})$ is set to be a random map, such that each piece of data $x_{i}$ points to a random float. The random floats are reset for each speed test, just as the dataset $\mathbf{X}$ is.

The speed-up was calculated by measuring the time taken to calculate the moment using the full dataset $\mathbf{X}^{\prime}$, and by using the update method given by (14). For each value of $N^{\prime}-N$ and moment, the timeit package was used to measure 100 runs of the algorithm.

As expected, the speed-up is less than 1 for high moments when the number of data points added ($N^{\prime}-N$) is large with respect to the size of the original dataset ($N$). One finds that the maximum speed up (for the lowest order moment, $n=2$) is close to the theoretical maximum of a factor 16 speed up between the $N=16$ and the $N=128$ set, as shown in Figure 5.

4.2 List of single floats

We now examine a list of single floats. Using boundaries at 0 and 1 for each float, and setting the list length to 4:

[TABLE]

The minimum update (in terms of $N\to N+1$) now adds 4 floats to the dataset. In principle, this could be reformulated using the same analysis as in Section 4.1 but updating the dataset in groups of 4. As such, we expect that the speed-up will look very similar, but will be affected by the 4 fold global increase in datasize.

Again, we use numpy.random to generate a random set of floats, and $f(x_{i})$ is a random scalar map. As shown in Figure 6, the speed-up for 4 sets of 16 floats is fairly similar to that observed with a single set of 16 floats. However, when we increase to 4 sets of 256 floats, the speed up is much greater than the case where we examine a single set of 256 floats (see Figure 5). Again, one finds that the maximum speed up (for the lowest order moment, $n=2$) is close to the theoretical maximum of a factor 16.

5 Conclusion

In conclusion, this algorithm and method allows one to swiftly update moments and arbitrary metrics of datasets using less memory, and less data storage. The resultant code is lightweight, and can easily be implemented using class-oriented programming in a language of the readers’ choice (only Python results are shown here).

This algorithm is released as copyleft under the GNU GPL v3 license.

Appendix : Integer $n$ updates (proof)

To be demonstrated:

[TABLE]

where $Z^{\prime}$ is given by:

[TABLE]

and $\bar{x}^{\prime}$ is given by:

[TABLE]

Proof.

From (1), one finds that after update:

[TABLE]

where $\bar{x}^{\prime}$ is given by:

[TABLE]

First, if one examines $Z^{\prime}$:

[TABLE]

such that the first term is the original value of the normalising function, and the second term represents the change in the value from adding new data.

If we split (A-34) into a sum over $1$ to $N$, and $N+1$ to $N^{\prime}$:

[TABLE]

One finds that the following holds true:

[TABLE]

where we have performed a binomial expansion of modified form of $(x_{i}-\bar{x}^{\prime})^{n}$. By substituting the above into (A-33):

[TABLE]

where we noted that $M_{1}=0$.