text:
\section{Missing Details}
\label{app:details}

\subsection{Proof of Proposition~\ref{prop:arma11}}
\label{prop:arma11:proof}

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\xspace_t - \alpha_1 \epsilon\xspace_{t-1}).
\]
Let $\vrinline{u_t} = \sigma^2_u$ and $\vrinline{\epsilon\xspace_t} = \sigma^2_\epsilon\xspace$. We have
\begin{align*}
\epinline{u_t \hat M_t} & = \epinline{u_t^2} + \epinline{u_t(\alpha_1 \hat M_{t-1}  + \beta_1 u_{t-1} + \epsilon\xspace_t - \alpha_1 \epsilon\xspace_{t-1})}
\\
& = \sigma^2_u + 0 = \sigma^2_u,
\end{align*}
from the fact that $u_t$ is dependent only on $\hat M_t$. Similarly,
\[
\epinline{\epsilon\xspace_t \hat M_t} = \sigma_\epsilon\xspace^2.
\]
Taking the variance of both sides of the model, we have
\begin{align*}
\vrinline{\hat M_t} = & \vrinline{u_t + \epsilon\xspace_t} + \vrinline{\alpha_1 \hat M_{t-1} + \beta_1 u_{t-1} - \alpha_1 \epsilon\xspace_{t-1}}
\\
= & \sigma^2_u + \sigma^2_\epsilon\xspace + \alpha_1^2 \vrinline{\hat M_{t-1}} + \beta_1^2 \sigma^2_u + \alpha_1^2 \sigma^2_\epsilon\xspace
\\
& + 2\alpha_1\beta_1\cvinline{\hat M_{t-1}, u_{t-1}} - 2\alpha_1^2\cvinline{\hat M_{t-1}, \epsilon\xspace_{t-1}}
\\
& - 2\beta_1\alpha_1\cvinline{u_{t-1}, \epsilon\xspace_{t-1}}.
\end{align*}
Since $u_{t}$ and $\epsilon\xspace_{t}$ are independent, $\cvinline{\hat M_{t}, u_{t}} = \epinline{u_t\hat M_t}$, and $\cvinline{\hat M_{t}, \epsilon\xspace_{t}} = \epinline{\epsilon\xspace_t \hat M_t}$, we have
\[
\vrinline{\hat M_t} = \alpha_1^2 \vrinline{\hat M_{t-1}}  + (1 + 2\alpha_1\beta_1 + \beta_1^2) \sigma_u^2 + (1-\alpha_1^2) \sigma_\epsilon\xspace^2.
\]
Assuming $\hat M_t$ is weakly stationary ($\vrinline{\hat M_t} = \vrinline{\hat M_{t-1}}$),
\begin{align} \label{equ:arma:var}
\vrinline{\hat M_t} & = \frac{(1+2\alpha_1\beta_1+\beta_1^2) \sigma_u^2 + (1-\alpha_1^2) \sigma_\epsilon\xspace^2}{1 - \alpha_1^2}
\\
& = \frac{1+2\alpha_1\beta_1+\beta_1^2}{1 - \alpha_1^2} \cdot \sigma_u^2 + \sigma_\epsilon\xspace^2 \nonumber
\end{align}
Therefore, the variances of forecasts are linear combinations of $\sigma_u^2$ and $\sigma_\epsilon\xspace^2$, with $\sigma_u^2$ decided by the underlying data distribution, and $\sigma_\epsilon\xspace^2$ decided by sampling scheme and sampling rate. {\hspace*{\fill}$\Box$\par}

\subsection{Proof of Theorem~\ref{thm:agg}}
\label{thm:agg:proof}
We can first calculate the variance of each ${\hat m_i}$ in \eqref{equ:gbin:variable:probability}:
\[
\vr{\hat m_i} = \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
\begin{equation} \label{equ:gbin:var}
\vr{\hat{M}} = \sum_{i=1}^n \vr{\hat m_i} = \sum_{i=1}^n\frac{\Delta m_i^2}{w_i}.
\end{equation}

Define $W = \sum_{i=1}^n w_i$. If ${\bf w}$ is $(\underline{\theta}\xspace, \overline{\theta}\xspace)$-consistent with ${\bf m}\xspace$,
\[
\vr{\hat{M}} = \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}\xspace} = \overline{\theta}\xspace \Delta M.
\]

From \eqref{equ:gbin:variable:probability} and $(\underline{\theta}\xspace, \overline{\theta}\xspace)$-consistency, the expected sample size is
\begin{equation} \label{equ:gbin:samplesize}
\!\ep{|{\mathcal S}_\Delta|} = \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}\xspace \Delta} = \frac{M}{\underline{\theta}\xspace\Delta}.
\end{equation}

From the fact that $\hat M$ is unbiased and the above two inequalities,
\begin{align}
\ep{\left(\frac{{\hat M}-M}{M}\right)^2} & \leq \frac{\overline{\theta}\xspace \Delta M}{M^2} = ({\overline{\theta}\xspace}/{\underline{\theta}\xspace}) \cdot \frac{\underline{\theta}\xspace\Delta}{M} \leq \frac{\overline{\theta}\xspace/\underline{\theta}\xspace}{\ep{|{\mathcal S}_{\Delta}|}}. \label{equ:gbin:re}
\end{align}

We can conclude the proof from the above inequality. {\hspace*{\fill}$\Box$\par}

\subsection{Finding the Optimal Weight}
\label{app:optweight}

We can choose $\bf w$ to minimize the variance of estimation in \eqref{equ:gbin:var}, while the expected sample size in \eqref{equ:gbin:samplesize} is bounded. That is, solving
\begin{align*}
& \min_{\bf w} \vr{\hat M} = \sum_{i=1}^n\frac{\Delta m_i^2}{w_i}
\\
& \hbox{s.t.} ~ \ep{|S_\Delta|} = \sum_{i=1}^n \frac{w_i}{\Delta + w_i} \leq B.
\end{align*}
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}$.

\subsection{Proof of Corollary~\ref{cor:agg:geo}}
\label{cor:agg:geo:proof}

If we use ${\bf w}^\times$ as the sampling weights for a particular measure ${\bf m}\xspace\supind{p}$, for any row $i \in [n]$, we have
\[
\frac{{m}\xspace\supind{p}_i}{w^\times_i} = \frac{\left(\prod_{j=1}^k {m}\xspace\supind{p}_{i}\right)^{1/k}}{\left(\prod_{j=1}^k {m}\xspace\supind{j}_{i}\right)^{1/k}}= \prod_{j:~j \neq p} \left(\frac{{m}\xspace\supind{p}_i}{{m}\xspace\supind{j}_i}\right)^{1/k},
\]
and thus,
\[
\max_i \frac{{m}\xspace\supind{p}_i}{w^\times_i} = \prod_{j:~j \neq p} \left(\max_i \frac{{m}\xspace\supind{p}_i}{{m}\xspace\supind{j}_i}\right)^{1/k} = \prod_{j:~j \neq p} \overline\rho^{1/k}_{p,j} \triangleq \overline{\theta}\xspace;
\]
similarly,
\[
\min_i \frac{{m}\xspace\supind{p}_i}{w^\times_i} = \prod_{j:~j \neq p} \underline\rho^{1/k}_{p,j} \triangleq \underline{\theta}\xspace.
\]
Therefore, ${\bf w}^\times$ is $(\underline{\theta}\xspace,\overline{\theta}\xspace)$-consistent with ${\bf m}\xspace\supind{p}$, where $\underline{\theta}\xspace$ and $\overline{\theta}\xspace$ are decided by $\underline\rho_{p,q}$'s and $\overline\rho_{p,q}$'s as above, respectively.
From Theorem~\ref{thm:agg}, we can derive the following error bound:
\[
{\rm RE}({\hat M}\supind{p}) \! \leq \! {\rm RSTD}({\hat M}\supind{p}) \! \leq \sqrt{\frac{\prod_{j: ~ j \neq p}\rho_{p,j}^{1/k}}{\ep{|{\mathcal S}_{\Delta}|}}} \leq \sqrt{\frac{(\rho^{1/k})^{k-1}}{\ep{|{\mathcal S}_{\Delta}|}}}
\]
which concludes the proof. {\hspace*{\fill}$\Box$\par}

\subsection{Proof of Proposition~\ref{prop:consistencyl2}}
\label{prop:consistencyl2:proof}

Since ${\bf w}$ $=$ $[w_i]$ are $(\underline{\theta}\xspace,\overline{\theta}\xspace)$-consistent with measure ${\bf m}\xspace = [{m}\xspace_i]$,
\[
{m}\xspace'_i = \frac{{m}\xspace_i}{\sum_j {m}\xspace_j} \leq \frac{\overline{\theta}\xspace w_i}{\sum_j \underline{\theta}\xspace w_j} = \frac{\overline{\theta}\xspace}{\underline{\theta}\xspace} \cdot w'_i = \theta \cdot w'_i,
\]
where $\theta = \overline{\theta}\xspace/\underline{\theta}\xspace$. Similarly,
\[
{m}\xspace'_i \geq \frac{1}{\theta} \cdot w'_i.
\]
Putting them together, we have
\[
|{m}\xspace'_i - w'_i| \leq \max\{(\theta - 1), (1 - 1/\theta)\} \cdot w'_i.
\]
Since $\theta \geq 1$, we have $\theta - 1 \geq 1 - 1/\theta$, and thus
\[
|{m}\xspace'_i - w'_i| \leq (\theta - 1) \cdot w'_i.
\]
Summing them up,
\[
\|{\bf m}\xspace' - {\bf w}'\|_1 = \sum_i |{m}\xspace'_i - w'_i| \leq (\theta - 1) \cdot \sum_i w'_i = (\theta - 1).
\]
\eop

\section{Deployment}
\label{sec:deploy}
\ifdefined
We have implemented and deployed {\sf {F\lowercase{lash}P}}\xspace in Alibaba's ads analytical platform.
Figure~\ref{fig:deploy} shows important components of {\sf {F\lowercase{lash}P}}\xspace
\fi

\stitle{Offline Sample Preprocessor} of {\sf {F\lowercase{lash}P}}\xspace is built on Alibaba's distributed data storage and analytics service, MaxCompute \cite{url:alibabacloud}. Time series of relations, partitioned by time, are stored in MaxCompute's data warehouse. Relations can be joined, {\em e.g.}\xspace, 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.

\stitle{Online Forecasting Service} first pulls Multi-layer Samples of Relations into Alibaba's in-memory OLAP engine, Hologres \cite{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 \eqref{equ:task:agg} 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 \cite{url:autoarima} (built on X-13ARIMA-SEATS \cite{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~\ref{fig:lstm:noisy}), which is implemented using Keras \cite{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.

