# FlashP: An Analytical Pipeline for Real-time Forecasting of Time-Series
Relational Data

Shuyuan Yan  Alibaba Group<EMAIL_ADDRESS>,  Bolin Ding  Alibaba
Group<EMAIL_ADDRESS>,  Wei Guo  Alibaba Group
<EMAIL_ADDRESS>,  Jingren Zhou  Alibaba Group
<EMAIL_ADDRESS>,  Zhewei Wei  Renmin University of China
<EMAIL_ADDRESS>,  Xiaowei Jiang  and  Sheng Xu  Alibaba Group
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Interactive response time is important in analytical pipelines for users to
explore a sufficient number of possibilities and make informed business
decisions. We consider a forecasting pipeline with large volumes of high-
dimensional time series data. Real-time forecasting can be conducted in two
steps. First, we specify the part of data to be focused on and the measure to
be predicted by slicing, dicing, and aggregating the data. Second, a
forecasting model is trained on the aggregated results to predict the trend of
the specified measure. While there are a number of forecasting models
available, the first step is the performance bottleneck. A natural idea is to
utilize sampling to obtain approximate aggregations in real time as the input
to train the forecasting model. Our scalable real-time forecasting system
FlashP (Flash Prediction) is built based on this idea, with two major
challenges to be resolved in this paper: first, we need to figure out how
approximate aggregations affect the fitting of forecasting models, and
forecasting results; and second, accordingly, what sampling algorithms we
should use to obtain these approximate aggregations and how large the samples
are. We introduce a new sampling scheme, called GSW sampling, and analyze
error bounds for estimating aggregations using GSW samples. We introduce how
to construct compact GSW samples with the existence of multiple measures to be
analyzed. We conduct experiments to evaluate our solution and compare it with
alternatives on real data.

PVLDB Reference Format:
PVLDB, 14(5): XXX-XXX, 2021.
doi:10.14778/3446095.3446096 ††This work is licensed under the Creative
Commons BY-NC-ND 4.0 International License. Visit
https://creativecommons.org/licenses/by-nc-nd/4.0/ to view a copy of this
license. For any use beyond those covered by this license, obtain permission
by emailing<EMAIL_ADDRESS>Copyright is held by the owner/author(s).
Publication rights licensed to the VLDB Endowment.
Proceedings of the VLDB Endowment, Vol. 14, No. 5 ISSN 2150-8097.
doi:10.14778/3446095.3446096

Age ${\bf a}^{(1)}$ | Gender ${\bf a}^{(2)}$ | Location ${\bf a}^{(3)}$ | Impression ${\bf m}^{(1)}$ | ViewTime ${\bf m}^{(2)}$ | TimeStamp $t$
---|---|---|---|---|---
30 | F | WA | 5 | 1.6min | 20200301
60 | M | WA | 1 | 1.8min | 20200301
20 | F | NY | 10 | 3.2min | 20200301
40 | M | NY | 20 | 6.3min | 20200302
Figure 1. A time series of relational data (yellow cells are relevant to the
task)

Figure 2. Processing a real-time forecasting task

Figure 3. Forecasting example

##  1\. Introduction

Large volumes of high-dimensional data are generated on eCommerce platforms
every day, from data sources about, e.g., sales and browsing activities of
online customers. Forecasting is among the most important BI analytical tasks
to utilize such data, as sound prediction of future trends helps make informed
business decisions.

Motivating scenario and challenges. An online advertising platform enables
advertisers to target certain user visits and submit their bids for displaying
advertisements on these visits during a time window (i.e., targeting ads
campaigns). To make the right decision about which customers to target and the
bid prices, it is imperative for an advertiser to access reliable forecasts of
important measures about the targeted customers and their activities, e.g.,
Impression, the number of times an advertisement can be showed to such
customers, and ViewTime, the time a customer spends each visit.

Time-series data about customers can be very high-dimensional, with tens or
hundreds attributes about customers’ profiles and activities, including their
demographics, devices (e.g., mobile/PC), machine-learned tags (e.g., their
preferences and intents of visits), and so on. An advertiser may decide to
target a group of customers for any combination of the attributes and values,
based on the merchandise in the advertisement and fine-grained forecasts after
different ways of slicing and dicing in the attribute space. For example, she
may decide to target 20-30 year old females interested in sports and located
in some cities for skirt sets in the advertisement.

Consider a time series of relation data in Figure 3. A forecasting task, e.g.,
the one in Figure 3, asks to predict the total number of impressions by female
customers under age 30.

Unlike in traditional forecasting applications such as airline planning, real-
time response is critical in our scenario. First, it is an exploratory process
to find the right attribute combination for targeting. An advertiser may try a
number of combinations and read the forecasting results within a short period
of time in order to support the decision. High latency can easily make
advertisers impatient during the exploration. Second, real-time bidding for
targeting ads campaigns is very competitive and dynamic market. Slow response
may result in loss of displaying opportunities and higher prices. Thus, it is
important to make our platform interactive.

It is challenging to provide real-time response to such forecasting tasks. The
volume of input time-series data to the analytical pipeline is huge, with
hundreds of millions of rows per day, and commonly months of historical data
is used to forecast a future point. Moreover, as a result of the high-
dimensionality of our data, there could be trillions of possible
attribute/value combinations; a forecasting task is given online with any
combination, and it would be too expensive to precompute and store all
possible combinations and the corresponding time series during an offline
phase.

The task and solution phases. Our forecasting system FlashP is designed to
process real-time forecasting tasks, in which we specify: i) a constraint
specifying the portion of data to be focused on (e.g., targeting customers
“${\sf Age}$ $<=$ $30$ ${\color[rgb]{0,0,1}\tt AND}$ ${\sf Gender}$ $=$ ${\rm
F}$” in the example in Figure 3); ii) the aggregated measure to be predicted
(e.g., ${\tt SUM}({\sf Impression})$); iii) historical data points to be used
to train a forecasting model (e.g., from $20200101$ to $20200331$).

We give an overview of the solution. Consider the example in Figure 3. A
forecasting task is specified in a SQL-like language, with the goal to predict
the total number of impressions by female customers under age 30. To prepare
training data for a forecasting model, we need to know the total number of
impressions on each day ($M_{t}$), which can be written as an aggregation
query. As the time-series data is partitioned on time, we can process these 91
aggregation queries with one scan of the data. After that, we have 91 data
points to fit (train) the forecasting model. Processing the 91 aggregation
queries (or, equivalently, a query with ${\tt{\color[rgb]{0,0,1}\tt
GROUP}~{}{\color[rgb]{0,0,1}\tt BY}}$) is the bottleneck. A natural idea is to
utilize sampling (e.g., (jacm:DuffieldLT07, )) or sample-based approximate
query processing (e.g., (sigmod:DingHCC016, ; sigmod:ParkMSW18, )) to estimate
their answers. Figure 3 shows an example of prediction for the next 7 days:
the red line shows estimated aggregations; the estimations are used to train
the model, which then produces forecasts (green line) with confidence
intervals (green dashed lines).

Contributions and organization. FlashP is implemented in Alibaba’s advertising
system. The first technical question is how the errors in sample-based
estimations affect the fitting of forecasting models (Section 3). Error bounds
of estimating aggregations using samples are well studied e.g., for priority
sampling (stoc:Szegedy06, ) which has been shown to be optimal for aggregating
$\tt SUM$. However, it is unclear what the implications of such error bounds
are in prediction results. We provide both analytical and experimental
evidence showing that aggregation errors add up to the forecasting model’s
noise (from historical data points), and they together decide how confident we
could be with the prediction. We give formal analysis for a concrete
forecasting model (i.e., ${\rm ARMA}(1,1)$, defined later), and experimental
results for more complicated models (e.g., LSTM).

The second question is about the space budget needed for samples (Section 4).
Uniform sampling provides unbiased estimations for aggregations, but the
estimation error is proportional to the range of a measure ($\max-\min$)
(sigmod:Hellerstein97, ). Weighted sampling schemes offer optimal estimations,
with much better error bounds that are independent on the range
(stoc:Szegedy06, ); however, as the sampling distribution is decided by the
measure values (jacm:DuffieldLT07, ; pods:AlonDLT05, ), we have to draw one
weighted sample per measure independently. When there are a number of measures
(as in our scenario), the total space consumption could be prohibitive. We
propose a new sampling scheme called GSW (Generalized Smoothed Weighted)
sampling, with sampling distributions that can be arbitrarily specified. We
analyze estimation error bounds by quantify the correlation between the
sampling distribution and measure values. We introduce how to use this scheme
to generate a compact sample which takes care of multiple measures, using a
sampling distribution that averages distributions of different measures, and
analyze its estimation error bounds.

We describe how FlashP is implemented in Section 5. We report experimental
results on real datasets in Section 6, and discuss related work in Section 7.
All missing proofs are deferred to Appendix A.

##  2\. Solution Overview

Time-series data model. The input to our forecasting pipeline is a time series
of relation $T$, i.e., a sequence of observed rows at successive time. We
assume that time is a discrete variable here. A row in the table $T$ is
specified by a pair $(i,t)$: an item $i$ and a time stamp $t$. An item is
associated with multiple dimensions, each denoted by ${\bf a}$, which are used
to filter data (e.g., $\sf Age$ and $\sf Location$), and multiple measures,
each denoted by ${\bf m}$, which we want to analyze and forecast (e.g., $\sf
Impression$ and $\sf ViewTime$). We use ${a}_{i,t}$ and ${m}_{i,t}$ to denote
the values of a dimension and a measure, respectively, on each row; or,
simply, ${a}_{i}$ and ${m}_{i}$, if the time stamp $t$ is clear from the
context or not important. The schema of $T$ is $({\bf a}^{(1)},{\bf
a}^{(2)},\ldots,{\bf a}^{(d_{a})};{\bf m}^{(1)},{\bf m}^{(2)},\ldots,{\bf
m}^{(d_{m})};t)$.

Real-time forecasting task. A forecasting tasks is specified as:

(1) $\displaystyle{\color[rgb]{0,0,1}\tt FORECAST}~{}{\tt SUM}({\bf
m})~{}{\color[rgb]{0,0,1}\tt FROM}~{}T~{}{\color[rgb]{0,0,1}\tt WHERE}~{}{\sf
C}~{}{\color[rgb]{0,0,1}\tt USING}~{}(t_{s},t_{e})$
$\displaystyle{\color[rgb]{0,0,1}\tt OPTION}~{}({\tt
MODEL}={\tt{}^{\prime}model\\_x^{\prime}},{\tt FORE\\_PERIOD}={\tt
t\\_future})$

Here, ${\bf m}$ is the measure we want to forecast from data source $T$ on a
given set of rows satisfying the constraint ${\sf C}$. FlashP allows ${\sf C}$
to be any logical expression on the dimension values of ${\bf
a}^{(1)},\ldots,{\bf a}^{(d_{a})}$. We want to use historical data from time
stamp $t_{s}$ to $t_{e}$ to fit the forecasting model. We can also specify the
forecasting model we want to use in the ${\color[rgb]{0,0,1}\tt OPTION}$
clause with the parameter $\tt MODEL$, and the number of future time stamps we
want to predict with the parameter $\tt FORE\\_PERIOD$ (e.g., 7 days). In most
of our application scenarios, we care about $\tt SUM$ aggregation, e.g., total
number of impressions from a certain group users; FlashP can also support $\tt
COUNT$ and $\tt AVG$.

FlashP facilitates the following class of forecasting models. Let $M_{t}$ be
the value of metric $M$ ($={\tt SUM}({\bf m})$ in our case) at time $t$. A
model is specified by a time-dependent function $f_{t}$ with order $K$:

(2) $M_{t}=f_{t}(M_{t-1},M_{t-2},\ldots,M_{t-K}).$

The model $f_{t}$ is fitted on historical data $M_{1},M_{2},\ldots,M_{t_{0}}$
with training data tuples, each in the form of
$(M_{t};M_{t-1},\ldots,M_{t-K})$ with $M_{t-1},\ldots,M_{t-K}$ as the input to
$f_{t}$ and $M_{t}$ as the output, for $t=t_{0},t_{0}-1,\ldots,K+1$. The
fitted model can then be used to forecast future values
$M_{t_{0}+1},M_{t_{0}+2},\ldots$ as $\hat{M}_{t_{0}+h|t_{0}}$ for
$h=1,2,\ldots$, iteratively (i.e., $\hat{M}_{t_{0}+1|t_{0}}$ can be used
forecast $M_{t_{0}+2}$).

For the forecasting models supported by FlashP, we now give brief introduction
to both classic ones, e.g., ARMA (book:Hamilton94, ), and those based on
recurrent neural networks and LSTM (neco:HochreiterS97, ).

Forecasting using ARMA. An ARMA (autoregressive moving average)
(book:Hamilton94, ) model uses stochastic processes to model how $M_{t}$ is
generated and evolves over time. It assumes that $M_{t}$ is a noisy linear
combination of the previous $p$ values; each $u_{t}$ is an independent
identically distributed zero-mean random noise at time $t$, and historical
noise at previous $q$ time stamps impacts $M_{t}$ too:

(3) ${\rm
ARMA}(p,q):~{}M_{t}=\sum_{i=1}^{p}\alpha_{i}M_{t-i}+u_{t}+\sum_{i=1}^{q}\beta_{i}u_{t-i}.$

For example, an ARMA model is: $M_{t}=0.8M_{t-1}+0.2M_{t-2}+u_{t}+0.1u_{t-1}$.
It falls into the form of (2), and is parameterized by
$\alpha_{1}\ldots\alpha_{p}$ and $\beta_{1}\ldots\beta_{q}$. The model can be
fitted on $M_{1},\ldots,M_{t_{0}}$, using, e.g., least squares regression, to
find the values of $\alpha_{i}$ and $\beta_{j}$ which minimize the error term,
and used to forecast future values. We can estimate forecast intervals, i.e.,
confidence intervals for forecasts $\hat{M}_{t_{0}+h|t_{0}}$: with a
confidence level (probability) $\gamma$, the true future value $M_{t_{0}+h}$
is within
$[\hat{M}_{t_{0}+h|t_{0}}-l_{\gamma},\hat{M}_{t_{0}+h|t_{0}}+r_{\gamma}]$.

When there are deterministic trends over time, the differential method can be
used: the first order difference is defined $\triangledown
M_{t}=M_{t}-M_{t-1}$, and the second order
$\triangledown^{2}M_{t}=\triangledown M_{t}-\triangledown M_{t-1}$, and in
general, the $d$-th order
$\triangledown^{d}M_{t}=\triangledown^{d-1}M_{t}-\triangledown^{d-1}M_{t-1}$.
If $\\{\triangledown^{d}M_{t}\\}_{t}$ is an ${\rm ARMA}(p,q)$ model,
$\\{M_{t}\\}_{t}$ is an ${\rm ARIMA}(p,d,q)$ model.

Forecasting using LSTM. LSTM (long short-term memory) (neco:HochreiterS97, )
is a network architecture that extends the memory of recurrent neural networks
using a cell. It is natural to use LSTM to learn and memorize trends and
patterns of time series for the purpose of forecasting. In a typical LSTM-
based forecasting model, e.g., (ijcai:SchmidhuberWG05, ), an LSTM unit with
dimensionality $d$ in the output space takes $(M_{t-1},$ $\ldots,$ $M_{t-K})$
as the input; the LSTM unit is then connected to a $d\times 1$ fully-connected
layer which outputs the forecast of $M_{t}$. This model also falls into the
general form (2). Here, the cell state is evolving over time and, at each time
stamp $t$, is encoded in $f_{t}$, which can be learned from training data
tuples $(M_{t};$ $M_{t-1},$ $\ldots,$ $M_{t-K})$ in order.

###  2.1. Overview of Our Approach

Our system FlashP works in two online phases to process a forecasting task:
preparing training data points by issuing aggregation queries, and fitting the
forecasting model using the training data.

  * $\bullet$

Aggregation query. In a forecasting task, specified in (1), to predict ${\tt
SUM}({\bf m})$, we have $t_{e}-t_{s}+1$ historical data points:

$\displaystyle M_{t_{s}}$ $\displaystyle={\color[rgb]{0,0,1}\tt SELECT}~{}{\tt
SUM}({\bf m})~{}{\color[rgb]{0,0,1}\tt FROM}~{}T~{}{\color[rgb]{0,0,1}\tt
WHERE}~{}{\sf C}~{}{\color[rgb]{0,0,1}\tt AND}~{}t=t_{s}$
$\displaystyle\ldots$ $\displaystyle\ldots$ (4) $\displaystyle M_{t_{e}}$
$\displaystyle={\color[rgb]{0,0,1}\tt SELECT}~{}{\tt SUM}({\bf
m})~{}{\color[rgb]{0,0,1}\tt FROM}~{}T~{}{\color[rgb]{0,0,1}\tt WHERE}~{}{\sf
C}~{}{\color[rgb]{0,0,1}\tt AND}~{}t=t_{e}$

Each data point is given by an aggregation query with constraint ${\sf C}$,
which can be any logical expression on the dimension values of ${\bf
a}^{(1)},$ $\ldots,$ ${\bf a}^{(d_{a})}$. We would compute them in the online
phase.

  * $\bullet$

Forecasting. In the next online phase, we use $M_{t_{s}},\ldots,M_{t_{e}}$ as
training data to fit the forecasting model. And we use the model to predict
future aggregations, $M_{t_{e}+1},\ldots,M_{t_{e}+{\tt t\\_future}}$.

Performance bottleneck. Suppose we have $N$ rows in $T$ for each time stamp,
and use a history of $t_{0}=t_{e}-t_{s}+1$ time stamps to train a forecasting
model with size (number of weights) $s$. The total cost of processing a
forecasting task is ${\rm O}(t_{0}\cdot N)+{\rm Train}(t_{0},s)$, where ${\rm
O}(t_{0}\cdot N)$ is the cost of processing $t_{0}$ aggregation queries by
scanning the table $T$, and ${\rm Train}(t_{0},s)$ is the time needed to train
a forecasting model with size $s$ using $t_{0}$ training data tuples. ${\rm
Train}(t_{0},s)$ is in the form of, e.g., $(t_{0}\cdot s\cdot{\rm iter})$
where ${\rm iter}$ is the number of iterations for the model training to
converge. We typically have $t_{0}$ in hundreds ($t_{0}=365$ if we use one
year’s history for training with one data point per day), $N$ (number rows in
$T$ per day) in tens or hundreds of millions in our application, and $s$ in
tens. Therefore, as $t_{0},s\ll N$, processing of aggregation queries is the
performance bottleneck (even in comprison to training a complex model).

Real-time forecasting on approximate aggregations. In order to process a
forecasting task in an interactive way in FlashP, we propose to estimate
$M_{t_{s}},$ $\ldots,$ $M_{t_{e}}$ as $\hat{M}_{t_{s}},$ $\ldots,$
$\hat{M}_{t_{e}}$ from offline samples drawn from $T$, and use these estimates
to form training data tuples and to fit the model (2). There are several
questions to be answered. i) If estimates $\hat{M}_{i}$, instead of $M_{i}$,
are used to fit the forecasting model, how much the prediction would deviate.
ii) How to draw these samples efficiently (preferably in a distributed
manner). iii) How much space we need to store these samples for multiple
measures.

##  3\. Real-Time Forecasting