\eat{
Discuss the real system architecture

In our system, there are four modules:
1.Offline sampling Module
2.Online aggregation Module
3.Online forecast module
4.Online application module

\subsection{Offline sampling Module}
Offline sampling module is built on \textit{ODPS} which is a distributed  compute system like \textit{Hive} in Alibaba.
There are two datasets including \textit{UserProfile} which with schema (user\_id, age, gender, ...) and \textit{UserAction} which with schema (date, user\_id, impression, click, favorite, cart)
Offline sampling module joins \textit{UserProfile} and \textit{UserAction} with
\textit{user\_id} in offline per day and use a sampler(e.g. Optimal Bernoulli Sampling) to take a multi-layer sample and load the sample dataset into the Online aggregation Module.

\subsection{Online aggregation Module}
Online aggregation Module is built on \textit{Hologres} which is an OLAP engine like \textit{Greenplum} in Alibaba. There are 30 servers in the cluster each one with 100 CPU and 500G memory. Online aggregation Module load the latest sample dataset every day and answer \textit{SUM} aggregation query which submit from the application server.

\subsection{Online forecast Module}
Online forecast Module is built on a \textit{python} server, which provide a \textit{auto\_arima} model to forecast. Online forecast Module receive a historical time-series data point and forecast period, and fit them into the \textit{auto\_arima} model. After \textit{auto\_arima} model provide the forecast result, Online forecast Module reply the answer to application module.

\subsection{Online Application Module}
Online Application Module process the query which submitted by the User via WEB-UI and transform the query to an Aggregation SQL and submit to the Online Aggregation Module. After receiving the historical result computed by the Online Aggregation Module, online application module fit the historical result into Online Forecast Module and wait the forecast result. After the result forecast, Online Application Module return it to the User.
}
\section{Experimental Evaluation}
\label{sec:exp}
We evaluate our system {\sf {F\lowercase{lash}P}}\xspace under the implementation and hardware specified in Section~\ref{sec:deploy},
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: {\bf Uniform} sampling, which is also used in \cite{sigmod:AgarwalCLSV10}, is the baseline; {\bf Priority}, the optimal weighted sampler \cite{jacm:DuffieldLT07}; our {\bf Optimal {\rm GSW}\xspace} and {\bf Arithmetic/Geometric compressed {\rm GSW}\xspace} introduced in Section~\ref{sec:sampler}. We also compare our sampling based methods with {\bf PIM} (Partwise Independence Model) \cite{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 {\em 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 {\em relative aggregation errors} (average of the 150 days), {\em relative forecast error}, and {\em 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).

\stitle{Exp-I: A summary of results.}
We first give a brief summary of experimental results. Table~\ref{fig:exp:methods} 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 {\rm GSW}\xspace and compressed {\rm GSW}\xspace in {\sf {F\lowercase{lash}P}}\xspace 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.

\begin{table}[t]
\small\centering
\caption{{A summary of results (0.1\% sample, Opt-GSW = Optimal GSW, C-GSW = Arithmetic compressed GSW)}}
\label{fig:exp:methods}
\vspace{-0.4cm}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 & Full & PIM & Uniform & Opt-GSW & C-GSW
\\ \hline \hline
$\sf Favorite$ & 0.105  & 0.695 & 0.248 & 0.131 & 0.196
\\ \hline
$\sf Impression$ & 0.140 & 0.374 & 0.147 & 0.142 & 0.144
\\ \hline
$\sf Click$ & 0.157 & 0.681 & 0.161 & 0.151 & 0.153
\\ \hline
$\sf Cart$ & 0.704 & 1.931 & 0.718 & 0.704 & 0.709
\\ \hline
\end{tabular}
\vspace{-0.4cm}
\end{table}

\begin{wrapfigure}{r}{0.18\textwidth}
 	\vspace{-0.6cm}
 	\begin{center}
 		\hspace{-0.9cm}\includegraphics[width=0.23\textwidth]{fig-new/vldb-fig5-5.eps}
 	\end{center}
 	\vspace{-0.4cm}
 	\caption{E2E response time (with ARIMA)}
 	\label{fig:exp:scalability:time:05}
 	\vspace{-0.4cm}
\end{wrapfigure}

\stitle{{Exp-II: Real-time response.}}
The end-to-end response time is reported in Figure~\ref{fig:exp:scalability:time:05}, 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.

\ifdefined
\stitle{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 {\sf Impression}. The average forecast errors, with one standard deviation, for varying sampling rate using our optimal {\rm GSW}\xspace are plotted in Figure~\ref{fig:exp:scalability}. 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.
\else
\stitle{Exp-III: Varying number of time stamps used in training.}
We consider different numbers of training data points used to fit the ARIMA and LSTM models. Please refer to the full version \cite{url:tr} for more details.
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.
\fi

\ifdefined
\begin{figure}[ht]
\vspace{-0.2cm}
\begin{minipage}[t]{0.479\linewidth}
\subfigure[ARIMA]{
\centering
\includegraphics[width=1\linewidth]{fig-new/vldb-fig3-fix-sample-500w_3-period-pv-acc.eps}
\label{fig:exp:scalability:error:05}
}
\end{minipage}
~~~
\begin{minipage}[t]{0.479\linewidth}
\subfigure[LSTM]{
\centering
\includegraphics[width=1\linewidth]{fig-paper-lstm/vldb-fig3-fix-sample-500w_3-period-pv-acc.eps}
\label{fig:exp:scalability:error:05:lstm}
}
\end{minipage}
\vspace{-0.5cm}
\caption{Number of time stamps used in training v.s. Forecast error (selectivity $5\%$, $\sf Impression$)}
\vspace{-0.3cm}
\label{fig:exp:scalability}
\end{figure}
\fi

\begin{figure}[ht]
\vspace{-0.2cm}
\subfigure[Selectivity 0.5\%]{
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-50w_3-150-col-aqp.eps}
\end{minipage}
}
\subfigure[Selectivity 5\%]{
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-500w_3-150-col-aqp.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Aggregation error of different sampling methods for varying selectivity and sampling rate ({\sf Favorite})}
\vspace{-0.2cm}
\label{fig:exp:aerror:favorite}
\end{figure}

\begin{figure}[t]
\vspace{-0.2cm}
\subfigure[ARIMA]{\label{fig:exp:ferror:favorite:arima}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-50w_3-150-col-acc.eps}
\end{minipage}
}
\subfigure[LSTM]{\label{fig:exp:ferror:favorite:lstm}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-paper-lstm/vldb-fig4-fix-sample-50w_3-150-col-acc-legend.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Forecast error of different sampling methods for varying sampling rate (selectivity 0.5\%, {\sf Favorite})}
\vspace{-0.4cm}
\label{fig:exp:ferror:favorite}
\end{figure}

\ifdefined
\begin{figure}[t]
\vspace{-0.2cm}
\subfigure[ARIMA]{\label{fig:exp:ferror:favorite:5:arima}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-500w_3-150-col-acc.eps}
\end{minipage}
}
\subfigure[LSTM]{\label{fig:exp:ferror:favorite:5:lstm}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-paper-lstm/vldb-fig4-fix-sample-500w_3-150-col-acc-legend.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Forecast error of different sampling methods for varying sampling rate (selectivity 5\%, {\sf Favorite})}
\vspace{-0.4cm}
\label{fig:exp:ferror:favorite:5}
\end{figure}
\fi

\stitle{Exp-IV: Varying sampling rate and selectivity.} We compare different samplers in {\sf {F\lowercase{lash}P}}\xspace, for varying sampling rate and selectivity. Note that for Arithmetic/Geometric compressed {\rm GSW}\xspace, one sample suffices for the relation; Priority are Optimal {\rm GSW}\xspace 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 {\rm GSW}\xspace for a fixed sampling rate.

\ifdefined
	For tasks on $\sf Favorite$ with selectivity $0.5\%$-$5\%$, Figure~\ref{fig:exp:aerror:favorite} reports aggregation errors, and Figures~\ref{fig:exp:ferror:favorite}-\ref{fig:exp:ferror:favorite:5} report forecast errors when using ARIMA and LSTM as forecasting models; the results on $\sf Impression$ are plotted in Figures~\ref{fig:exp:ferror:impression}-\ref{fig:exp:ferror:impression:5}.
\else
	For tasks on $\sf Favorite$ with selectivity $0.5\%$-$5\%$, Figure~\ref{fig:exp:aerror:favorite} reports aggregation errors, and Figure~\ref{fig:exp:ferror:favorite} reports forecast errors when using ARIMA and LSTM as forecasting models; the results on $\sf Impression$ are plotted in Figure~\ref{fig:exp:ferror:impression} (those with selectivity $5\%$ can be found in the full version \cite{url:tr}).
\fi
In terms of both aggregation errors and forecasting errors, Priority and Optimal {\rm GSW}\xspace are very close and better than the others (indeed, at the cost of storing four samples).
In some cases, Optimal {\rm GSW}\xspace 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 \cite{sigmod:Hellerstein97}. Arithmetic/Geometric compressed {\rm GSW}\xspace needs only one sample, too; they are better than Uniform, and get very close to Priority and Optimal {\rm GSW}\xspace 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~\ref{fig:exp:finterval:favorite:05} 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 {\rm GSW}\xspace give the narrowest ones. Figure~\ref{fig:exp:finterval:favorite:05:case} 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 ({\em e.g.}\xspace, when sampling rate = $1\%$ in Figures~\ref{fig:exp:aerror:favorite}-\ref{fig:exp:ferror:impression}), 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 ({\em e.g.}\xspace, $u_t$ in \eqref{equ:arma} for ARIMA). Forecast errors and forecast intervals have similar trends as aggregation errors. Both observations are consistent with our analytical results in Section~\ref{sec:rtforecast}.

\begin{figure}[t]
\vspace{-0.2cm}
\subfigure[Varying sampling rate]{
\label{fig:exp:finterval:favorite:05}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-50w_3-150-col-gap.eps}
\end{minipage}
}
\subfigure[Forecast intervals of a real query on 0.02\% sampling]{
\label{fig:exp:finterval:favorite:05:case}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig9-fix-sample-6492231-5000-col-case.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Forecast intervals (ARIMA) with different sampling methods for varying sampling rate (selectivity 0.5\%, {\sf Favorite})}
\vspace{-0.2cm}
\label{fig:exp:finterval:favorite}
\end{figure}

\begin{figure}[t]
\vspace{-0.2cm}
\subfigure[ARIMA]{\label{fig:exp:ferror:impression:arima}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-50w_3-150-pv-acc.eps}
\end{minipage}
}
\subfigure[LSTM]{\label{fig:exp:ferror:impression:lstm}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-paper-lstm/vldb-fig4-fix-sample-50w_3-150-pv-acc-legend.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Forecast error of different sampling methods for varying sampling rate (selectivity 0.5\%, {\sf Impression})}
\vspace{-0.2cm}
\label{fig:exp:ferror:impression}
\end{figure}

\ifdefined
\begin{figure}[t]
\vspace{-0.2cm}
\subfigure[ARIMA]{\label{fig:exp:ferror:impression:5:arima}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig4-fix-sample-500w_3-150-pv-acc.eps}
\end{minipage}
}
\subfigure[LSTM]{\label{fig:exp:ferror:impression:5:lstm}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-paper-lstm/vldb-fig4-fix-sample-500w_3-150-pv-acc-legend.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Forecast error of different sampling methods for varying sampling rate (selectivity 5\%, {\sf Impression})}
\vspace{-0.4cm}
\label{fig:exp:ferror:impression:5}
\end{figure}
\fi

\stitle{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 {\rm GSW}\xspace perform similarly, we focus on Optimal {\rm GSW}\xspace and compare it with Arithmetic compressed {\rm GSW}\xspace.
We fix the sample size of Arithmetic compressed {\rm GSW}\xspace (from $0.02\%$ to $1\%$). And for each measure, we choose the size of an Optimal {\rm GSW}\xspace sample so that it gives approximately the same aggregation error as Arithmetic compressed {\rm GSW}\xspace does. In Figure~\ref{fig:exp:fixacc:samplesize}, we report the total size of the four Optimal {\rm GSW}\xspace samples (the portions for different measures are stacked and labeled with different colors), and the size of the Arithmetic compressed {\rm GSW}\xspace sample; $x$-axis is the max aggregation error in Arithmetic compressed {\rm GSW}\xspace with the parameter $\Delta$ in brackets. It shows that, under the same accuracy requirement, the total size of Optimal {\rm GSW}\xspace samples is around 1.8 times of the size of Arithmetic compressed {\rm GSW}\xspace.
Figure~\ref{fig:exp:fixacc:ferror} reports forecast errors (of ARIMA) using the samples chosen above on different measures. Due to the way how we choose the sizes of Optimal {\rm GSW}\xspace samples, Optimal {\rm GSW}\xspace and Arithmetic compressed {\rm GSW}\xspace give very close forecast errors.

\begin{figure}[t]
\vspace{-0.2cm}
\subfigure[Optimal v.s. Compressed]{
\label{fig:exp:fixacc:samplesize}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig11-fix-aqp-sample.eps}
\end{minipage}
}
\subfigure[Forecast error {(selectivity 5\%)}]{
\label{fig:exp:fixacc:ferror}
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1.1\linewidth]{fig-new/vldb-fig13-fix-aqp-500w_3-acc.eps}
\end{minipage}
}
\vspace{-0.3cm}
\caption{Space cost for given accuracy requirement}
\vspace{-0.4cm}
\label{fig:exp:fixacc}
\end{figure}

\eat{
\subsection{Efficiency}
Experiments on speedup via AQP

\cfig1: fix the dataset and the sampling rate; (randomly) pick a few queries; x-axis is the selectivity ({\em i.e.}\xspace, number of users satisfying the predicates in the queries); y-axis is the (average) AQP time + time used by the prediction model (using, for example, stacked histograms)

\cfig2: try a public dataset and repeat \cfig1?

\cfig3: fix the dataset and the sampling rate; (randomly) pick a few queries; x-axis is the number of constraints in the query; y-axis is the (average) AQP time + time used by the prediction model (using, for example, stacked histograms)

\cfig4: try a public dataset and repeat \cfig3?

\cfig5: fix the dataset; (randomly) pick a few queries; x-axis is the sampling rate (from 0.001 to 0.1); y-axis is the (average) AQP time + time used by the prediction model (using, for example, stacked histograms)

\cfig6: try a public dataset and repeat \cfig5?

\subsection{Prediction Precision}
How AQP impacts prediction precision

For each experiment, try different datasets (for different ads), and for each dataset, pick a few queries

\cfig7: fix the dataset; x-axis is the sampling rate; y-axis is the predication error (group by days in the prediction)

\cfig8: fix the dataset; x-axis is the number of training time stamps; y-axis is the predication error (group by days in the prediction)

\cfig9: fix the dataset; x-axis is the selectivity of queries; y-axis is the predication error (group by days in the prediction)
}
\section{Experimental Evaluation}

\subsection{Efficiency}
Experiments on speedup via AQP

\cfig1: fix the dataset and the sampling rate; (randomly) pick a few queries; x-axis is the selectivity ({\em i.e.}\xspace, number of users satisfying the predicates in the queries); y-axis is the (average) AQP time + time used by the prediction model (using, for example, stacked histograms)
\begin{figure*}[ht]
\begin{minipage}[t]{0.33\linewidth}
\subfigure[Sampling rate 0.1\%]{
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig1-1000.eps}
}
\end{minipage}
\begin{minipage}[t]{0.33\linewidth}
\subfigure[Sampling rate 1\%]{
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig1-100.eps}
}
\end{minipage}
\begin{minipage}[t]{0.33\linewidth}
\subfigure[Sampling rate 10\%]{
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig1-10.eps}
}
\end{minipage}
\caption{**AQP and ARIMA Performance with different selectivity at fixed sampling rate}
\end{figure*}

\cfig2: try a public dataset and repeat \cfig1?

\cfig3: fix the dataset and the sampling rate; (randomly) pick a few queries; x-axis is the number of constraints in the query; y-axis is the (average) AQP time + time used by the prediction model (using, for example, stacked histograms)

\cfig4: try a public dataset and repeat \cfig3?

\cfig5: fix the dataset; (randomly) pick a few queries; x-axis is the sampling rate (from 0.001 to 0.1); y-axis is the (average) AQP time + time used by the prediction model (using, for example, stacked histograms)

\begin{figure*}[t]
\begin{minipage}[t]{0.33\linewidth}
\subfigure[Selectivity 0.5\%]{
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig5-500000.eps}
}
\end{minipage}
\begin{minipage}[t]{0.33\linewidth}
\subfigure[Selectivity 1\%]{
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig5-1000000.eps}
}
\end{minipage}
\begin{minipage}[t]{0.33\linewidth}
\subfigure[Selectivity 5\%]{
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig5-5000000.eps}
}
\end{minipage}
\caption{AQP and ARIMA Performance with different sampling rate at fixed selectivity (how many rows and days in total? ms?)}
\end{figure*}

\cfig6: try a public dataset and repeat \cfig5?

\subsection{Prediction Precision}
How AQP impacts prediction precision

For each experiment, try different datasets (for different ads), and for each dataset, pick a few queries

\cfig7: fix the dataset; x-axis is the sampling rate; y-axis is the predication error (group by days in the prediction)

\begin{figure*}[ht]
\subfigure[Selectivity 0.5\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig7-50w_3366106_2_180.eps}
\end{minipage}
}
\subfigure[Selectivity 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig7-100w_202210022_2_180.eps}
\end{minipage}
}
\subfigure[Selectivity 5\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig7-500w_718978045_2_180.eps}
\end{minipage}
}
\caption{**ARIMA predict error on impression at different sampling rate}
\end{figure*}

\cfig8: fix the dataset; x-axis is the number of training time stamps; y-axis is the predication error (group by days in the prediction)

\begin{figure*}[hb]
\subfigure[Selectivity 0.5\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig9-500w_718978045_2_100.eps}
\end{minipage}
}
\subfigure[Selectivity 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig9-100w_202210022_2_100.eps}
\end{minipage}
}
\subfigure[Selectivity 5\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig9-500w_718978045_2_100.eps}
\end{minipage}
}
\caption{ARIMA predict error on impression at different training time stamps (add sampling rate 0.1\%)}
\end{figure*}

\cfig9: fix the dataset; x-axis is the selectivity of queries; y-axis is the predication error (group by days in the prediction)

\begin{figure*}[ht]
\subfigure[Average prediction error]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig40-50w_129370505_3-30-pv.eps}
\end{minipage}
}
\subfigure[Average prediction interval]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig41-50w_129370505_3-30-pv.eps}
\end{minipage}
}
\caption{ARIMA predict error and interval at selectivity 0.5\% (how many queries?) on impression via different sampling methods (add sampling rate 0.005\%? how about AQP error?)}
\end{figure*}

\begin{figure*}[ht]
\subfigure[Average prediction error]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig40-50w_129370505_3-30-clk.eps}
\end{minipage}
}
\subfigure[Average prediction interval]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig41-50w_129370505_3-30-clk.eps}
\end{minipage}
}
\caption{ARIMA predict error and interval at selectivity 0.5\% on click via different sampling methods}
\end{figure*}

\begin{figure*}[ht]
\subfigure[Average prediction error]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig40-50w_129370505_3-30-col.eps}
\end{minipage}
}
\subfigure[Average prediction interval]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig41-50w_129370505_3-30-col.eps}
\end{minipage}
}
\caption{ARIMA predict error and interval at selectivity 0.5\% on favorite via different sampling methods}
\end{figure*}

\begin{figure*}[ht]
\subfigure[Average prediction error]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig40-50w_129370505_3-30-cart.eps}
\end{minipage}
}
\subfigure[Average prediction interval]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=3.3in, height=3in]{fig/vldb-fig41-50w_129370505_3-30-cart.eps}
\end{minipage}
}
\caption{ARIMA predict error and interval at selectivity 0.5\% on cart via different sampling methods}
\end{figure*}

\begin{figure*}[hb]
\subfigure[Sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-5000-pv.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-1000-pv.eps}
\end{minipage}
}
\subfigure[Sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-100-pv.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\% (why not 0.5\%? -- inconsistent with 0.5\% above) on impression via different sampling methods (so the predictions in different sampling methods are almost the same?)}
\end{figure*}

\begin{figure*}[hb]
\subfigure[Sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-5000-clk.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-1000-clk.eps}
\end{minipage}
}
\subfigure[Sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-100-clk.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\% on click via different sampling methods }
\end{figure*}

\begin{figure*}[hb]
\subfigure[Sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-5000-col.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-1000-col.eps}
\end{minipage}
}
\subfigure[Sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-100-col.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\%  on favorite via different sampling methods }
\end{figure*}

\begin{figure*}[hb]
\subfigure[Sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-5000-cart.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-1000-cart.eps}
\end{minipage}
}
\subfigure[Sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig42-9248224-100-cart.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\%  on cart via different sampling methods }
\end{figure*}

\begin{figure*}[hb]
\subfigure[Sampling rate (what is x-axis? AQP error?)]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig1-sample.eps}
\end{minipage}
}
\subfigure[AQP error bound (only one compressed sample, right?)]{
\begin{minipage}[t]{0.5\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig2-3-aqp.eps}
\end{minipage}
}
\caption{Sampling rate via different sampling methods with nearly aqp error bound (can we add uniform and geometric-mean here?)}
\end{figure*}

\begin{figure*}[ht]
\subfigure[Sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2in]{fig/vldb-fig43-500w_9248224_3-100-acc.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2in]{fig/vldb-fig43-500w_9248224_3-1000-acc.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2in]{fig/vldb-fig43-500w_9248224_3-5000-acc.eps}
\end{minipage}
}
\caption{ARIMA predict error at selectivity 5\%  via different sampling methods with nearly aqp error bound (can we add uniform and geometric-mean here? almost no change for different sampling rates? how can we guarantee that aqp error bounds are nearly the same given the same sampling rate?)}
\end{figure*}

\begin{figure*}[ht]
\subfigure[Sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2in]{fig/vldb-fig43-500w_9248224_3-100-gap.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2in]{fig/vldb-fig43-500w_9248224_3-1000-gap.eps}
\end{minipage}
}
\subfigure[Sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2in]{fig/vldb-fig43-500w_9248224_3-5000-gap.eps}
\end{minipage}
}
\caption{ARIMA predict interval at selectivity 5\% via different sampling methods with nearly aqp error bound (can we add uniform and geometric-mean here? almost no change for different sampling rates? how can we guarantee that aqp error bounds are nearly the same given the same sampling rate?)}
\end{figure*}

\begin{figure*}[hb]
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-5000-pv.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-1000-pv.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-100-pv.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval on impression at selectivity 5\%  via different sampling methods with nearly AQP error bound }
\end{figure*}

\begin{figure*}[hb]
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-5000-clk.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-1000-clk.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-100-clk.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\%  on click via different sampling methods with nearly AQP error bound }
\end{figure*}

\begin{figure*}[hb]
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-5000-col.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-1000-col.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-100-col.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\%  on favorite via different sampling methods with nearly AQP error bound }
\end{figure*}

\begin{figure*}[hb]
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.02\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-5000-cart.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 0.1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-1000-cart.eps}
\end{minipage}
}
\subfigure[AQP average error nearly Compressed-arithmetic-mean sampling with sampling rate 1\%]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=2.2in, height=2.0in]{fig/vldb-fig43-9248224-100-cart.eps}
\end{minipage}
}
\caption{ARIMA average predict error and interval at selectivity 5\%  on cart via different sampling methods with nearly AQP error bound }
\end{figure*}

\section{Introduction}
Large volumes of high-dimensional data are generated on eCommerce platforms every day, from data sources about, {\em e.g.}\xspace, 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.