We first analyze how sampling and approximate aggregations impact model
fitting and, resultingly, forecasts. We give an analytical result for a
special case of ARMA model, and will conduct experimental study for more
complex models, e.g., LSTM-based ones in Section 6. It is not surprising that
there is a tradeoff between sampling rate (or, quality of estimated
aggregations) and forecast accuracy.

Required properties of sampling and estimates. In FlashP, we have no access to
the accurate value of $M_{t}$; but instead, we have $\hat{M}_{t}$ estimated
from offline samples. We require the estimates to have two essential
properties for the fitting of forecast models:

  * $\bullet$

(Unbiasedness) $\hat{M}_{t}=M_{t}+\epsilon_{t}$ with ${\bf
E}[\epsilon_{t}]=0$;

  * $\bullet$

(Independence) $\epsilon_{t}$’s for different time stamps $t$ are independent.

GSW sampling introduced in Section 4 offers unbiased estimates, with bounded
variance of $\epsilon_{t}$. GSW samplers are run independently on the data for
each $t$—that is how we get independence.

Impact on ${\rm ARMA}(p,q)$. When an ${\rm ARMA}(p,q)$ model has to be trained
only on noisy metric values $\\{\hat{M}_{t}\\}$, we rewrite (3) as

(5) $\displaystyle\hat{M}_{t}$
$\displaystyle=\sum_{i=1}^{p}\alpha_{i}\hat{M}_{t-i}+(u_{t}+\epsilon_{t})+$
$\displaystyle+\begin{cases}\sum_{i=1}^{p}(\beta_{i}u_{t-i}-\alpha_{i}\epsilon_{t-i})+\sum_{i=p}^{q}\beta_{i}u_{t-i}&(1<p\leq
q)\\\
\sum_{i=1}^{q}(\beta_{i}u_{t-i}-\alpha_{i}\epsilon_{t-i})-\sum_{i=q}^{p}\alpha_{i}\epsilon_{t-i}&(1<q\leq
p)\end{cases}$

In both cases, it is a combination of an autoregressive model (5) of order
$p$, and a moving average model of order $\max\\{p,q\\}$. The model differs
from the ARMA model in (3) only on the additional zero-mean error terms
$\epsilon_{t}$, which are independent on other terms in the model and have
known variance (for fixed sampling and the estimation methods). Thus, the
model can be fitted on $\\{\hat{M}_{t}\\}$ using, e.g., maximum likelihood
estimator, as normal ARMA models. $\epsilon_{t}$ increases the model’s
uncertainty, and, together with the model noise terms $u_{t}$, decides the
confidence of the model prediction.

For a fixed confidence level, the narrower the forecast intervals are, the
more confident we are about the prediction. Their widths are proportional to
the standard deviations of forecasts of $\hat{M}_{t}$, which in turn depends
on the variance of noise terms $u_{t}$’s and $\epsilon_{t}$’s.

In comparison to the original ARMA model in (3), the additional noise term
$\epsilon_{t}$ will indeed incur wider forecast intervals. However, with a
proper sampling-estimation scheme, if $\epsilon_{t}$’s variance is negligible
in comparison to $u_{t}$’s, $\epsilon_{t}$ will have little impact on the
forecast error/interval, which will be demonstrated in our experiments later.
Here, we give a formal analysis for the ${\rm ARMA(1,1)}$ model:

###### Proposition 1.

Suppose we have a time series $\\{M_{t}\\}$ satisfies ${\rm ARMA}(1,1)$. Let
$\hat{M}_{t}=M_{t}+\epsilon_{t}$ be an estimation of $M_{t}$ satisfying
unbiasedness and independence. Then, ${\bf
Var}[\hat{M}_{t}]=a\cdot\sigma_{u}^{2}+\sigma_{\epsilon}^{2}$, where
$\sigma_{u}^{2}={\bf Var}\left[u_{t}\right]$, $\sigma_{\epsilon}^{2}={\bf
Var}\left[\epsilon_{t}\right]$, and
$a=(1+2\alpha_{1}\beta_{1}+\beta_{1}^{2})/(1-\alpha_{1}^{2})$ is a constant
decided by parameters in ${\rm ARMA}(1,1)$.

We can use normal random variables to approximate $u_{t}$ and $\epsilon_{t}$,
and estimate forecast intervals for a given confidence level.

Impact on LSTM-based model. LSTM can be naturally applied for forecasting
tasks thanks to its ability to memorize trends and patterns of time series.
Figure 4 depicts such a forecasting model and illustrates where noise in the
estimates impacts the model fitting.

At time stamp $t$, we want to forecast $M_{t}$ with the previous $K$ metric
values $M_{t-1},\ldots,M_{t-K}$ and the “memory” $({\bf c}_{t-1},{\bf
h}_{t-1})$. However, we have only estimates
$\hat{M}_{t-i}=M_{t-i}+\epsilon_{t-i}$ available; these estimates are fed into
the LSTM unit as inputs. LSTM then generates an output vector ${\bf h}_{t}$
and update its memory cell from ${\bf c}_{t-1}$ to ${\bf c}_{t}$. ${\bf
h}_{t}$ can be interpreted as the current status of LSTM and it is used as the
input to a fully-connected layer for forecasting $M_{t}$. Again, we have only
$\hat{M}_{t}=M_{t}+\epsilon_{t}$ available as an approximation, which we learn
to fit with LSTM and the fully-connected layer. ${\bf c}_{t}$ and ${\bf
h}_{t}$ are used to deliver memory to the next time stamp.

Noise terms $\epsilon_{t}$’s may affect how the weight matrices in LSTM and
the fully-connected layer are learned and the values of ${\bf c}_{t}$ and
${\bf h}_{t}$ derived (via linear transformations and activation functions).
We conjecture that the LSTM-based forecasting model performs well as long as
the estimates $\hat{M}_{t}$’s are accurate enough (or $\epsilon_{t}$’s are
small enough). It is difficult to derive any formal analytical result here,
but we will evaluate it experimentally in Section 6.

Figure 4. An LSTM-based forecasting model with noisy inputs

##  4\. Generalized Weighted Sampler

We now focus on how to draw samples for estimating results of aggregation
queries. An aggregation query in (4) is essentially to compute the sum of a
subset of measure values in a relation $T$ at time $t$; the subset is decided
online by the constraint ${\sf C}$ specified in the forecasting task.
W.l.o.g., suppose the subset of rows satisfying ${\sf C}$ is
$[n]=\\{1,\ldots,n\\}$, we want to estimate the metric
$M=\sum_{i=1}^{n}m_{i}$, for a measure ${\bf m}=[{m}_{i}]$. The (offline)
sampling algorithm should be independent on ${\sf C}$, but could depend on
${\bf m}$.

Existing samplers. There are two categories of sampling schemes for estimating
subset sums: uniform sampling and weighted sampling. In uniform sampling, each
row in the relation is drawn into the sample with equal probability; an
unbiased estimation for $M$ is simply the rescaled sum of values in the
sample. It has been used in online aggregation extensively, but the main
deficiency is that the error bound is proportional to the difference between
the maximum and minimum (or, the range of) measure values
(sigmod:Hellerstein97, ).

In weighted sampling, each row $i$ is drawn with probability proportional to
$m_{i}$, so that we can remove the dependency of the estimation’s error bound
on the range of the measure values. More concretely, (sigmod:chaudhuri99, )
and (sigmod:DingHCC016, ) give efficient implementations of such a sampling
distribution: for some fixed constant $\tau$ (deciding the sampling size), if
$m_{i}<\tau$, the probability of drawing a row $i$ is $m_{i}/\tau$, and if
$m_{i}\geq\tau$, multiple (roughly $\tau/m_{i}$) copies of $i$ are included
into the sample. Threshold sampling (tit:DuffieldLT05, ) and priority sampling
(pods:AlonDLT05, ; jacm:DuffieldLT07, ) differ from the above one in that, if
$m_{i}\geq\tau$, exactly one copy of $i$ is included. Priority sampling has
been shown to be optimal (stoc:Szegedy06, ) in terms of the sampling
efficiency with relative standard deviation $\sqrt{{\bf Var}[\rm
estimation]}/M=\sqrt{1/({\rm sample\\_size-1}})$.

Requirements for the sampler in FlashP. Weighted sampling offers much better
sampling efficiency than uniform sampling, especially for heavy tailed
distributions, which is common in practice. However, in all the above weighted
sampling schemes, the sampling distribution is decided by the measure’s
values. In FlashP, we have $d_{m}$ different measures to forecast; and thus,
we have to draw $d_{m}$ such weighted samples independently; when $d_{m}$ is
large (e.g., $10$), the storage cost of all the samples (e.g., even with
sampling rate $1\%$) is prohibitive in memory. Therefore, the requirement here
is: whether we can generate a compact sample which takes care of multiple
measures and still have accuracy guarantees.

To this end, we propose GSW (Generalized Smoothed Weighted) sampling. We
introduce its sampling distribution and analyze its sampling efficiency in
Section 4.1. We utilizes the “generality” of its sampling distribution to
generate a compact sample which takes care of multiple measures, and analyzes
under which conditions it gives provable accuracy guarantees in Section 4.2.

###  4.1. Generalized Smoothed Sampling

The GSW sampling scheme is parameterized by a positive constant $\Delta$ and
positive sampling weights ${\bf w}=[w_{i}]$: each row $i$ in the given
relation $T$ is drawn into the sample $S_{\Delta}$ with probability
$\frac{w_{i}}{\Delta+w_{i}}$, independently. For fixed sampling weights,
$\Delta$ decides the sample size; and the choice of $\bf w$ decides the
estimation accuracy.