\stitle{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 ({\em i.e.}\xspace, {\em 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, {\em e.g.}\xspace, {\sf Impression}, the number of times an advertisement can be showed to such customers, and {\sf 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 ({\em e.g.}\xspace, mobile/PC), machine-learned tags ({\em e.g.}\xspace, 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~\ref{fig:motivation}. A forecasting task, {\em e.g.}\xspace, the one in Figure~\ref{fig:execution:example}, 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.

\stitle{The task and solution phases.}
Our forecasting system {\sf {F\lowercase{lash}P}}\xspace is designed to process {\em real-time forecasting tasks}, in which we specify: i) a constraint specifying the portion of data to be focused on ({\em e.g.}\xspace, targeting customers ``${\sf Age}$ $<=$ $30$ $\code{AND}$ ${\sf Gender}$ $=$ ${\rm F}$'' in the example in Figure~\ref{fig:execution:example}); ii) the aggregated measure to be predicted ({\em e.g.}\xspace, ${\tt SUM}({\sf Impression})$); iii) historical data points to be used to train a forecasting model ({\em e.g.}\xspace, from $20200101$ to $20200331$).

We give an overview of the solution.
Consider the example in Figure~\ref{fig:execution:example}. 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 \code{GROUP} ~ \code{BY}}$) is the bottleneck. A natural idea is to utilize sampling ({\em e.g.}\xspace, \cite{jacm:DuffieldLT07}) or sample-based approximate query processing ({\em e.g.}\xspace, \cite{sigmod:DingHCC016,sigmod:ParkMSW18}) to estimate their answers.
Figure~\ref{fig:execution:result} 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).