Inspired by the Horvitz–Thompson estimator (jasa:HorvitzT52, ), we define the
calibrated measure to be $\hat{m}_{i}=m_{i}(\Delta+w_{i})/w_{i}$ if row $i$ is
drawn into the sample $S_{\Delta}$, and $\hat{m}_{i}=0$ otherwise. We would
associate $\hat{m}_{i}$ with each sample row $i$ and store it in $S_{\Delta}$.
Formally, we have

(6) $\begin{cases}{\bf Pr}[i\in
S_{\Delta}\wedge\hat{m}_{i}=\frac{m_{i}(\Delta+w_{i})}{w_{i}}]=\frac{w_{i}}{\Delta+w_{i}}\\\
{\bf Pr}[i\notin
S_{\Delta}\wedge\hat{m}_{i}=0]=1-\frac{w_{i}}{\Delta+w_{i}}=\frac{\Delta}{\Delta+w_{i}}\end{cases}.$

We can estimate $M$ as $\hat{M}=\sum_{i\in S_{\Delta}}\hat{m}_{i}$, which is
unbiased since

${\bf E}\left[{\hat{M}}\right]=\sum_{i=1}^{n}{\bf
E}\left[\hat{m}_{i}\right]=\sum_{i=1}^{n}\frac{m_{i}(\Delta+w_{i})}{w_{i}}\cdot\frac{w_{i}}{\Delta+w_{i}}=M.$

Simple and efficient implementations. A GSW sampler can be easily implemented
in a distributed manner: each row $i$ generates uniformly random number
$p_{i}$ from $[0,1]$, independently; according to (6), if
$p_{i}\leq\frac{w_{i}}{\Delta+w_{i}}$, the row $i$ is put into the sample
$S_{\Delta}$.

A GSW sample $S_{\Delta}$ can be maintained in an incremental way. Suppose we
have drawn $S_{\Delta}$ from rows $[n]$ for some fixed $\Delta$. If more rows
$n+1,\ldots,n+k$ are coming, we want to increase $\Delta$ to $\Delta^{\prime}$
and obtain a GSW sample $S_{\Delta^{\prime}}$ from $[n+k]$ with size
comparable to $|S_{\Delta}|$. Suppose rows in $S_{\Delta}$ are sorted by
$(\frac{1}{p_{i}}-1)w_{i}$ in an ascending order; we only need to delete those
with $\Delta\leq(\frac{1}{p_{i}}-1)w_{i}<\Delta^{\prime}$ from $S_{\Delta}$,
and insert those with $\Delta^{\prime}\leq(\frac{1}{p_{i}}-1)w_{i}$, for
$i=n+1,\ldots,n+k$. In this way, we update $S_{\Delta}$ to
$S_{\Delta^{\prime}}$ without touching any row in $[n]-S_{\Delta}$.

####  4.1.1. Accuracy Guarantee on Aggregations

We now analyze estimation errors of the class of GSW-based estimators
$\hat{M}$, instantiated by $(\Delta,{\bf w})$. We consider relative standard
deviation ($\rm RSTD$) and relative error ($\rm RE$). As $\hat{M}$ is
unbiased,

$\displaystyle{\rm RSTD}({\hat{M}})\\!\triangleq\\!\\!\sqrt{{\bf
E}\left[\left(\frac{{\hat{M}}-M}{M}\right)^{2}\right]}\geq{\bf
E}\left[\frac{|{\hat{M}}-M|}{M}\right]\\!\\!\triangleq\\!{\rm RE}({\hat{M}}).$

The goal of our analysis is to establish a relationship between the sample
size and the (expected) $\rm RSTD$ and $\rm RE$. Intuitively, when the
sampling weight $w_{i}$ is “consistent” with the measure $m_{i}$, the
estimation error is minimum (for fixed sample size). We introduce the
following notation to quantify the “consistency”.

###### Definition 2.

($(\underline{\theta},\overline{\theta})$-consistency) Sampling weights ${\bf
w}$ $=$ $[w_{i}]$ are $(\underline{\theta},\overline{\theta})$-consistent with
measure ${\bf m}=[{m}_{i}]$, iff $\underline{\theta}=\min_{i}(m_{i}/w_{i})$
and $\overline{\theta}=\max_{i}(m_{i}/w_{i})$.
$\theta\triangleq\overline{\theta}/\underline{\theta}$ is called the
consistency scale.

The above notation about “consistency” says, for any row $i$,
$\underline{\theta}\leq m_{i}/w_{i}\leq\overline{\theta}$. It allows ${\bf m}$
and ${\bf w}$ to differ in scale (e.g., both $\underline{\theta}$ and
$\overline{\theta}$ could be large) but measures their similarity in patterns
and trends. For example, suppose ${\bf m}$ $=$ $[100,$ $100,$ $200,$ $400]$
and ${\bf w}$ $=$ $[10,$ $10,$ $20,$ $50]$. We have
$\underline{\theta}=400/50=8$, $\overline{\theta}=10$, and thus
$\theta=10/8=1.25$. In general, we have $\theta\geq 1$, and the following
theorem shows that, the relative error has an upper bound that is proportional
to $\sqrt{\theta}$.

###### Theorem 3.

(Sampling efficiency of GSW) Suppose the sampling weights $\bf w$ used in GSW
sampling are $(\underline{\theta},\overline{\theta})$-consistent with measure
values $\bf m$, the estimate $\hat{M}$ is unbiased and has expected relative
standard deviation and error bounded by:
($\theta\triangleq\overline{\theta}/\underline{\theta}$)

${\rm RE}({\hat{M}})\leq{\rm
RSTD}({\hat{M}})\leq\sqrt{\frac{\overline{\theta}/\underline{\theta}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}=\sqrt{\frac{\theta}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}.$

An open question raised by Alon et al. when introducing priority sampling
(pods:AlonDLT05, ) is: whether we can provide any error bound if subset sum on
a measure (${\bf m}$) is estimated using a priority sample drawn based on a
different measure ($\bf w$). Theorem 3 answers this question in GSW sampling
by specifying under what condition
($(\underline{\theta},\overline{\theta})$-consistency)) there is an error
bound.

####  4.1.2. A Special Case: Optimal GSW Sampler

We can choose $\bf w$ to minimize the variance of estimation, while the
expected sample size is bounded. Please refer to Appendix A.3 for a precise
formulation and the optimal solution. From Theorem 3, a nearly-optimal
solution is ${\bf w}={\bf m}$, as in this case, $\bf w$ is $(1,1)$-consistent
with $\bf m$ ($\theta=1$). We call it optimal GSW sampler. Directly from
Theorem 3, we have

###### Corollary 4.

(Optimal GSW sampler) If we use ${\bf w}={\bf m}$ ($w_{i}=m_{i}$ for each row
$i$) as the sampling weights, we have

${\rm RE}({\hat{M}})\leq{\rm RSTD}({\hat{M}})\leq\sqrt{\frac{1}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}.$

The optimal GSW sampler has sampling efficiency that is comparable to the best
known sampler for estimating subset sums, e.g., priority sampling
(pods:AlonDLT05, ) with ${\rm RSTD}=\sqrt{1/({\rm sample\\_size-1})}$. We will
compare their empirical performance in Section 6.

###  4.2. Compact Sample for Multiple Measures

The size of GSW sample can be controlled by the parameter $\Delta$. When there
is only one measure ${\bf m}$, we draw an optimal GSW sample (setting ${\bf
w}={\bf m}$). A more common scenario is that we have multiple measures (e.g.,
in Figure 3) in one relation. We can draw one optimal GSW sample per each
measure, which, however, increases the space consumption significantly. The
question is whether we can use one sample to take care multiple measures.

Suppose there are $k$ measures in a relation to be aggregated and predicted,
${\bf m}^{(1)},$ $\ldots,$ ${\bf m}^{(k)}$. For each measure $j$, we want to
estimate $M^{(j)}=\sum_{i=1}^{n}{m}^{(j)}_{i}$ for rows in a set $[n]$
(satisfying the constraint ${\sf C}$ in a forecasting task). A GSW sample
$S_{\Delta}$ is drawn using weights ${\bf w}=[w_{i}]$. The calibrated measure
on each sample row $i$ for each measure $j$ is
$\hat{m}^{(j)}_{i}=m^{(j)}_{i}(\Delta+w_{i})/w_{i}$. From Theorem 3,
$\hat{M}^{(j)}=\sum_{i\in S_{\Delta}}\hat{m}^{(j)}_{i}$ is an unbiased
estimation of $M^{(j)}$.

Intuitively, if the chosen sampling weight vector $\bf w$ centers around ${\bf
m}^{(1)},$ $\ldots,$ ${\bf m}^{(k)}$ and is not too far away from any of the
$k$, from Theorem 3, the error can be better bounded. We can find such centers
by taking the average of measures. For example, for ${\bf
m}^{(1)}=[100,100,200,400]$ and ${\bf m}^{(2)}=[1,1,2,1]$, the geometric mean
is ${\bf w}^{\times}=[\sqrt{100\cdot 1}=10,10,20,20]$, and the arithmetic mean
is ${\bf w}^{+}=[(100+1)/2=50.5,50.5,101,200.5]$. We now analyze how the error
can be bounded for these two choices.

Geometric compressed GSW sampling. We can use the geometric mean of the $k$
measures as the sampling weights:

(7) $w^{\times}_{i}=\left(\prod_{j=1}^{k}{m}^{(j)}_{i}\right)^{1/k}.$

Among the $k$ measures to be grouped, define the trend deviation between any
two measures ${\bf m}^{(p)}$, ${\bf m}^{(q)}$ (for $p,q\in[k]$):

(8)
$\overline{\rho}_{p,q}\triangleq\max_{i\in[n]}\frac{{m}^{(p)}_{i}}{{m}^{(q)}_{i}},~{}~{}~{}\underline{\rho}_{p,q}\triangleq\min_{i\in[n]}\frac{{m}^{(p)}_{i}}{{m}^{(q)}_{i}},~{}~{}~{}\hbox{and}~{}\rho_{p,q}\triangleq\frac{\overline{\rho}_{p,q}}{\underline{\rho}_{p,q}},$

which measures how similar the trends (instead of their absolute values) of
the two measures are111Or, equivalently, ${\bf m}^{(q)}$ is
$(\underline{\rho}_{p,q},\overline{\rho}_{p,q})$-consistent with ${\bf
m}^{(p)}$.. The smaller $\rho_{p,q}$ is, the more similar ${\bf m}^{(p)}$ and
${\bf m}^{(q)}$ are. For example, if ${\bf m}^{(p)}=c\cdot{\bf m}^{(q)}$ for a
constant $c$, $\rho_{p,q}=c/c=1$. Define
$\rho\triangleq\max_{p,q\in[k]}\rho_{p,q}$ to be the maximum deviation among
the $k$ measures. From Theorem 3,

###### Corollary 5.

If we use ${\bf w}^{\times}=[w^{\times}_{i}]$ as the sampling weights for a
relation with $k$ measures, with a GSW sample $S_{\Delta}$, we can estimate
$M^{(p)}$ as $\hat{M}^{(p)}$ for each measure $p$ with error

${\rm RE}({\hat{M}}^{(p)})\leq{\rm
RSTD}({\hat{M}}^{(p)})\leq\sqrt{\frac{\prod_{j:~{}j\neq
p}\rho_{p,j}^{1/k}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}\leq\sqrt{\frac{\rho^{\frac{k-1}{k}}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}.$

Arithmetic compressed GSW sampling. Another choice is to use the arithmetic
mean as the sampling weights:

(9) $w^{+}_{i}=\frac{1}{k}\sum_{j=1}^{k}{m}^{(j)}_{i}.$

Define the range deviation $\delta$ among $k$ measures: for each row $i$,
consider the ratio between the maximum measure and the minimum one; $\delta$
is the maximum ratio among all rows:

(10)
$\delta\triangleq\max_{i\in[n]}\left(\frac{\max_{j\in[k]}{m}^{(j)}_{i}}{\min_{j\in[k]}{m}^{(j)}_{i}}\right).$

From the definitions, for any measure $p$, we have $1/\delta$ $\leq$
${m}^{(p)}_{i}/w^{+}_{i}$ $\leq$ $\delta$. Thus, directly from Theorem 3, we
have the following bound.

###### Corollary 6.

If we use ${\bf w}^{+}=[w^{+}_{i}]$ as the sampling weights for a relation
with $k$ measures, with a GSW sample $S_{\Delta}$, we can estimate $M^{(p)}$
as $\hat{M}^{(p)}$ for each measure $p$ with error

${\rm RE}({\hat{M}}^{(p)})\leq{\rm
RSTD}({\hat{M}}^{(p)})\leq\sqrt{\frac{\delta^{2}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}.$

How to group measures? In the above, we have shown that, for chosen sampling
weights, how the error can be bounded with some parameters ($\rho$ and
$\delta$) decided by the data about a group of measures. When there are many
(e.g., $k=100$) measures, $\rho$ and $\delta$ could be huge and thus the above
error bounds are not informative. We want to partition measures into small
groups of size 4-5 based on their correlation, so that within each group one
GSW sample gives good estimations for the measures. To this end, we establish
a relationship between $(\underline{\theta},\overline{\theta})$-consistency
and $L_{1}$ distance as follows.

###### Proposition 7.

(consistency and $L_{1}$ distance) Suppose sampling weights ${\bf w}$ $=$
$[w_{i}]$ are $(\underline{\theta},\overline{\theta})$-consistent with measure
${\bf m}=[{m}_{i}]$. We normalize $\bf w$ as ${\bf
w}^{\prime}=[w^{\prime}_{i}=w_{i}/\sum_{j}w_{j}]_{i}$, and similarly, ${\bf
m}$ as ${\bf m}^{\prime}$. Let $\theta=\overline{\theta}/\underline{\theta}$.
We have $\|{\bf m}^{\prime}-{\bf w}^{\prime}\|_{1}\leq(\theta-1)$.

$(\underline{\theta},\overline{\theta})$-consistency is a “worst-case” notion
about all the rows. It is not a good metric for grouping measures, because
first, it is expensive to compute $\theta$, and second, an aggregation query
may not touch all the rows. Instead, we use $L_{1}$ distance (after
normalization) as in Proposition 7 to quantify the correlation between two
measures. Intuitively, it quantifies how two measures are similar in trends
and patterns regardless of their absolute values.

We consider a formulation based on the KCenter problem (book:Harpeled2011, ).
The goal is to partition the measures into $g$ groups so that the max $L_{1}$
distance (after normalization) between any measure to the center of the group
it belongs to is minimized. Here, the $L_{1}$ distance between two measures
can be estimated using a sample of rows. We can apply the standard greedy
algorithm (book:Harpeled2011, ) to find a 2-approximation. Within each group,
we use the geometric/arithmetic mean as the sampling weight vector. It is a
heuristic strategy without any formal guarantee, but from Proposition 7 and
the triangle inequality in $L_{1}$, at least we know that we’d better not
group two measures that are far way (e.g., $>2(\theta-1)$) together, as there
is no sampling weight vector that is consistent with both of them at the same
time. We conduct a preliminary evaluation on the above correlation metric,
$L_{1}$ distance (after normalization). Using the same datasets and workloads
(with sensitivity ranging from $0.5\%$ to $10\%$) as Exp-I in Section 6, the
experimental results are plotted in Figure 5. Here, the four measures,
Impression, Click, Favorite, and Cart, are partitioned into two equal-size
groups (there are a total of three ways of partitioning). For example, the
first way in Figure 5 is group1 = {Impression, Click} and group2 = {Favorite,
Cart}. For each group, we pick the sampling weight vector as the arithmetic
mean of measure vectors in this group. Figure 5(b) plots the $L_{1}$ distance
between each measure vector and the sampling weight vector, and Figure 5(a)
plots the aggregation error using this sampling weight vector in GSW. It can
be seen that aggregation errors and $L_{1}$ distances have similar trends.
That is, if a sampling weight vector is close to a measure vector, it gives
better estimation; thus, if two measure vectors are close in $L_{1}$, it is
better to put them in one group (their mean can be used as the sampling weight
vector). In Section 6, we focus on experimental evaluation after the grouping
is fixed. We will leave the development and analysis of better grouping
strategies as future work.

(a) Three ways to group the four measures: Imp = Impression, clk = Click, fav
= Favorite, cart = Cart

(b) $L_{1}$ distance between measure vectors and the sampling weight vector
(the arithmetic mean of measure vectors in a group)

Figure 5. Aggregation error and $L_{1}$ distance Figure 6. Deployment of the
system FlashP

##  5\. Deployment

We have implemented and deployed FlashP in Alibaba’s ads analytical platform.
Figure 6 shows important components of FlashP.

Offline Sample Preprocessor of FlashP is built on Alibaba’s distributed data
storage and analytics service, MaxCompute (url:alibabacloud, ). Time series of
relations, partitioned by time, are stored in MaxCompute’s data warehouse.
Relations can be joined, e.g., tables $\sf UserProfile$ and $\sf AdTraffic$
are joined on $\sf UserID$. GSW sampler is implemented as UDFs (user-defined
functions), and draws samples from one relation or the view of joined
relation. A set of samples of different sizes (with increasing values of
$\Delta$) are drawn and stored as Multi-layer Samples of Relations for
response time-prediction accuracy tradeoff.

Online Forecasting Service first pulls Multi-layer Samples of Relations into
Alibaba’s in-memory OLAP engine, Hologres (url:alibabacloud, ), to enable
real-time response. There are 30 servers in the cluster to support this
service, each with 96 CPU cores and 512G memory. Sample data is partitioned by
time in the OLAP engine.

Users submit forecasting tasks to an Application Server via a Web UI. For a
task, aggregation queries in (4) are generated from Query Rewriter and
processed on samples in the OLAP engine to obtain estimated answers. These
estimations are used as training data to fit the Forecasting Model, which is
built on a Python server.

Two models are implemented. One is ARIMA. An open-source library
(url:autoarima, ) (built on X-13ARIMA-SEATS (url:x13arima, )) that trains
ARIMA and automatically tunes for the best values of parameters $p,d,q$ is
used. We also support an LSTM-based model (in Figure 4), which is implemented
using Keras (url:keras, ). We use the LSTM and fully-connected layers in Keras
with $K=7$ and $d=4$ as the default parameter setting. Other forecasting
models can be plugged in here, too. After fitting the model, we send forecasts
back to users.

##  6\. Experimental Evaluation

We evaluate our system FlashP under the implementation and hardware specified
in Section 5, on a real-life dataset with 11 dimensions about users’ profiles
and 4 numeric measures to be predicted, including $\sf Favorite$, $\sf
Impression$, $\sf Click$, and $\sf Cart$. There are around 15 million rows per
day, and 200 days of data.

We compare different samplers: Uniform sampling, which is also used in
(sigmod:AgarwalCLSV10, ), is the baseline; Priority, the optimal weighted
sampler (jacm:DuffieldLT07, ); our Optimal GSW and Arithmetic/Geometric
compressed GSW introduced in Section 4. We also compare our sampling based
methods with PIM (Partwise Independence Model) (sigmod:AgarwalCLSV10, ) based
on a Bayesian model assuming partially independence.

In the following, forecasting tasks are randomly picked with different
measures to be predicted and different combinations of dimensions in their
constraints, with some (approximately) fixed selectivity (the fraction of rows
satisfying the constraint). By default, we use 150 days’ data to fit the model
and predict the next 7 days; we report relative aggregation errors (average of
the 150 days), relative forecast error, and forecast intervals (average of the
next 7 days), taking the average of 400 independent runs of different tasks,
together with one standard deviation, for each measure and for each value of
selectivity (on independent samples).

Exp-I: A summary of results. We first give a brief summary of experimental
results. Table 1 reports the average forecast errors on 20 random tasks with
selectivity from $0.5\%$ to $10\%$ using ARIMA. “Full” stands for the result
when we use the full data to process aggregation queries for training. With a
sampling rate $0.1\%$, our optimal GSW and compressed GSW in FlashP perform
consistently better than Uniform and PIM in terms of forecast errors, and
sometimes are very close to Full (the best we can do). These two also offer
interactive response time (less than $100$ms). In the rest part, we will
report more detailed results about response time and performance of different
sampling-based methods.

Table 1. A summary of results (0.1% sample, Opt-GSW = Optimal GSW, C-GSW = Arithmetic compressed GSW) | Full | PIM | Uniform | Opt-GSW | C-GSW
---|---|---|---|---|---
$\sf Favorite$ | 0.105 | 0.695 | 0.248 | 0.131 | 0.196
$\sf Impression$ | 0.140 | 0.374 | 0.147 | 0.142 | 0.144
$\sf Click$ | 0.157 | 0.681 | 0.161 | 0.151 | 0.153
$\sf Cart$ | 0.704 | 1.931 | 0.718 | 0.704 | 0.709
Figure 7. E2E response time (with ARIMA)

Exp-II: Real-time response. The end-to-end response time is reported in Figure
7, partitioned into the portion for processing (estimating) aggregation
queries, and the portion for forecasting (model fitting + prediction using
ARIMA). It can be seen that the portion for aggregation queries is the
bottleneck, but with sampling, the response time can be reduced from around
20sec on the full data to 30ms on a $0.02\%$ sample, which still gives
reasonable prediction as will be shown later. If LSTM is used, the model
fitting is much more expensive, but we still have an interactive response time
around 1sec if a $1\%$ sample is used.

Exp-III: Varying number of time stamps used in training. We consider different
numbers of training data points used to fit the model. We focus on tasks with
selectivity $5\%$ to forecast Impression. The average forecast errors, with
one standard deviation, for varying sampling rate using our optimal GSW are
plotted in Figure 8. It shows that the number of time stamps used to fit model
has an obvious impact on the forecast error, with 150 (days) giving the most
accurate and stable prediction for both ARIMA and LSTM. It motivates us to
speedup the processing of aggregation queries, as more time stamps mean more
aggregation queries. We also test different selectivities on other measures.
The trends are similar.

(a) ARIMA

(b) LSTM

Figure 8. Number of time stamps used in training v.s. Forecast error
(selectivity $5\%$, $\sf Impression$)

(a) Selectivity 0.5%

(b) Selectivity 5%

Figure 9. Aggregation error of different sampling methods for varying
selectivity and sampling rate (Favorite)

(a) ARIMA

(b) LSTM

Figure 10. Forecast error of different sampling methods for varying sampling
rate (selectivity 0.5%, Favorite)

(a) ARIMA

(b) LSTM

Figure 11. Forecast error of different sampling methods for varying sampling
rate (selectivity 5%, Favorite)

Exp-IV: Varying sampling rate and selectivity. We compare different samplers
in FlashP, for varying sampling rate and selectivity. Note that for
Arithmetic/Geometric compressed GSW, one sample suffices for the relation;
Priority are Optimal GSW are measure-dependent, and thus using either of them
we need four samples (one per measure), with the total space consumption four
times of Uniform and compressed GSW for a fixed sampling rate.

For tasks on $\sf Favorite$ with selectivity $0.5\%$-$5\%$, Figure 9 reports
aggregation errors, and Figures 10-11 report forecast errors when using ARIMA
and LSTM as forecasting models; the results on $\sf Impression$ are plotted in
Figures 13-14. In terms of both aggregation errors and forecasting errors,
Priority and Optimal GSW are very close and better than the others (indeed, at
the cost of storing four samples). In some cases, Optimal GSW is even slightly
better than Priority (theoretically optimal). This is because, in Priority, if
the measure is above some threshold, a row is included in the sample
deterministically, which favors the long tail; however, the long tail may or
may not satisfy the constraint specified online in the forecasting task.
Uniform is the worst one which is consistent with its analytical error bound
(sigmod:Hellerstein97, ). Arithmetic/Geometric compressed GSW needs only one
sample, too; they are better than Uniform, and get very close to Priority and
Optimal GSW when the sampling rate is close to $1\%$. For larger selectivity,
with more rows satisfying the constraint included in the samples, every
sampler gets better.

Figure 12(a) reports widths of forecast intervals of ARIMA with a confidence
level of $90\%$ on measure $\sf Favorite$ for varying sampling rate. With
larger sampling rate, every sampler gives narrower forecast intervals, meaning
predictions with more confidence. Uniform gives the widest one and Priority
and Optimal GSW give the narrowest ones. Figure 12(b) shows the forecast
intervals (dashed lines) for one particular task using different samplers.

In terms of forecasting, LSTM performs consistently better than ARIMA, at the
cost of longer response time (as discussed in Exp-II). It is observed that,
with increasing sampling rates, when the aggregation error is small enough
(e.g., when sampling rate = $1\%$ in Figures 9-13), it will have little impact
on the forecast error/interval in comparison to the case when we use the full
data (sampling rate = $100\%$), because it is negligible in comparison to the
model and data’s noise (e.g., $u_{t}$ in (3) for ARIMA). Forecast errors and
forecast intervals have similar trends as aggregation errors. Both
observations are consistent with our analytical results in Section 3.

(a) Varying sampling rate

(b) Forecast intervals of a real query on 0.02% sampling

Figure 12. Forecast intervals (ARIMA) with different sampling methods for
varying sampling rate (selectivity 0.5%, Favorite)

(a) ARIMA

(b) LSTM

Figure 13. Forecast error of different sampling methods for varying sampling
rate (selectivity 0.5%, Impression)

(a) ARIMA

(b) LSTM

Figure 14. Forecast error of different sampling methods for varying sampling
rate (selectivity 5%, Impression)

Exp-V: Space cost under the same accuracy requirement. We evaluate the space
cost needed to achieve the same accuracy for different samplers. Since
Priority and Optimal GSW perform similarly, we focus on Optimal GSW and
compare it with Arithmetic compressed GSW. We fix the sample size of
Arithmetic compressed GSW (from $0.02\%$ to $1\%$). And for each measure, we
choose the size of an Optimal GSW sample so that it gives approximately the
same aggregation error as Arithmetic compressed GSW does. In Figure 15(a), we
report the total size of the four Optimal GSW samples (the portions for
different measures are stacked and labeled with different colors), and the
size of the Arithmetic compressed GSW sample; $x$-axis is the max aggregation
error in Arithmetic compressed GSW with the parameter $\Delta$ in brackets. It
shows that, under the same accuracy requirement, the total size of Optimal GSW
samples is around 1.8 times of the size of Arithmetic compressed GSW. Figure
15(b) reports forecast errors (of ARIMA) using the samples chosen above on
different measures. Due to the way how we choose the sizes of Optimal GSW
samples, Optimal GSW and Arithmetic compressed GSW give very close forecast
errors.

(a) Optimal v.s. Compressed

(b) Forecast error (selectivity 5%)

Figure 15. Space cost for given accuracy requirement

##  7\. Related Work

Approximate query processing and other samplers. An orthogonal line of work is
approximate query processing (AQP), which has been studied extensively during
last few decades. One can refer to (sigmod:ChaudhuriDK17, ) for a
comprehensive review. There are two major lines of AQP techniques. i) Online
aggregation (sigmod:Hellerstein97, ; nsdi:CondieCAHES10, ; pvldb:PansareBJC11,
) and online sampling-based AQP (sigmod:KandulaSVOGCD16, ) either assume that
the data is randomly ordered, or need to draw random rows from the data table
via random I/O accesses which is inefficient in our setting. ii) Offline
sampling-based AQP draws offline samples before queries come: some are based
on historical workloads (sigmod:Acharya00, ; VLDB:ganti00, ;
sigmod:Chaudhuri01, ; TODS:Chaudhuri07, ; cidr:OlstonBEJR09, ;
cidr:SidirourgosKB11, ; Eurosys:Agarwal13, ), and some are workload-
independent (sigmod:Acharya99, ; sigmod:Acharya00, ; ICDE:Chaudhuri01, ;
sigmod:BabcockCD03, ; sigmod:DingHCC016, ; sigmod:ParkMSW18, ).

The most relevant part in the line of AQP techniques are the samplers proposed
to estimate aggregations. We have reviewed such samplers at the beginning of
Section 4. There are other samplers such as universe (hashed) sampling and
stratified sampling introduced in AQP systems (Eurosys:Agarwal13, ;
sigmod:KandulaSVOGCD16, ; sigmod:ParkMSW18, ). These samplers were proposed to
handle orthogonal aspects such as missing groups in $\tt GroupBy$ and joins.
They can also be used in our system if we want to extend the task class and
data schema we want to support.

Precomputing aggregations. Another choice is to precompute aggregations or
summaries using techniques such as view materialization (VLDB:Agrawal00, ),
datacube (Springer:Gray1997data, ), histograms (PODS:Gilbert01, ), and
wavelets (sigmod:VitterW99, ; VLDB:Chakrabarti01, ; PODS:Gilbert01, ). These
techniques either have too large space overhead (super-linear) to be
applicable for enterprise-scale high-dim datasets, or cannot support complex
constraints in our forecasting tasks.

Aggregation-forecasting analytics. (sigmod:AgarwalCLSV10, ) studies a very
similar problem of forecasting multi-dimensional time series. It considers
capturing correlations across the high-dimensional space using either Bayesian
models or uniform sampling, and shows that the one based on uniform sampling
(which is also evaluated in our experimentes) offers much better forecast
accuracy. Another relevant line of work is about aggregation-disaggregation
techniques (jof:TobiasZ00, ; book:WestH97, ) in the forecasting literature.
These techniques share some similarity with the Bayesian model in
(sigmod:AgarwalCLSV10, ), but focus on one-time offline analysis with smaller
scale and lower dimensional data.

For multi-dimensional time series, there are works about how to conduct fast
similarity search (dasfaa:GongXHCLH15, ), but they are less relevant here.

##  8\. Conclusions

We introduce a real-time forecasting system FlashP for analyzing time-series
relations. There are two core technical contributions, which enable FlashP to
handle forecasting tasks on enterprise-scale high-dim time series
interactively. First, we analyze how sampling and estimated aggregations on
time-series relations enables interactive and accurate predictions. Second, we
propose GSW sampling, a new sampling scheme; we establish a connection from
the difference between GSW’s sampling weights and measure values to the
accuracy of aggregation estimations. Using this scheme, we can effectively
reduce the space (memory) consumption of weighted samples by letting multiple
measures share one GSW sample with the means of these measures’ values as
sampling weights.

## References

  * (1) https://www.alibabacloud.com.  [Online; accessed 1/15/2021].
  * (2) https://pypi.org/project/pmdarima/.  [Online; accessed 1/15/2021].
  * (3) https://www.statsmodels.org/stable/generated/statsmodels.tsa.x13.x13_arima_analysis.html.  [Online; accessed 1/15/2021].
  * (4) https://keras.io/.  [Online; accessed 1/15/2021].
  * (5) S. Acharya, P. B. Gibbons, and V. Poosala.  Congressional samples for approximate answering of group-by queries.  In SIGMOD, pages 487–498, 2000.
  * (6) S. Acharya, P. B. Gibbons, V. Poosala, and S. Ramaswamy.  The aqua approximate query answering system.  In SIGMOD, pages 574–576, 1999.
  * (7) D. Agarwal, D. Chen, L. Lin, J. Shanmugasundaram, and E. Vee.  Forecasting high-dimensional data.  In SIGMOD, pages 1003–1012, 2010.
  * (8) S. Agarwal, B. Mozafari, A. Panda, H. Milner, S. Madden, and I. Stoica.  Blinkdb: queries with bounded errors and bounded response times on very large data.  In Eurosys, pages 29–42, 2013.
  * (9) S. Agrawal, S. Chaudhuri, and V. R. Narasayya.  Automated selection of materialized views and indexes in sql databases.  In VLDB, pages 496–505, 2000.
  * (10) N. Alon, N. G. Duffield, C. Lund, and M. Thorup.  Estimating arbitrary subset sums with few probes.  In PODS, pages 317–325, 2005.
  * (11) B. Babcock, S. Chaudhuri, and G. Das.  Dynamic sample selection for approximate query processing.  In SIGMOD, pages 539–550, 2003.
  * (12) K. Chakrabarti, M. Garofalakis, R. Rastogi, and K. Shim.  Approximate query processing using wavelets.  VLDBJ, 10(2-3):199–223, 2001.
  * (13) S. Chaudhuri, G. Das, M. Datar, R. Motwani, and V. Narasayya.  Overcoming limitations of sampling for aggregation queries.  In ICDE, pages 534–542, 2001.
  * (14) S. Chaudhuri, G. Das, and V. Narasayya.  A robust, optimization-based approach for approximate answering of aggregate queries.  In SIGMOD, pages 295–306, 2001.
  * (15) S. Chaudhuri, G. Das, and V. Narasayya.  Optimized stratified sampling for approximate query processing.  TODS, 32(2):9, 2007.
  * (16) S. Chaudhuri, B. Ding, and S. Kandula.  Approximate query processing: No silver bullet.  In SIGMOD, pages 511–519, 2017.
  * (17) S. Chaudhuri, R. Motwani, and V. Narasayya.  On random sampling over joins.  In SIGMOD, pages 263–274, 1999.
  * (18) T. Condie, N. Conway, P. Alvaro, J. M. Hellerstein, K. Elmeleegy, and R. Sears.  Mapreduce online.  In NSDI, pages 313–328, 2010.
  * (19) B. Ding, S. Huang, S. Chaudhuri, K. Chakrabarti, and C. Wang.  Sample + seek: Approximating aggregates with distribution precision guarantee.  In SIGMOD, pages 679–694, 2016.
  * (20) N. G. Duffield, C. Lund, and M. Thorup.  Learn more, sample less: control of volume and variance in network measurement.  IEEE Trans. Information Theory, 51(5):1756–1775, 2005.
  * (21) N. G. Duffield, C. Lund, and M. Thorup.  Priority sampling for estimation of arbitrary subset sums.  J. ACM, 54(6):32, 2007.
  * (22) V. Ganti, M.-L. Lee, and R. Ramakrishnan.  Icicles: Self-tuning samples for approximate query answering.  In VLDB, page 187, 2000.
  * (23) A. C. Gilbert, Y. Kotidis, S. Muthukrishnan, and M. J. Strauss.  Optimal and approximate computation of summary statistics for range aggregates.  In PODS, pages 227–236, 2001.
  * (24) X. Gong, Y. Xiong, W. Huang, L. Chen, Q. Lu, and Y. Hu.  Fast similarity search of multi-dimensional time series via segment rotation.  In DASFAA, pages 108–124, 2015.
  * (25) J. Gray, S. Chaudhuri, A. Bosworth, A. Layman, D. Reichart, M. Venkatrao, F. Pellow, and H. Pirahesh.  Data cube: A relational aggregation operator generalizing group-by, cross-tab, and sub-totals.  Data Min. Knowl. Discov., 1(1):29–53, 1997.
  * (26) J. D. Hamilton.  Time Series Analysis.  Princeton University Press, 1994.
  * (27) S. Har-peled.  Geometric Approximation Algorithms.  American Mathematical Society, 2011.
  * (28) J. M. Hellerstein, P. J. Haas, and H. J. Wang.  Online aggregation.  SIGMOD, pages 171–182, 1997.
  * (29) S. Hochreiter and J. Schmidhuber.  Long short-term memory.  Neural Comput., 9(8):1735–1780, 1997.
  * (30) D. G. Horvitz and D. J. Thompson.  A generalization of sampling without replacement from a finite universe.  Journal of the American Statistical Association, 47(260):663–685, 1952.
  * (31) S. Kandula, A. Shanbhag, A. Vitorovic, M. Olma, R. Grandl, S. Chaudhuri, and B. Ding.  Quickr: Lazily approximating complex adhoc queries in bigdata clusters.  In SIGMOD, pages 631–646, 2016.
  * (32) C. Olston, E. Bortnikov, K. Elmeleegy, F. Junqueira, and B. Reed.  Interactive analysis of web-scale data.  In CIDR, 2009.
  * (33) N. Pansare, V. R. Borkar, C. Jermaine, and T. Condie.  Online aggregation for large mapreduce jobs.  PVLDB, 4(11):1135–1145, 2011.
  * (34) Y. Park, B. Mozafari, J. Sorenson, and J. Wang.  Verdictdb: Universalizing approximate query processing.  In SIGMOD, pages 1461–1476, 2018.
  * (35) J. Schmidhuber, D. Wierstra, and F. J. Gomez.  Evolino: Hybrid neuroevolution/optimal linear search for sequence learning.  In IJCAI, pages 853–858, 2005.
  * (36) L. Sidirourgos, M. L. Kersten, and P. A. Boncz.  Sciborq: Scientific data management with bounds on runtime and quality.  In CIDR, pages 296–301, 2011.
  * (37) M. Szegedy.  The DLT priority sampling is essentially optimal.  In STOC, pages 150–158, 2006.
  * (38) J. Tobias and A. Zellner.  A note on aggregation,disaggregation and forecasting performance.  Journal of Forecasting, 19(5):457–469, 2000.
  * (39) J. S. Vitter and M. Wang.  Approximate computation of multidimensional aggregates of sparse data using wavelets.  In SIGMOD, pages 193–204, 1999.
  * (40) M. West and J. Harrison.  Bayesian Forecasting and Dynamic Models.  Springer-Verlag, 1997.

##  Appendix A Missing Details

###  A.1. Proof of Proposition 1

We can write the noisy version of ${\rm ARMA(1,1)}$ as

$\hat{M}_{t}=\alpha_{1}\hat{M}_{t-1}+(u_{t}+\beta_{1}u_{t-1})+(\epsilon_{t}-\alpha_{1}\epsilon_{t-1}).$

Let ${\bf Var}[u_{t}]=\sigma^{2}_{u}$ and ${\bf
Var}[\epsilon_{t}]=\sigma^{2}_{\epsilon}$. We have

$\displaystyle{\bf E}[u_{t}\hat{M}_{t}]$ $\displaystyle={\bf
E}[u_{t}^{2}]+{\bf
E}[u_{t}(\alpha_{1}\hat{M}_{t-1}+\beta_{1}u_{t-1}+\epsilon_{t}-\alpha_{1}\epsilon_{t-1})]$
$\displaystyle=\sigma^{2}_{u}+0=\sigma^{2}_{u},$

from the fact that $u_{t}$ is dependent only on $\hat{M}_{t}$. Similarly,

${\bf E}[\epsilon_{t}\hat{M}_{t}]=\sigma_{\epsilon}^{2}.$

Taking the variance of both sides of the model, we have

$\displaystyle{\bf Var}[\hat{M}_{t}]=$ $\displaystyle{\bf
Var}[u_{t}+\epsilon_{t}]+{\bf
Var}[\alpha_{1}\hat{M}_{t-1}+\beta_{1}u_{t-1}-\alpha_{1}\epsilon_{t-1}]$
$\displaystyle=$
$\displaystyle\sigma^{2}_{u}+\sigma^{2}_{\epsilon}+\alpha_{1}^{2}{\bf
Var}[\hat{M}_{t-1}]+\beta_{1}^{2}\sigma^{2}_{u}+\alpha_{1}^{2}\sigma^{2}_{\epsilon}$
$\displaystyle+2\alpha_{1}\beta_{1}{\bf
Cov}[\hat{M}_{t-1},u_{t-1}]-2\alpha_{1}^{2}{\bf
Cov}[\hat{M}_{t-1},\epsilon_{t-1}]$ $\displaystyle-2\beta_{1}\alpha_{1}{\bf
Cov}[u_{t-1},\epsilon_{t-1}].$

Since $u_{t}$ and $\epsilon_{t}$ are independent, ${\bf
Cov}[\hat{M}_{t},u_{t}]={\bf E}[u_{t}\hat{M}_{t}]$, and ${\bf
Cov}[\hat{M}_{t},\epsilon_{t}]={\bf E}[\epsilon_{t}\hat{M}_{t}]$, we have

${\bf Var}[\hat{M}_{t}]=\alpha_{1}^{2}{\bf
Var}[\hat{M}_{t-1}]+(1+2\alpha_{1}\beta_{1}+\beta_{1}^{2})\sigma_{u}^{2}+(1-\alpha_{1}^{2})\sigma_{\epsilon}^{2}.$

Assuming $\hat{M}_{t}$ is weakly stationary (${\bf Var}[\hat{M}_{t}]={\bf
Var}[\hat{M}_{t-1}]$),

(11) $\displaystyle{\bf Var}[\hat{M}_{t}]$
$\displaystyle=\frac{(1+2\alpha_{1}\beta_{1}+\beta_{1}^{2})\sigma_{u}^{2}+(1-\alpha_{1}^{2})\sigma_{\epsilon}^{2}}{1-\alpha_{1}^{2}}$
$\displaystyle=\frac{1+2\alpha_{1}\beta_{1}+\beta_{1}^{2}}{1-\alpha_{1}^{2}}\cdot\sigma_{u}^{2}+\sigma_{\epsilon}^{2}$

Therefore, the variances of forecasts are linear combinations of
$\sigma_{u}^{2}$ and $\sigma_{\epsilon}^{2}$, with $\sigma_{u}^{2}$ decided by
the underlying data distribution, and $\sigma_{\epsilon}^{2}$ decided by
sampling scheme and sampling rate. $\Box$

###  A.2. Proof of Theorem 3

We can first calculate the variance of each ${\hat{m}_{i}}$ in (6):

${\bf
Var}\left[\hat{m}_{i}\right]=\frac{m_{i}^{2}(\Delta+w_{i})^{2}}{w_{i}^{2}}\cdot\frac{w_{i}}{\Delta+w_{i}}\cdot\frac{\Delta}{\Delta+w_{i}}=\frac{\Delta
m_{i}^{2}}{w_{i}}.$

Since random variables $\hat{m}_{i}$’s are independent, we have

(12) ${\bf Var}\left[\hat{M}\right]=\sum_{i=1}^{n}{\bf
Var}\left[\hat{m}_{i}\right]=\sum_{i=1}^{n}\frac{\Delta m_{i}^{2}}{w_{i}}.$

Define $W=\sum_{i=1}^{n}w_{i}$. If ${\bf w}$ is
$(\underline{\theta},\overline{\theta})$-consistent with ${\bf m}$,

${\bf Var}\left[\hat{M}\right]=\sum_{i=1}^{n}\frac{\Delta
m_{i}^{2}}{w_{i}}\leq\sum_{i=1}^{n}\frac{\Delta
m_{i}^{2}}{m_{i}/\overline{\theta}}=\overline{\theta}\Delta M.$

From (6) and $(\underline{\theta},\overline{\theta})$-consistency, the
expected sample size is

(13) $\\!{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]=\sum_{i=1}^{n}\frac{w_{i}}{\Delta+w_{i}}\leq\frac{\sum_{i=1}^{n}w_{i}}{\Delta}\leq\frac{\sum_{i=1}^{n}m_{i}}{\underline{\theta}\Delta}=\frac{M}{\underline{\theta}\Delta}.$

From the fact that $\hat{M}$ is unbiased and the above two inequalities,

(14) $\displaystyle{\bf E}\left[\left(\frac{{\hat{M}}-M}{M}\right)^{2}\right]$
$\displaystyle\leq\frac{\overline{\theta}\Delta
M}{M^{2}}=({\overline{\theta}}/{\underline{\theta}})\cdot\frac{\underline{\theta}\Delta}{M}\leq\frac{\overline{\theta}/\underline{\theta}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}.$

We can conclude the proof from the above inequality. $\Box$

###  A.3. Finding the Optimal Weight

We can choose $\bf w$ to minimize the variance of estimation in (12), while
the expected sample size in (13) is bounded. That is, solving

$\displaystyle\min_{\bf w}{\bf
Var}\left[\hat{M}\right]=\sum_{i=1}^{n}\frac{\Delta m_{i}^{2}}{w_{i}}$
$\displaystyle\hbox{s.t.}~{}{\bf
E}\left[|S_{\Delta}|\right]=\sum_{i=1}^{n}\frac{w_{i}}{\Delta+w_{i}}\leq B.$

Using the method of Lagrange multipliers, we can solve the above problem and
obtain the optimal solution $w^{*}_{i}$ satisfying:

$-\frac{\Delta
m_{i}^{2}}{{w^{*}_{i}}^{2}}+\frac{\lambda\Delta}{(\Delta+w^{*}_{i})^{2}}=0,$

where $\lambda$ is the Lagrange multiplier. If we consider a more restrive
constraint instead

$\sum_{i=1}^{n}\frac{w_{i}}{\Delta}\leq B,$

we have

$-\frac{\Delta m_{i}^{2}}{{w^{*}_{i}}^{2}}+\frac{\lambda}{\Delta}=0,$

and thus $w^{*}_{i}=\Delta m_{i}/\sqrt{\lambda}$.

###  A.4. Proof of Corollary 5

If we use ${\bf w}^{\times}$ as the sampling weights for a particular measure
${\bf m}^{(p)}$, for any row $i\in[n]$, we have

$\frac{{m}^{(p)}_{i}}{w^{\times}_{i}}=\frac{\left(\prod_{j=1}^{k}{m}^{(p)}_{i}\right)^{1/k}}{\left(\prod_{j=1}^{k}{m}^{(j)}_{i}\right)^{1/k}}=\prod_{j:~{}j\neq
p}\left(\frac{{m}^{(p)}_{i}}{{m}^{(j)}_{i}}\right)^{1/k},$

and thus,

$\max_{i}\frac{{m}^{(p)}_{i}}{w^{\times}_{i}}=\prod_{j:~{}j\neq
p}\left(\max_{i}\frac{{m}^{(p)}_{i}}{{m}^{(j)}_{i}}\right)^{1/k}=\prod_{j:~{}j\neq
p}\overline{\rho}^{1/k}_{p,j}\triangleq\overline{\theta};$

similarly,

$\min_{i}\frac{{m}^{(p)}_{i}}{w^{\times}_{i}}=\prod_{j:~{}j\neq
p}\underline{\rho}^{1/k}_{p,j}\triangleq\underline{\theta}.$

Therefore, ${\bf w}^{\times}$ is
$(\underline{\theta},\overline{\theta})$-consistent with ${\bf m}^{(p)}$,
where $\underline{\theta}$ and $\overline{\theta}$ are decided by
$\underline{\rho}_{p,q}$’s and $\overline{\rho}_{p,q}$’s as above,
respectively. From Theorem 3, we can derive the following error bound:

${\rm RE}({\hat{M}}^{(p)})\\!\leq\\!{\rm
RSTD}({\hat{M}}^{(p)})\\!\leq\sqrt{\frac{\prod_{j:~{}j\neq
p}\rho_{p,j}^{1/k}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}\leq\sqrt{\frac{(\rho^{1/k})^{k-1}}{{\bf
E}\left[|{\mathcal{S}}_{\Delta}|\right]}}$

which concludes the proof. $\Box$

###  A.5. Proof of Proposition 7

Since ${\bf w}$ $=$ $[w_{i}]$ are
$(\underline{\theta},\overline{\theta})$-consistent with measure ${\bf
m}=[{m}_{i}]$,

${m}^{\prime}_{i}=\frac{{m}_{i}}{\sum_{j}{m}_{j}}\leq\frac{\overline{\theta}w_{i}}{\sum_{j}\underline{\theta}w_{j}}=\frac{\overline{\theta}}{\underline{\theta}}\cdot
w^{\prime}_{i}=\theta\cdot w^{\prime}_{i},$

where $\theta=\overline{\theta}/\underline{\theta}$. Similarly,

${m}^{\prime}_{i}\geq\frac{1}{\theta}\cdot w^{\prime}_{i}.$

Putting them together, we have

$|{m}^{\prime}_{i}-w^{\prime}_{i}|\leq\max\\{(\theta-1),(1-1/\theta)\\}\cdot
w^{\prime}_{i}.$

Since $\theta\geq 1$, we have $\theta-1\geq 1-1/\theta$, and thus

$|{m}^{\prime}_{i}-w^{\prime}_{i}|\leq(\theta-1)\cdot w^{\prime}_{i}.$

Summing them up,

$\|{\bf m}^{\prime}-{\bf
w}^{\prime}\|_{1}=\sum_{i}|{m}^{\prime}_{i}-w^{\prime}_{i}|\leq(\theta-1)\cdot\sum_{i}w^{\prime}_{i}=(\theta-1).$

$\Box$