\stitle{Contributions and organization.}
{\sf {F\lowercase{lash}P}}\xspace 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~\ref{sec:rtforecast}). Error bounds of estimating aggregations using samples are well studied {\em e.g.}\xspace, for priority sampling \cite{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 ({\em i.e.}\xspace, ${\rm ARMA}(1,1)$, defined later), and experimental results for more complicated models ({\em e.g.}\xspace, LSTM).

The second question is about the space budget needed for samples (Section~\ref{sec:sampler}).
Uniform sampling provides unbiased estimations for aggregations, but the estimation error is proportional to the range of a measure ($\max - \min$) \cite{sigmod:Hellerstein97}.
Weighted sampling schemes offer optimal estimations, with much better error bounds that are independent on the range \cite{stoc:Szegedy06}; however, as the sampling distribution is decided by the measure values \cite{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 {\rm GSW}\xspace ({\em \underline{G}eneralized \underline{S}moothed \underline{W}eighted}) 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 {\sf {F\lowercase{lash}P}}\xspace is implemented in Section~\ref{sec:deploy}. We report experimental results on real datasets in Section~\ref{sec:exp}, and discuss related work in Section~\ref{sec:related}.
\ifdefined
	All missing proofs are deferred to Appendix~\ref{app:details}.
\else
	All missing proofs are in the full version \cite{url:tr}.
\fi

\section{Solution Overview}
\label{sec:overview}

\stitle{Time-series data model.} The input to our forecasting pipeline is a {\em time series of relation} $T$, {\em i.e.}\xspace, a sequence of {\em 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 {\em dimensions}, each denoted by ${\bf a}\xspace$, which are used to filter data ({\em e.g.}\xspace, $\sf Age$ and $\sf Location$), and multiple {\em measures}, each denoted by ${\bf m}\xspace$, which we want to analyze and forecast ({\em e.g.}\xspace, $\sf Impression$ and $\sf ViewTime$). We use ${a}\xspace_{i,t}$ and ${m}\xspace_{i,t}$ to denote the values of a dimension and a measure, respectively, on each row; or, simply, ${a}\xspace_{i}$ and ${m}\xspace_{i}$, if the time stamp $t$ is clear from the context or not important.
The schema of $T$ is $({\bf a}\xspace\supind{1}, {\bf a}\xspace\supind{2}, \ldots, {\bf a}\xspace\supind{d_a}; {\bf m}\xspace\supind{1}, {\bf m}\xspace\supind{2}, \ldots, {\bf m}\xspace\supind{d_m}; t)$.

\stitle{Real-time forecasting task.} A {\em forecasting tasks} is specified as:
\begin{align}
& \code{FORECAST} ~ {\tt SUM}({\bf m}\xspace) ~ \code{FROM} ~ T ~ \code{WHERE} ~ {\sf C}\xspace ~ \code{USING} ~ (t_s, t_e) \label{equ:task}
\\
& \code{OPTION} ~ ({\tt MODEL} = {\tt 'model\_x'}, {\tt FORE\_PERIOD} = {\tt t\_future}) \nonumber
\end{align}

Here, ${\bf m}\xspace$ is the measure we want to forecast from data source $T$ on a given set of rows satisfying the constraint ${\sf C}\xspace$. {\sf {F\lowercase{lash}P}}\xspace allows ${\sf C}\xspace$ to be any logical expression on the dimension values of ${\bf a}\xspace\supind{1}, \ldots, {\bf a}\xspace\supind{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 $\code{OPTION}$ clause with the parameter $\tt MODEL$, and the number of future time stamps we want to predict with the parameter $\tt FORE\_PERIOD$ ({\em e.g.}\xspace, 7 days).
In most of our application scenarios, we care about $\tt SUM$ aggregation, {\em e.g.}\xspace, total number of impressions from a certain group users; {\sf {F\lowercase{lash}P}}\xspace can also support $\tt COUNT$ and $\tt AVG$.

{\sf {F\lowercase{lash}P}}\xspace facilitates the following class of forecasting models.
Let $M_t$ be the value of metric $M$ ($={\tt SUM}({\bf m}\xspace)$ in our case) at time $t$.
A model is specified by a time-dependent function $f_t$ with order $K$:
\begin{equation}\label{equ:forecastmodel}
M_t = f_t(M_{t-1}, M_{t-2}, \ldots, M_{t-K}).
\end{equation}
The model $f_t$ is {\em fitted} on historical data $M_1, M_2, \ldots, M_{t_0}$ with {\em 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 {\em 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 ({\em i.e.}\xspace, $\hat M_{t_0+1|t_0}$ can be used forecast $M_{t_0+2}$).

For the forecasting models supported by {\sf {F\lowercase{lash}P}}\xspace, we now give brief introduction to both classic ones, {\em e.g.}\xspace, ARMA \cite{book:Hamilton94}, and those based on recurrent neural networks and LSTM \cite{neco:HochreiterS97}.

\stitle{Forecasting using ARMA.}
An ARMA ({\em autoregressive moving average}) \cite{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:
\begin{equation} \label{equ:arma}
{\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}.
\end{equation}

For example, an ARMA model is: $M_t = 0.8 M_{t-1} + 0.2 M_{t-2} + u_t + 0.1 u_{t-1}$. It falls into the form of \eqref{equ:forecastmodel}, 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, {\em e.g.}\xspace, 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 {\em forecast intervals}, {\em i.e.}\xspace, {confidence intervals for forecasts} $\hat M_{t_0+h|t_0}$: with a {\em 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 {\em 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.

\stitle{Forecasting using LSTM.}
LSTM ({\em long short-term memory}) \cite{neco:HochreiterS97} is a network architecture that extends the memory of recurrent neural networks using a {\em 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, {\em e.g.}\xspace, \cite{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$ {\em fully-connected layer} which outputs the forecast of $M_t$.
This model also falls into the general form \eqref{equ:forecastmodel}. 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.

\subsection{Overview of Our Approach}
Our system {\sf {F\lowercase{lash}P}}\xspace works in two online phases to process a forecasting task: {\em preparing training data points by issuing aggregation queries}, and {\em fitting the forecasting model using the training data}.

\squishlist
\item {\em Aggregation query.} In a forecasting task, specified in \eqref{equ:task}, to predict ${\tt SUM}({\bf m}\xspace)$, we have $t_e - t_s + 1$ historical data points:
\begin{align}
M_{t_s} & = \code{SELECT} ~ {\tt SUM}({\bf m}\xspace) ~ \code{FROM} ~ T ~ \code{WHERE} ~ {\sf C}\xspace ~ \code{AND} ~ t = t_s \nonumber
\\
\ldots & \ldots \nonumber
\\
M_{t_e} & = \code{SELECT} ~ {\tt SUM}({\bf m}\xspace) ~ \code{FROM} ~ T ~ \code{WHERE} ~ {\sf C}\xspace ~ \code{AND} ~ t = t_e \label{equ:task:agg}
\end{align}
Each data point is given by an {\em aggregation query} with constraint ${\sf C}\xspace$, which can be any logical expression on the dimension values of ${\bf a}\xspace\supind{1},$ $\ldots,$ ${\bf a}\xspace\supind{d_a}$. We would compute them in the online phase.
\item {\em 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}}$.
\squishend

\stitle{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 $\bigohinline{t_0 \cdot N} + {\rm Train}(t_0, s)$, where $\bigohinline{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, {\em e.g.}\xspace, $(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).

\stitle{Real-time forecasting on approximate aggregations.}
In order to process a forecasting task in an interactive way in {\sf {F\lowercase{lash}P}}\xspace, 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 \eqref{equ:forecastmodel}. 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.

\section{Real-Time Forecasting}
\label{sec:rtforecast}

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, {\em e.g.}\xspace, LSTM-based ones in Section~\ref{sec:exp}.
It is not surprising that there is a tradeoff between sampling rate (or, quality of estimated aggregations) and forecast accuracy.

\eat{
Consider an {\em autoregressive–moving-average} (ARMA) model
\begin{equation}
{\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}
\end{equation}
which is a combination of an {\em auto-regression} model ${\rm AR}(p):$
$M^{\rm AR}_t$ $=$ $\sum_{i=1}^p \alpha_i M_{t-i},$
and a {\em moving average} model ${\rm MA}(q):$
$M^{\rm MA}_t$ $=$ $u_t + \sum_{i=1}^q \beta_i u_{t-i},$
where each $u_t$ at time stamp $t$ is an independent zero-mean random noise. It is a stochastic process that models the evolving of the metric $M$ over time. After fixing $p$ and $q$, the model is fitted by, {\em e.g.}\xspace, least squares regression, to find the values of $\alpha_1 \ldots \alpha_p$ and $\beta_1 \ldots \beta_q$ which minimize the error term.
}

\stitle{Required properties of sampling and estimates.}
In {\sf {F\lowercase{lash}P}}\xspace, 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:

\squishlist
\item ({\em Unbiasedness}) $\hat M_t = M_t + \epsilon\xspace_t$ with $\epinline{\epsilon\xspace_t} = 0$;
\item ({\em Independence}) $\epsilon\xspace_t$'s for different time stamps $t$ are independent.
\squishend
{\em {\rm GSW}\xspace sampling} introduced in Section~\ref{sec:sampler} offers unbiased estimates, with bounded variance of $\epsilon\xspace_t$. {\rm GSW}\xspace samplers are run independently on the data for each $t$---that is how we get independence.

\stitle{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 \eqref{equ:arma} as
\begin{align}
\hat M_t & = \sum_{i=1}^p \alpha_i \hat M_{t-i} + (u_t + \epsilon\xspace_t) + \label{equ:arma:noise}
\\
&
+ \begin{cases}
\sum_{i=1}^p (\beta_i u_{t-i} - \alpha_i \epsilon\xspace_{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\xspace_{t-i}) - \sum_{i = q}^p \alpha_i \epsilon\xspace_{t-i} & (1 < q \leq p)
\end{cases} \nonumber
\end{align}
\eat{
\begin{align}
& \hbox{i) $1 < p \leq q$:} ~~~~ \hat M_t = \sum_{i=1}^p \alpha_i \hat M_{t-i} + \label{equ:arma:1}
\\
& ~~~~~~ + (u_t + \epsilon\xspace_t) + \sum_{i=1}^p (\beta_i u_{t-i} - \alpha_i \epsilon\xspace_{t-i}) + \sum_{i = p}^q \beta_i u_{t-i}. \label{equ:arma:1:noise}
\end{align}
\begin{align}
& \hbox{ii) $1 < q \leq p$:} ~~~~ \hat M_t = \sum_{i=1}^p \alpha_i \hat M_{t-i} + \label{equ:arma:2}
\\
& ~~~~~~ + (u_t + \epsilon\xspace_t) + \sum_{i=1}^q (\beta_i u_{t-i} - \alpha_i \epsilon\xspace_{t-i}) - \sum_{i = q}^p \alpha_i \epsilon\xspace_{t-i}. \label{equ:arma:2:noise}
\end{align}
}
In both cases, it is a combination of an autoregressive model \eqref{equ:arma:noise} of order $p$, and a moving average model of order $\max\{p,q\}$.
The model differs from the ARMA model in \eqref{equ:arma} only on the additional zero-mean error terms $\epsilon\xspace_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, {\em e.g.}\xspace, maximum likelihood estimator, as normal ARMA models. $\epsilon\xspace_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\xspace_t$'s.

In comparison to the original ARMA model in \eqref{equ:arma}, the additional noise term $\epsilon\xspace_{t}$ will indeed incur wider forecast intervals.
However, with a proper sampling-estimation scheme,
if $\epsilon\xspace_t$'s variance is negligible in comparison to $u_t$'s, $\epsilon\xspace_{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:
\begin{proposition}\label{prop:arma11}
{Suppose} we have a time series $\{M_t\}$ satisfies ${\rm ARMA}(1,1)$. Let $\hat M_t = M_t + \epsilon\xspace_t$ be an estimation of $M_t$ satisfying unbiasedness and independence. Then, $\vrinline{\hat M_t} = a \cdot \sigma_u^2 + \sigma_\epsilon\xspace^2$, where $\sigma_u^2 = \vr{u_t}$, $\sigma_\epsilon\xspace^2 = \vr{\epsilon\xspace_t}$, 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)$.
\end{proposition}
We can use normal random variables to approximate $u_t$ and $\epsilon\xspace_t$, and estimate forecast intervals for a given confidence level.

\eat{
In particular, when we use the optimal Bernoulli sampling to estimate $M_t$, we have $\vrinline{\epsilon\xspace_t} = \Delta M_t$ and
\begin{corollary}
If we use optimal Bernoulli sampling to estimate $M_t = \sum_{i=1}^n m_i$ as $\hat M_t$, with sampling parameter $\Delta$, under the condition that $m_i = \bigthetainline{M_t/n}$, the noise term $\epsilon\xspace_t = \hat M_t - M_t$ converges in distribution to a normal random variable:
\[
\frac{1}{\sqrt{\Delta M_t}}\epsilon\xspace_t \xrightarrow{d} N(0, 1).
\]
\end{corollary}
The proof is directly from our analysis of the optimal Bernoulli sampling in Section~\ref{sec:sampler:opt} and Lyapunov Central Limit Theorem.
}

\stitle{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~\ref{fig:lstm:noisy} 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\xspace_{t-i}$ available; these estimates are fed into the LSTM unit as inputs. LSTM then generates an {\em output vector} ${\bf h}_t$ and update its memory {\em 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\xspace_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\xspace_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\xspace_t$'s are small enough). It is difficult to derive any formal analytical result here, but we will evaluate it experimentally in Section~\ref{sec:exp}.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{fig-paper/lstm_noisy}
\vspace{-0.4cm}
\caption{An LSTM-based forecasting model with noisy inputs}
\label{fig:lstm:noisy}
\vspace{-0.6cm}
\end{figure}

\eat{
\subsection{Time-Accuracy Tradeoff: an Example}
\label{sec:rtforecast:progressive}
It is clear that forecasting accuracy largely depends on the variance of the error term $\epsilon\xspace_t$ in the estimates from the above discussion. Meanwhile, processing of aggregation queries is the performance bottleneck in {\sf {F\lowercase{lash}P}}\xspace. Therefore, by varying sample sizes, we can tradeoff forecasting accuracy for faster response time.

In {\sf {F\lowercase{lash}P}}\xspace, samples with different sizes are maintained.
For scenarios with different levels of requirements on the response time, we can choose to use different samples to derive estimates $\hat M_t$'s as training data. We give a numeric example in Figure~\ref{fig:progressive:example}. We use a time series of relation with around 15 million rows per day (refer to Section~\ref{sec:exp} for more information about this dataset), and draw samples with different sizes, from $0.02\%$, $0.1\%$, to $1\%$ of the total data size, using the GSW sampler (with measure ${\sf Impression}$ as sampling weights) which is to be introduced in Section~\ref{sec:sampler}.
A forecasting task with selectivity around $0.5\%$ is to be processed. We use the last 150 days of data to fit the model and forecast the next 7 days (only 90 days' data is plotted as otherwise the 7-day window is too small to be read).
It can be seen that with $0.02\%$-$1\%$ sample size, we achieve interactive response time (30ms-126ms). While the estimations to aggregation queries (red lines) become more and more accurate, the forecast (green lines) gets closer and closer to the true values, with narrower and narrower forecast intervals.

\begin{figure}[t]
\centering
\subfigure[0.02\% sample (30ms)]{
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1\linewidth]{fig-new/vldb-fig9-fix-sample-272704924-5000-pv-regressive-example.eps}
\end{minipage}%
}%
\subfigure[1\% sample (126ms)]{
\begin{minipage}[t]{0.47\linewidth}
\centering
\includegraphics[width=1\linewidth]{fig-new/vldb-fig9-fix-sample-272704924-100-pv-regressive-example.eps}
\end{minipage}
}%
\centering
\vspace{-0.5cm}
\caption{{A numeric example on time-accuracy tradeoff}}
\vspace{-0.6cm}
\label{fig:progressive:example}
\end{figure}
}
\section{Priority v.s. Threshold}
Suppose we have $n$ rows with positive measures $m_1, \ldots, m_n$. Each row generates a random variable $\alpha_i \in (0,1)$. The priority \cite{tit:DuffieldLT05,pods:AlonDLT05,jacm:DuffieldLT07} of each row is calculated as $p_i = m_i/\alpha_i$. For a threshold $\tau > 0$, define a random variable $x_i$ for each row:
\begin{equation}\label{equ:threshold:variable}
x_i =
	\begin{cases}
	0 & \quad \text{if } \frac{m_i}{\alpha_i} < \tau \\
	\tau & \quad \text{if } m_i < \tau \leq \frac{m_i}{\alpha_i} \\
	m_i & \quad \text{if } \tau \leq m_i \\
	\end{cases}.
\end{equation}
In the priority sampling scheme \cite{pods:AlonDLT05,jacm:DuffieldLT07}, $\tau$ is picked to be the $k$th highest priority among all the rows (or a subset of rows).

Let's analyze the distribution of $x_i$. There are two cases:
\squishlist
\item Case i) $m_i < \tau$:
\begin{equation}\label{equ:threshold:variable:probability:1}
\begin{cases}
\pr{x_i = \tau} = \frac{m_i}{\tau}
\\
\pr{x_i = 0} = 1 - \frac{m_i}{\tau}
\end{cases}.
\end{equation}
\item Case ii) $m_i \geq \tau$:
\begin{equation}\label{equ:threshold:variable:probability:2}
\pr{x_i = m_i} = 1 ~\text{ or }~ \pr{x_i = \tau \cdot \frac{m_i}{\tau}} = 1.
\end{equation}
\squishend

We find that \eqref{equ:threshold:variable:probability:1}-\eqref{equ:threshold:variable:probability:2} are also equivalent to the measure-biased sampling \cite{sigmod:DingHCC016}. Intuitively, if $\tau$ is large enough, each row is {\em selected} ($x_i = \tau$) with probability proportional to $m_i$, as in the first line of \eqref{equ:threshold:variable:probability:1}. In particular, we can construct a measure-biased sample $\cal S$ from a priority sample $\{x_i\}$ as follows: if $x_i = \tau$, put one copy of row $i$ into $\cal S$; if $x_i = m_i$, put $\frac{m_i}{\tau}$ copies of row $i$ into $\cal S$.

\stitle{Estimating (subset) sum.} Define the estimator ${\hat M} = \sum_{i=1}^n x_i$. Indeed, it is unbiased
\begin{align*}
\ep{{\hat M}} & = \sum_{i:~m_i < \tau} \ep{x_i} + \sum_{i:~m_i \geq \tau} \ep{x_i}
\\
& = \sum_{i:~m_i < \tau} \tau \cdot \frac{m_i}{\tau} + \sum_{i:~m_i \geq \tau} m_i = \sum_{i=1}^n m_i = M.
\end{align*}

The variance of ${\hat M}$ depends on the value of $\tau$. And indeed, the smaller $\tau$ is, the larger the sample size is. We want to derive the relationship between the variance and the (expected) sample size.

For each $x_i$, in case i) $m_i < \tau$, we have $\vr{x_i} = \ep{x_i^2} - \ep{x_i}^2 = \tau^2 \cdot \frac{m_i}{\tau} - m_i^2 = (\tau - m_i)m_i$; and in case ii) $m_i \geq \tau$, we have $\vr{x_i} = 0$. Since $x_i$'s are independent, we have
\[
\vr{\hat{M}} = \sum_{i=1}^n \vr{x_i} = \sum_{i:~m_i<\tau} (\tau - m_i) m_i.
\]

Regarding the sample size, for a fixed $\tau$, let's define ${\cal S}_{\tau+} = \{i \mid m_i \geq \tau\}$ and ${\cal S}_{\tau-} = \{i \mid \alpha_i \tau \leq m_i < \tau \}$. The sample obtained from \eqref{equ:threshold:variable} is essentially ${\cal S}_\tau = {\cal S}_{\tau+} \cup {\cal S}_{\tau-}$. The set ${\cal S}_{\tau+}$ has a fixed size $|{\cal S}_{\tau+}| = \sum_{i:~m_i \geq \tau} 1$. And the expected size of ${\cal S}_{\tau-}$ is $\ep{|{\cal S}_{\tau-}|} = \sum_{i:~m_i < \tau} \frac{m_i}{\tau}$. Thus, the expected sample size is
\[
\ep{|{\cal S}_\tau|} = \sum_{i=1}^n \frac{\min\{m_i, \tau\}}{\tau} \leq \frac{M}{\tau}.
\]

We can also analyze the relative standard deviation of ${\hat M}$ as in \cite{pods:AlonDLT05}. Simply from its definition, we have
\begin{align*}
\ep{\left(\frac{{\hat M}-M}{M}\right)^2} & = \frac{\sum_{i:~m_i<\tau} (\tau - m_i)m_i}{M^2} \leq \frac{\sum_{i:~m_i<\tau} \tau m_i}{M^2}
\\
& = \frac{\tau^2 \ep{|{\cal S}_{\tau-}|}}{M^2} \leq \frac{\ep{|{\cal S}_{\tau-}|}}{\ep{|{\cal S}_{\tau}|}^2} \leq \frac{1}{\ep{|{\cal S}_{\tau}|}},
\end{align*}
which is consistent with the main result (\cthm1) of \cite{pods:AlonDLT05}.

\section{Related Work}
\label{sec:related}

{\noindent\bf 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 \cite{sigmod:ChaudhuriDK17} for a comprehensive review. There are two major lines of AQP techniques. i) Online aggregation \cite{sigmod:Hellerstein97,nsdi:CondieCAHES10,pvldb:PansareBJC11} and online sampling-based AQP \cite{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 \cite{sigmod:Acharya00,VLDB:ganti00,sigmod:Chaudhuri01,TODS:Chaudhuri07,cidr:OlstonBEJR09,cidr:SidirourgosKB11,Eurosys:Agarwal13}, and some are workload-independent \cite{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~\ref{sec:sampler}.
There are other samplers such as universe (hashed) sampling and stratified sampling introduced in AQP systems \cite{Eurosys:Agarwal13,sigmod:KandulaSVOGCD16,sigmod:ParkMSW18}. These samplers were proposed to handle orthogonal aspects such as missing groups in $\tt Group By$ 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.

\stitle{Precomputing aggregations.}
Another choice is to precompute aggregations or summaries using techniques such as view materialization \cite{VLDB:Agrawal00}, datacube \cite{Springer:Gray1997data}, histograms \cite{PODS:Gilbert01}, and wavelets \cite{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.

\stitle{Aggregation-forecasting analytics.}
\cite{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 \cite{jof:TobiasZ00,book:WestH97} in the forecasting literature. These techniques share some similarity with the Bayesian model in \cite{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 \cite{dasfaa:GongXHCLH15}, but they are less relevant here.
\section{Generalized Weighted Sampler}
\label{sec:sampler}

We now focus on how to draw samples for estimating results of aggregation queries. An aggregation query in \eqref{equ:task:agg} 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}\xspace$ specified in the forecasting task. W.l.o.g., suppose the subset of rows satisfying ${\sf C}\xspace$ is $[n] = \{1,\ldots,n\}$, we want to estimate the metric $M = \sum_{i=1}^n m_i$, for a measure ${\bf m}\xspace = [{m}\xspace_i]$. The (offline) sampling algorithm should be independent on ${\sf C}\xspace$, but could depend on ${\bf m}\xspace$.

\stitle{Existing samplers.}
There are two categories of sampling schemes for estimating subset sums: {\em uniform sampling} and {\em 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 {\em range} of) measure values \cite{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, \cite{sigmod:chaudhuri99} and \cite{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 \cite{tit:DuffieldLT05} and priority sampling \cite{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 \cite{stoc:Szegedy06} in terms of the sampling efficiency with relative standard deviation $\sqrt{\vrinline{\rm estimation}}/M = \sqrt{1/({\rm sample\_size - 1}})$.

\stitle{Requirements for the sampler in {\sf {F\lowercase{lash}P}}\xspace.}
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 {\sf {F\lowercase{lash}P}}\xspace, 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 ({\em e.g.}\xspace, $10$), the storage cost of all the samples ({\em e.g.}\xspace, even with sampling rate $1\%$) is prohibitive in memory.
Therefore, the requirement here is: {\em whether we can generate a compact sample which takes care of multiple measures and still have accuracy guarantees}.

To this end, we propose {\em {\rm GSW}\xspace} ({\em \underline{G}eneralized \underline{S}moothed \underline{W}eighted}) sampling. We introduce its sampling distribution and analyze its sampling efficiency in Section~\ref{sec:sampler:gb}. 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~\ref{sec:sampler:compact}.

\subsection{Generalized Smoothed Sampling}
\label{sec:sampler:gb}
The {\rm GSW}\xspace 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}$, {\em 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 \cite{jasa:HorvitzT52}, we define the {\em 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
\begin{equation}\label{equ:gbin:variable:probability}
\begin{cases}
\prinline{i \in S_\Delta \wedge \hat m_i = \frac{m_i(\Delta + w_i)}{w_i}} = \frac{w_i}{\Delta + w_i}
\\
\prinline{i \notin S_\Delta \wedge \hat m_i = 0} = 1 - \frac{w_i}{\Delta + w_i} = \frac{\Delta}{\Delta + w_i}
\end{cases}.
\end{equation}
We can estimate $M$ as $\hat M = \sum_{i \in S_\Delta} \hat m_i$, which is unbiased since
\[
\ep{{\hat M}} = \sum_{i=1}^n \ep{\hat m_i} = \sum_{i=1}^n \frac{m_i(\Delta + w_i)}{w_i} \cdot \frac{w_i}{\Delta + w_i} = M.
\]

\stitle{Simple and efficient implementations.}
A {\rm GSW}\xspace sampler can be easily implemented in a distributed manner: each row $i$ generates uniformly random number $p_i$ from $[0,1]$, independently; according to \eqref{equ:gbin:variable:probability}, if $p_i \leq \frac{w_i}{\Delta + w_i}$, the row $i$ is put into the sample $S_\Delta$.

A {\rm GSW}\xspace 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'$ and obtain a {\rm GSW}\xspace sample $S_{\Delta'}$ 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'$
from $S_\Delta$, and insert those with $\Delta' \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'}$ {\em without} touching any row in $[n] - S_\Delta$.

\subsubsection{Accuracy Guarantee on Aggregations}
\label{sec:sampler:gb:acc}
We now analyze estimation errors of the class of {\rm GSW}\xspace-based estimators $\hat M$, instantiated by $(\Delta,{\bf w})$.
We consider {\em relative standard deviation} ($\rm RSTD$) and {\em relative error} ($\rm RE$). As $\hat M$ is unbiased,
\begin{align*}
{\rm RSTD}({\hat M}) \! \triangleq \!\! \sqrt{\ep{\left(\frac{{\hat M}-M}{M}\right)^2}}
\geq \ep{\frac{|{\hat M} - M|}{M}} \!\! \triangleq \! {\rm RE}({\hat M}).
\end{align*}

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''.
\begin{definition}\label{def:consistency}
{\bf ($(\underline{\theta}\xspace,\overline{\theta}\xspace)$-consistency)}
Sampling weights ${\bf w}$ $=$ $[w_i]$ are $(\underline{\theta}\xspace,\overline{\theta}\xspace)$-consistent with measure ${\bf m}\xspace = [{m}\xspace_i]$, iff $\underline{\theta}\xspace = \min_i (m_i/w_i)$ and $\overline{\theta}\xspace = \max_i (m_i/w_i)$. $\theta \triangleq \overline{\theta}\xspace/\underline{\theta}\xspace$ is called the consistency scale.
\end{definition}

The above notation about ``consistency'' says, for any row $i$, $\underline{\theta}\xspace \leq m_i/w_i \leq \overline{\theta}\xspace$. It allows ${\bf m}\xspace$ and ${\bf w}$ to differ in scale ({\em e.g.}\xspace, both $\underline{\theta}\xspace$ and $\overline{\theta}\xspace$ could be large) but measures their similarity in patterns and trends. For example, suppose ${\bf m}\xspace$ $=$ $[100,$ $100,$ $200,$ $400]$ and ${\bf w}$ $=$ $[10,$ $10,$ $20,$ $50]$. We have $\underline{\theta}\xspace = 400/50 = 8$, $\overline{\theta}\xspace = 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}$.

\begin{theorem}\label{thm:agg} {\bf (Sampling efficiency of {\rm GSW}\xspace)}
Suppose the sampling weights $\bf w$ used in {\rm GSW}\xspace sampling are $(\underline{\theta}\xspace,\overline{\theta}\xspace)$-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}\xspace/\underline{\theta}\xspace$)
\[
{\rm RE}({\hat M}) \leq {\rm RSTD}({\hat M}) \leq \sqrt{\frac{\overline{\theta}\xspace/\underline{\theta}\xspace}{\ep{|{\mathcal S}_{\Delta}|}}} = \sqrt{\frac{\theta}{\ep{|{\mathcal S}_{\Delta}|}}}.
\]
\end{theorem}

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

\subsubsection{A Special Case: Optimal GSW Sampler}
\label{sec:sampler:opt}
We can choose $\bf w$ to minimize the variance of estimation, while the expected sample size is bounded.
\ifdefined
	Please refer to Appendix~\ref{app:optweight} for a precise formulation and the optimal solution.
\else
	Please refer to \cite{url:tr} for a precise formulation and the optimal solution.
\fi
From Theorem~\ref{thm:agg}, a nearly-optimal solution is ${\bf w} = {\bf m}\xspace$, as in this case, $\bf w$ is $(1,1)$-consistent with $\bf m$ ($\theta = 1$). We call it {\em optimal {\rm GSW}\xspace sampler}. Directly from Theorem~\ref{thm:agg}, we have
\begin{corollary}\label{cor:agg:opt}
{\bf (Optimal {\rm GSW}\xspace 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}{\ep{|{\mathcal S}_{\Delta}|}}}.
\]
\end{corollary}

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

\eat{
Suppose we have $n$ rows with positive measures $m_1, \ldots, m_n$. For each row, let's consider the following distribution:
\begin{equation}\label{equ:bin:variable:probability}
\begin{cases}
\pr{x_i = \Delta + m_i} = \frac{m_i}{\Delta + m_i}
\\
\pr{x_i = 0} = 1 - \frac{m_i}{\Delta + m_i}
\end{cases}.
\end{equation}

\stitle{Estimating (subset) sum.} Define the estimator ${\hat M} = \sum_{i=1}^n x_i$. Indeed, it is unbiased
\begin{align*}
\ep{{\hat M}} & = \sum_{i=1}^n \ep{x_i} = \sum_{i=1}^n \frac{(\Delta + m_i) m_i}{\Delta + m_i} = \sum_{i=1}^n m_i = M.
\end{align*}

The variance of ${\hat M}$ depends on the value of $\Delta$. And indeed, the smaller $\Delta$ is, the larger the sample size is. We want to derive the relationship between the variance and the (expected) sample size.

For each $x_i$, we have
\[
\vr{x_i} = (\Delta + m_i)^2 \cdot \frac{m_i}{\Delta + m_i} \cdot \left(1 - \frac{m_i}{\Delta + m_i}\right) = \Delta m_i.
\]

Since $x_i$'s are independent, we have
\[
\vr{\hat{M}} = \sum_{i=1}^n \vr{x_i} = \Delta M.
\]

Regarding the sample size, for a fixed $\Delta$, let's define ${\cal S}_{\Delta} = \{i \mid x_i = \Delta + m\}$, which is essentially the sample obtained from \eqref{equ:bin:variable:probability}. The expected sample size is
\[
\ep{|{\cal S}_\Delta|} = \sum_{i=1}^n  \frac{m_i}{\Delta + m_i} \leq \frac{M}{\Delta}.
\]

We can also analyze the relative standard deviation of ${\hat M}$ as in \cite{pods:AlonDLT05}. Simply from its definition, we have
\begin{align*}
\ep{\left(\frac{{\hat M}-M}{M}\right)^2} & = \frac{\Delta M}{M^2} = \frac{\Delta}{M} \leq \frac{1}{\ep{|{\cal S}_{\Delta}|}}.
\end{align*}
}

\eat{
\subsubsection{Comparison with Other Samplers}
\label{sec:sampler:compare}
Uniform sampling is a simple but effective sampling scheme for estimating aggregations \cite{sigmod:Hellerstein97}. Each row in the subset $[n]$ is drawn into a sample $S$ with equal probability $p$, an unbiased estimation for $M$ is $\frac{1}{p}\sum_{i \in S} \frac{m_i}{p}$. The main advantage is that its sampling distribution is independent on the measure's distribution; and thus no matter how many measures there are in the relation, one sample is sufficient.
However, the main deficiency is that the estimation has large error for heavy tailed distributions, which is common in practice. The analytical error bound is proportional to the difference between the maximum and minimum measure values \cite{sigmod:Hellerstein97}.

Weighted sampling schemes give estimators for ${\tt SUM}$ aggregations with (relative) error bounds dependent on only the sample size. In the simplest form, for some fixed constant $B$ (deciding the sample size), each row $i$ is drawn into the sample with probability $m_i/B$; if $m_i > B$, multiple copies of $i$ are drawn into the sample. \cite{sigmod:chaudhuri99} and \cite{sigmod:DingHCC016} give efficient implementations of this scheme. One disadvantage of this sampling scheme is that, for heavy tailed distribution, the sample is dominated by copies of rows at the tail.

Threshold sampling \cite{tit:DuffieldLT05} is a kind of Poisson sampling that resolves the above disadvantage by including each row in the sample at most once. Priority sampling \cite{pods:AlonDLT05,jacm:DuffieldLT07} is its variant in the streaming and online-aggregation setting. Priority sampling has been shown to be optimal \cite{stoc:Szegedy06} in terms of the sampling efficiency with ${\rm RSTD} = \sqrt{1/{\rm sample\_size - 1}}$. In comparison, our optimal Bernoulli sampler introduced in Section~\ref{sec:sampler:opt} has nearly-optimal sampling efficiency with ${\rm RSTD} = \sqrt{1/\epinline{\rm sample\_size}}$

There are also other samplers such as universe (hashed) sampling and stratified sampling introduced in the literatures of approximate query processing \cite{Eurosys:Agarwal13,sigmod:KandulaSVOGCD16,sigmod:ParkMSW18}. These samplers were proposed to handle orthogonal aspects such as missing groups in $\tt Group By$ 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.
}

\subsection{Compact Sample for Multiple Measures}
\label{sec:sampler:compact}
The size of {\rm GSW}\xspace sample can be controlled by the parameter $\Delta$. When there is only one measure ${\bf m}\xspace$, we draw an optimal {\rm GSW}\xspace sample (setting ${\bf w} = {\bf m}\xspace$). A more common scenario is that we have multiple measures ({\em e.g.}\xspace, in Figure~\ref{fig:motivation}) in one relation. We can draw one optimal {\rm GSW}\xspace 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}\xspace\supind{1},$ $\ldots,$ ${\bf m}\xspace\supind{k}$. For each measure $j$, we want to estimate $M\supind{j} = \sum_{i=1}^n {m}\xspace\supind{j}_i$ for rows in a set $[n]$ (satisfying the constraint ${\sf C}\xspace$ in a forecasting task).
A {\rm GSW}\xspace sample $S_\Delta$ is drawn using weights ${\bf w} = [w_i]$.
The {\em calibrated measure} on each sample row $i$ for each measure $j$ is $\hat m\supind{j}_i = m\supind{j}_i(\Delta + w_i)/w_i$. From Theorem~\ref{thm:agg}, $\hat M\supind{j} = \sum_{i \in S_\Delta} \hat m\supind{j}_i$ is an unbiased estimation of $M\supind{j}$.

Intuitively, if the chosen sampling weight vector $\bf w$ centers around ${\bf m}\xspace\supind{1},$ $\ldots,$ ${\bf m}\xspace\supind{k}$ and is not too far away from any of the $k$, from Theorem~\ref{thm:agg}, the error can be better bounded. We can find such centers by taking the average of measures. For example, for ${\bf m}\xspace\supind{1} = [100, 100, 200, 400]$ and ${\bf m}\xspace\supind{2} = [1, 1, 2, 1]$, the {\em geometric mean} is ${\bf w}^\times = [\sqrt{100 \cdot 1} = 10, 10, 20, 20]$, and the {\em 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.

\stitle{Geometric compressed {\rm GSW}\xspace sampling.} We can use the geometric mean of the $k$ measures as the sampling weights:
\begin{equation}
w^\times_i = \left(\prod_{j=1}^k {m}\xspace\supind{j}_{i}\right)^{1/k}.
\end{equation}

Among the $k$ measures to be grouped, define the {\em trend deviation} between any two measures ${\bf m}\xspace\supind{p}$, ${\bf m}\xspace\supind{q}$ (for $p,q \in [k]$):
\begin{equation}
\overline\rho_{p,q} \triangleq \max_{i\in[n]} \frac{{m}\xspace\supind{p}_{i}}{{m}\xspace\supind{q}_{i}},
~~~
\underline\rho_{p,q} \triangleq \min_{i\in[n]} \frac{{m}\xspace\supind{p}_{i}}{{m}\xspace\supind{q}_{i}},
~~~ \hbox{and} ~ \rho_{p,q} \triangleq \frac{\overline\rho_{p,q}}{\underline\rho_{p,q}},
\end{equation}
which measures how similar the trends (instead of their absolute values) of the two measures are\footnote{Or, equivalently, ${\bf m}\xspace\supind{q}$ is $(\underline\rho_{p,q}, \overline\rho_{p,q})$-consistent with ${\bf m}\xspace\supind{p}$.}. The smaller $\rho_{p,q}$ is, the more similar ${\bf m}\xspace\supind{p}$ and ${\bf m}\xspace\supind{q}$ are. For example, if ${\bf m}\xspace\supind{p} = c \cdot {\bf m}\xspace\supind{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 {\em maximum deviation} among the $k$ measures. From Theorem~\ref{thm:agg},

\begin{corollary}\label{cor:agg:geo}
If we use ${\bf w}^\times = [w^\times_i]$ as the sampling weights for a relation with $k$ measures, with a {\rm GSW}\xspace sample $S_\Delta$, we can estimate $M\supind{p}$ as $\hat M\supind{p}$ for each measure $p$ with error
\[
{\rm RE}({\hat M}\supind{p}) \leq {\rm RSTD}({\hat M}\supind{p}) \leq \sqrt{\frac{\prod_{j: ~ j \neq p}\rho_{p,j}^{1/k}}{\ep{|{\mathcal S}_{\Delta}|}}} \leq \sqrt{\frac{\rho^{\frac{k-1}{k}}}{\ep{|{\mathcal S}_{\Delta}|}}}.
\]
\end{corollary}

\stitle{Arithmetic compressed {\rm GSW}\xspace sampling.}
Another choice is to use the arithmetic mean as the sampling weights:
\begin{equation}
w^+_i = \frac{1}{k} \sum_{j=1}^k {m}\xspace\supind{j}_{i}.
\end{equation}

Define the {\em 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:
\begin{equation}
\delta \triangleq \max_{i\in[n]}\left(\frac{\max_{j\in[k]} {m}\xspace\supind{j}_{i}}{\min_{j\in[k]} {m}\xspace\supind{j}_{i}}\right).
\end{equation}

From the definitions, for any measure $p$, we have $1/\delta$ $\leq$ ${m}\xspace\supind{p}_{i}/w^+_i$ $\leq$ $\delta$. Thus, directly from Theorem~\ref{thm:agg}, we have the following bound.
\begin{corollary}\label{cor:agg:ari}
If we use ${\bf w}^+ = [w^+_i]$ as the sampling weights for a relation with $k$ measures, with a {\rm GSW}\xspace sample $S_\Delta$, we can estimate $M\supind{p}$ as $\hat M\supind{p}$ for each measure $p$ with error
\[
{\rm RE}({\hat M}\supind{p}) \leq {\rm RSTD}({\hat M}\supind{p}) \leq \sqrt{\frac{\delta^2}{\ep{|{\mathcal S}_{\Delta}|}}}.
\]
\end{corollary}

\stitle{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 ({\em e.g.}\xspace, $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 {\rm GSW}\xspace sample gives good estimations for the measures. To this end, we establish a relationship between $(\underline{\theta}\xspace,\overline{\theta}\xspace)$-consistency and $L_1$ distance as follows.
\begin{proposition}\label{prop:consistencyl2}
{\bf (consistency and $L_1$ distance)}
Suppose sampling weights ${\bf w}$ $=$ $[w_i]$ are $(\underline{\theta}\xspace,\overline{\theta}\xspace)$-consistent with measure ${\bf m}\xspace = [{m}\xspace_i]$. We normalize $\bf w$ as ${\bf w}' = [w'_i = w_i/\sum_j w_j]_i$, and similarly, ${\bf m}\xspace$ as ${\bf m}\xspace'$. Let $\theta = \overline{\theta}\xspace/\underline{\theta}\xspace$. We have $\|{\bf m}\xspace' - {\bf w}'\|_1 \leq (\theta - 1)$.
\end{proposition}

$(\underline{\theta}\xspace,\overline{\theta}\xspace)$-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~\ref{prop:consistencyl2} to quantify the {\em 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 {{\sc KCenter} problem} \cite{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 \cite{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~\ref{prop:consistencyl2} and the triangle inequality in $L_1$, at least we know that we'd better not group two measures that are far way ({\em e.g.}\xspace, $> 2(\theta-1)$) together, as there is no sampling weight vector that is consistent with both of them at the same time.
\ifdefined
	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~\ref{sec:exp}, the experimental results are plotted in Figure~\ref{fig:exp:correlation}. Here, the four measures, {\sf Impression}, {\sf Click}, {\sf Favorite}, and {\sf Cart}, are partitioned into two equal-size groups (there are a total of three ways of partitioning). For example, the first way in Figure~\ref{fig:exp:correlation} is group1 = \{{\sf Impression}, {\sf Click}\} and group2 = \{{\sf Favorite}, {\sf Cart}\}. For each group, we pick the sampling weight vector as the arithmetic mean of measure vectors in this group. Figure~\ref{fig:exp:correlation:l1} plots the $L_1$ distance between each measure vector and the sampling weight vector, and Figure~\ref{fig:exp:correlation:aerr} plots the aggregation error using this sampling weight vector in {\rm GSW}\xspace. 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).
\fi
\ifdefined
	In Section~\ref{sec:exp}, we focus on experimental evaluation after the grouping is fixed. We will leave the development and analysis of better grouping strategies as future work.
\else
	A preliminary evaluation of this grouping strategy can be found in \cite{url:tr}.

\fi

\ifdefined
\begin{figure}[t]
\begin{minipage}[t]{0.46\linewidth}
\subfigure[Three ways to group the four measures: {\sf Imp} = {\sf Impression}, {\sf clk} = {\sf Click}, {\sf fav} = {\sf Favorite}, {\sf cart} = {\sf Cart}]{
\centering
\includegraphics[width=1\linewidth]{fig-paper-lstm/vldb-fig4-fix-sample-100-group.eps}
\label{fig:exp:correlation:aerr}
}
\end{minipage}
\hspace{0.26cm}
\begin{minipage}[t]{0.46\linewidth}
\subfigure[$L_1$ distance between measure vectors and the sampling weight vector (the arithmetic mean of measure vectors in a group)]{
\centering
\includegraphics[width=1\linewidth]{fig-paper-lstm/vldb-fig4-fix-sample-theta.eps}
\label{fig:exp:correlation:l1}
}
\end{minipage}
\vspace{-0.5cm}
\caption{Aggregation error and $L_1$ distance}
\vspace{-0.3cm}
\label{fig:exp:correlation}
\end{figure}
\fi

\eat{
\subsection{Optimal Bernoulli Sampling}
\label{sec:sampler:opt}
Suppose we have $n$ rows with positive measures $m_1, \ldots, m_n$. For each row, let's consider the following distribution:
\begin{equation}\label{equ:bin:variable:probability}
\begin{cases}
\pr{x_i = \Delta + m_i} = \frac{m_i}{\Delta + m_i}
\\
\pr{x_i = 0} = 1 - \frac{m_i}{\Delta + m_i}
\end{cases}.
\end{equation}

\stitle{Estimating (subset) sum.} Define the estimator ${\hat M} = \sum_{i=1}^n x_i$. Indeed, it is unbiased
\begin{align*}
\ep{{\hat M}} & = \sum_{i=1}^n \ep{x_i} = \sum_{i=1}^n \frac{(\Delta + m_i) m_i}{\Delta + m_i} = \sum_{i=1}^n m_i = M.
\end{align*}

The variance of ${\hat M}$ depends on the value of $\Delta$. And indeed, the smaller $\Delta$ is, the larger the sample size is. We want to derive the relationship between the variance and the (expected) sample size.

For each $x_i$, we have
\[
\vr{x_i} = (\Delta + m_i)^2 \cdot \frac{m_i}{\Delta + m_i} \cdot \left(1 - \frac{m_i}{\Delta + m_i}\right) = \Delta m_i.
\]

Since $x_i$'s are independent, we have
\[
\vr{\hat{M}} = \sum_{i=1}^n \vr{x_i} = \Delta M.
\]

Regarding the sample size, for a fixed $\Delta$, let's define ${\cal S}_{\Delta} = \{i \mid x_i = \Delta + m\}$, which is essentially the sample obtained from \eqref{equ:bin:variable:probability}. The expected sample size is
\[
\ep{|{\cal S}_\Delta|} = \sum_{i=1}^n  \frac{m_i}{\Delta + m_i} \leq \frac{M}{\Delta}.
\]

We can also analyze the relative standard deviation of ${\hat M}$ as in \cite{pods:AlonDLT05}. Simply from its definition, we have
\begin{align*}
\ep{\left(\frac{{\hat M}-M}{M}\right)^2} & = \frac{\Delta M}{M^2} = \frac{\Delta}{M} \leq \frac{1}{\ep{|{\cal S}_{\Delta}|}}.
\end{align*}
}

\eat{
\begin{theorem}[Lindberg-Levy Central Limit Theorem]
Let $\left\{X_{1}, ..., X_{n}\right\}$ be a simple random sample of size $n$, i.e. a sequence of $n$
i.i.d random variables, drawn from a distribution $X$ with $E(X)=\mu$ and $Var(X)=\sigma^2<\infty$, there are:
\begin{displaymath}
 \lim \limits_{n \to \infty}P\left \{\frac{\sum_{k=1}^{n}X_{k}-n\cdot \mu}{\sqrt{n}\cdot \sigma }\leq x\right \}=\phi (x)
\end{displaymath}
\end{theorem}

\begin{displaymath}
\lim \limits_{n \to \infty}P\left \{\frac{\frac{\sum_{k=1}^{n}\tilde{M}_{k}}{n}-M}{M}\leq x\cdot \sqrt{\frac{\Delta}{M\cdot n}}\right \}=\phi (x) \\
\end{displaymath}
\begin{equation}
\epsilon \sim z\cdot \sqrt{\frac{\Delta}{M}}
\end{equation}
}

\eat{
\subsection{Comparison to Known Sampling Schemes}
Suppose we have $n$ rows with positive measures $m_1, \ldots, m_n$. Each row generates a random variable $\alpha_i \in (0,1)$. The priority \cite{tit:DuffieldLT05,pods:AlonDLT05,jacm:DuffieldLT07} of each row is calculated as $p_i = m_i/\alpha_i$. For a threshold $\tau > 0$, define a random variable $x_i$ for each row:
\begin{equation}\label{equ:threshold:variable}
x_i =
	\begin{cases}
	0 & \quad \text{if } \frac{m_i}{\alpha_i} < \tau \\
	\tau & \quad \text{if } m_i < \tau \leq \frac{m_i}{\alpha_i} \\
	m_i & \quad \text{if } \tau \leq m_i \\
	\end{cases}.
\end{equation}
In the priority sampling scheme \cite{pods:AlonDLT05,jacm:DuffieldLT07}, $\tau$ is picked to be the $k$th highest priority among all the rows (or a subset of rows).

Let's analyze the distribution of $x_i$. There are two cases:
\squishlist
\item Case i) $m_i < \tau$:
\begin{equation}\label{equ:threshold:variable:probability:1}
\begin{cases}
\pr{x_i = \tau} = \frac{m_i}{\tau}
\\
\pr{x_i = 0} = 1 - \frac{m_i}{\tau}
\end{cases}.
\end{equation}
\item Case ii) $m_i \geq \tau$:
\begin{equation}\label{equ:threshold:variable:probability:2}
\pr{x_i = m_i} = 1 ~\text{ or }~ \pr{x_i = \tau \cdot \frac{m_i}{\tau}} = 1.
\end{equation}
\squishend

We find that \eqref{equ:threshold:variable:probability:1}-\eqref{equ:threshold:variable:probability:2} are also equivalent to the measure-biased sampling \cite{sigmod:DingHCC016}. Intuitively, if $\tau$ is large enough, each row is {\em selected} ($x_i = \tau$) with probability proportional to $m_i$, as in the first line of \eqref{equ:threshold:variable:probability:1}. In particular, we can construct a measure-biased sample $\cal S$ from a priority sample $\{x_i\}$ as follows: if $x_i = \tau$, put one copy of row $i$ into $\cal S$; if $x_i = m_i$, put $\frac{m_i}{\tau}$ copies of row $i$ into $\cal S$.

\stitle{Estimating (subset) sum.} Define the estimator ${\hat M} = \sum_{i=1}^n x_i$. Indeed, it is unbiased
\begin{align*}
\ep{{\hat M}} & = \sum_{i:~m_i < \tau} \ep{x_i} + \sum_{i:~m_i \geq \tau} \ep{x_i}
\\
& = \sum_{i:~m_i < \tau} \tau \cdot \frac{m_i}{\tau} + \sum_{i:~m_i \geq \tau} m_i = \sum_{i=1}^n m_i = M.
\end{align*}

The variance of ${\hat M}$ depends on the value of $\tau$. And indeed, the smaller $\tau$ is, the larger the sample size is. We want to derive the relationship between the variance and the (expected) sample size.

For each $x_i$, in case i) $m_i < \tau$, we have $\vr{x_i} = \ep{x_i^2} - \ep{x_i}^2 = \tau^2 \cdot \frac{m_i}{\tau} - m_i^2 = (\tau - m_i)m_i$; and in case ii) $m_i \geq \tau$, we have $\vr{x_i} = 0$. Since $x_i$'s are independent, we have
\[
\vr{\hat{M}} = \sum_{i=1}^n \vr{x_i} = \sum_{i:~m_i<\tau} (\tau - m_i) m_i.
\]

Regarding the sample size, for a fixed $\tau$, let's define ${\cal S}_{\tau+} = \{i \mid m_i \geq \tau\}$ and ${\cal S}_{\tau-} = \{i \mid \alpha_i \tau \leq m_i < \tau \}$. The sample obtained from \eqref{equ:threshold:variable} is essentially ${\cal S}_\tau = {\cal S}_{\tau+} \cup {\cal S}_{\tau-}$. The set ${\cal S}_{\tau+}$ has a fixed size $|{\cal S}_{\tau+}| = \sum_{i:~m_i \geq \tau} 1$. And the expected size of ${\cal S}_{\tau-}$ is $\ep{|{\cal S}_{\tau-}|} = \sum_{i:~m_i < \tau} \frac{m_i}{\tau}$. Thus, the expected sample size is
\[
\ep{|{\cal S}_\tau|} = \sum_{i=1}^n \frac{\min\{m_i, \tau\}}{\tau} \leq \frac{M}{\tau}.
\]

We can also analyze the relative standard deviation of ${\hat M}$ as in \cite{pods:AlonDLT05}. Simply from its definition, we have
\begin{align*}
\ep{\left(\frac{{\hat M}-M}{M}\right)^2} & = \frac{\sum_{i:~m_i<\tau} (\tau - m_i)m_i}{M^2} \leq \frac{\sum_{i:~m_i<\tau} \tau m_i}{M^2}
\\
& = \frac{\tau^2 \ep{|{\cal S}_{\tau-}|}}{M^2} \leq \frac{\ep{|{\cal S}_{\tau-}|}}{\ep{|{\cal S}_{\tau}|}^2} \leq \frac{1}{\ep{|{\cal S}_{\tau}|}},
\end{align*}
which is consistent with the main result (\cthm1) of \cite{pods:AlonDLT05}.
}

\eat{
\subsection{Generalized Bernoulli Sampling}
In the optimal Bernoulli sampling, the sampling probability \eqref{equ:bin:variable:probability} is proportional to the measure $m_i$ (when $\Delta$ is large). Now the question is, if we use a sampling weight $w_i$ instead of $m_i$, how good the estimation could be. Consider the following distribution:
\begin{equation}\label{equ:gbin:variable:probability}
\begin{cases}
\pr{x_i = m_i(\Delta + w_i)/w_i} = \frac{w_i}{\Delta + w_i}
\\
\pr{x_i = 0} = 1 - \frac{w_i}{\Delta + w_i}
\end{cases}.
\end{equation}

\stitle{Estimating (subset) sum.} Again, we use ${\hat M} = \sum_{i=1}^n x_i$ to estimate $M$. Indeed, it is unbiased
\begin{align*}
\ep{{\hat M}} & = \sum_{i=1}^n \ep{x_i} = \sum_{i=1}^n m_i = M.
\end{align*}
The variance of $x_i$ is
\[
\vr{x_i} = \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 $x_i$'s are independent, we have
\[
\vr{\hat{M}} = \sum_{i=1}^n \vr{x_i} = \sum_{i=1}^n\frac{\Delta m_i^2}{w_i}.
\]

Regarding the sample size, for a fixed $\Delta$, let's define ${\cal S}_{\Delta} = \{i \mid x_i = m_i(\Delta + w_i)/w_i\}$, which is essentially the sample obtained from \eqref{equ:gbin:variable:probability}. The expected sample size is
\[
\ep{|{\cal S}_\Delta|} = \sum_{i=1}^n  \frac{w_i}{\Delta + w_i} \leq \frac{W}{\Delta}.
\]

Intuitively, the closer $w_i$ and $m_i$ are, the more accurate the estimator is. If $\alpha \leq m_i/w_i \leq \beta$, we can derive
\[
\vr{\hat{M}} \leq \beta \Delta M
\]
and
\[
\ep{|{\cal S}_\Delta|} \leq \frac{W}{\Delta} \leq \frac{M}{\alpha \Delta}.
\]
Therefore,
\begin{align}
\ep{\left(\frac{{\hat M}-M}{M}\right)^2} & \leq \frac{\beta \Delta M}{M^2} = \frac{\beta\Delta}{M} \leq \frac{\beta/\alpha}{\ep{|{\cal S}_{\Delta}|}}. \label{equ:gbin:re}
\end{align}
In particular, if $\alpha = \beta$, the sampling efficiency is the same as the optimal Bernoulli sample.
}

\section{Conclusions}
We introduce a real-time forecasting system {\sf {F\lowercase{lash}P}}\xspace for analyzing time-series relations. There are two core technical contributions, which enable {\sf {F\lowercase{lash}P}}\xspace 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 {\rm GSW}\xspace sampling, a new sampling scheme; we establish a connection from the difference between {\rm GSW}\xspace'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 {\rm GSW}\xspace sample with the means of these measures' values as sampling weights.
\else
\subsubsection*{Acknowledgments}
We thank the anonymous reviewers for their helpful comments that improved the quality of the
paper. Zhewei Wei was supported by NSFC No.61832017, No. 61972401 and No. 61932001, by the Fundamental Research Funds for the Central Universities and the Research Funds of Renmin University of China under Grant 18XNLG21, by Beijing Outstanding Young Scientist Program NO. BJJWZYJH01201910002009, and by Alibaba Group through Alibaba Innovative Research Program.
\fi

redpajama_set_name: RedPajamaArXiv