% most useful Latex command \newcommand{\comment}[1]{} \documentclass[10pt,twocolumn]{article} \usepackage{latex8} \usepackage{times} %\usepackage{fullpage} \usepackage{epsf} %\usepackage{changebar} %\usepackage{lscape} %\input{restore} \setlength{\floatsep}{0pt} \setlength{\textfloatsep}{2pt} \setlength{\columnsep}{14pt} \renewcommand{\dbltopfraction}{.9} \renewcommand{\topfraction}{.9} \renewcommand{\textfraction}{.1} \renewcommand{\floatpagefraction}{.9} %\renewcommand{\dbltopnumber}{1} \renewcommand{\dblfloatpagefraction}{.9} \def\normspacing {\renewcommand{\baselinestretch}{0.9}} \def\tinyspacing {\renewcommand{\baselinestretch}{0.7}} %\def\normspacing {} %\def\tinyspacing {} \setlength{\oddsidemargin}{-.231in} \setlength{\evensidemargin}{-.231in} %\setlength{\oddsidemargin}{-.296in} %\setlength{\evensidemargin}{-.296in} %\setlength{\hoffset}{1.1in} \normspacing \def\thesection{\arabic{section}} \def\thefigure{\arabic{figure}} %\newcounter {subsubsubsection}[subsubsection] %\def\thesubsubsubsection {\thesubsubsection.\arabic{subsubsubsection}} %\def\subsubsubsection {\subsubsection*} \makeatletter\long\def\unmarkedfootnote#1{{\long\def\@makefntext##1{##1}\footnotetext{#1} }} \def\pagefigwidth {\epsfxsize=5in} \def\colfigwidth {\epsfxsize=2.5in} \def\figfontsize {\tiny} \def\leadfigwidth{\epsfxsize=2in} %\def\absfigsize {\epsfxsize=3.0in} %\def\relfigsize {\epsfxsize=3.0in} \newcounter{ecount} \newenvironment{ecompact}{ \begin{list}% {\arabic{ecount}}{\usecounter{ecount} \parsep 0pt plus 1pt \partopsep 0pt plus 1pt \topsep 2pt plus 2pt minus 1pt \itemsep 0pt plus 1pt}}% {\end{list}} \newenvironment{icompact}{ \begin{list}{$\bullet$}{ \parsep 0pt plus 1pt \partopsep 0pt plus 1pt \topsep 2pt plus 2pt minus 1pt \itemsep 0pt plus 1pt}}% {\end{list}} \begin{document} \bibliographystyle{acm} \title{An Evaluation of Linear Models for Host Load Prediction} \author{Peter A. Dinda \hspace{.5in} David R. O'Hallaron \\ Carnegie Mellon University\\ 5000 Forbes Avenue \\ Pittsburgh, PA 15213 \\ \{pdinda,droh\}@cs.cmu.edu\\ } \maketitle \unmarkedfootnote{ Effort sponsored in part by the Advanced Research Projects Agency and Rome Laboratory, Air Force Materiel Command, USAF, under agreement number F30602-96-1-0287, in part by the National Science Foundation under Grant CMS-9318163, and in part by a grant from the Intel Corporation. } \begin{abstract} This paper evaluates linear models for predicting the Digital Unix five-second host load average from 1 to 30 seconds into the future. A detailed statistical study of a large number of long, fine grain load traces from a variety of real machines leads to consideration of the Box-Jenkins models (AR, MA, ARMA, ARIMA), and the ARFIMA models (due to self-similarity.) These models, as well as a simple windowed-mean scheme, are then rigorously evaluated by running a large number of randomized testcases on the load traces and data-mining their results. The main conclusions are that load is consistently predictable to a very useful degree, and that the simpler models such as AR are sufficient for doing this prediction. % % This says everything important but is still too long % %The running time of a task is strongly related to the host load it %experiences, as measured by the Unix one-minute load average. High %quality predictions of that load signal can greatly reduce the %uncertainty in the expected running time of tasks, sometimes by %factors of five or more. We studied the prediction of load for 1 to %30 seconds in the future given 1 Hz samples of it's past behavior. We %began by collecting and analyzing a large number of long load traces. %The analysis suggested that the Box-Jenkins linear time series models %(AR, MA, ARMA, ARIMA) might be appropriate for predicting load but %that such models might require an infeasible number of parameters. We %also found self-similarity, which suggested that expensive ARFIMA %models might be appropriate. We next evaluated the Box-Jenkins and %ARFIMA models, as well as a simple ad hoc scheme. We ran a large %number of randomized testcases on the traces and used the results to %rigorously compare the models' mean squared error levels as well as %their consistency in providing this level of quality. We found that %load is consistently predictable to a very useful degree from past %behavior and that simple linear models are sufficient. Although there %are statistically significant differences in the performance of the %different models, these differences are not large enough to warrant %using more the sophisticated models given their much higher run-time %costs and the possibility of instability. Our recommendation is to use %AR models of order 16 or better. \end{abstract} \Section{Introduction} \label{sec:intro} Consider an application program that wants to schedule a compute-bound soft real-time task in a typical distributed computing environment~\cite{DINDA-CASE-BERT-WPDRTS-99}. By using mechanisms such as CORBA~\cite{CORBA-20-ARCH-SPEC-TECHREPORT} or Java RMI~\cite{JAVA-RMI-SPEC}, the task can be executed on any of the host machines in the network. Given this freedom, the application can choose the host on which the task's deadline is most likely to be met. If the application could predict the exact running time of the task on each of the hosts, choosing a host would be easy. However, such predictions are unlikely to be exact due to the dynamic nature of the distributed system. Each of the hosts is acting independently, its vendor-supplied operating system scheduling tasks initiated by other users, paying no special attention to the task the application is trying to schedule. The computational load on each of the hosts can vary drastically over time. Because of this dynamic nature, the application's predictions of the task's running time on each of the hosts have confidence intervals associated with them. Better predictions lead to smaller confidence intervals, which makes it easier for the application to choose between the hosts, or to decide how likely a particular host is to meet the deadline. The running time of a compute-bound task on a host is strongly related to the computational load on the host. If we could predict the load that a task would encounter on some host, we could predict the running time on that host. Better load predictions lead to better predictions of running time and thus to smaller confidence intervals. For this reason, we concentrate on host load prediction here. \comment{ \begin{figure} \centerline{ \colfigwidth \epsfbox{unix16-simple.epsf} } \caption{Benefits of prediction: server machine with 40 users and long term load average of 10.} \label{fig:extreme} \end{figure} Better load predictions can indeed lead to drastically smaller confidence intervals on real hosts, even heavily loaded ones. For example, Figure~\ref{fig:extreme} plots, for one such host, the length of the confidence interval for the running time of a one second task as a function of how far ahead the load predictions are made. The confidence intervals for a predictive AR($9$) model (described in Section~\ref{sec:ltsmodels}) and for the raw variance of the load signal itself are represented. Notice that for up to 15 seconds into the future, the AR($9$) model provides a confidence interval that is less than half the length of that of the raw load signal. } \begin{figure} \centerline{ \colfigwidth \epsfbox{axpfea-simple.epsf} } \caption{Benefits of prediction: length of 95\% confidence interval for running time of a one second task on an interactive cluster machine with long term load average of 0.17.} \label{fig:lessextreme} \end{figure} Better load predictions can indeed lead to drastically smaller confidence intervals on real hosts, even lightly loaded ones. For example, Figure~\ref{fig:lessextreme} plots, for one such host, the length of the confidence interval for the running time of a one second task as a function of how far ahead the load predictions are made. The confidence intervals for a predictive AR($18$) model (described in Section~\ref{sec:ltsmodels}) and for the raw variance of the load signal itself are represented. Notice that for up to 25 seconds into the future, the AR($18$) model provides a confidence interval of less than 200 ms while the raw load signal provides a confidence interval of over two seconds. The {\em existence} of clear benefits of this kind motivates the following questions: Is host load {\em consistently} predictable or are examples like Figure~\ref{fig:lessextreme} merely flukes? If host load is indeed consistently predictable, what classes of predictive models are appropriate for predicting it? What are the differences between these different classes? Finally, what class can be recommended for use in a real system? This paper describes the results of a large scale, real world study to provide statistically rigorous answers to these questions. We found that host load is, in fact, consistently predictable to a very useful degree from past behavior, and that simple, practical linear time series models are sufficiently powerful load predictors. These results are somewhat surprising because load has complex behavior and exhibits properties such as self-similarity and epochal behavior that suggest that more complex models would be more appropriate. As far as we are aware, this is the first study to identify these properties of host load and then to rigorously evaluate the practical predictive power of linear time series models that both can and cannot capture them. Furthermore, our evaluation approach, while much more computationally intensive than those of traditional time series analysis, is unbiased and therefore lets us realistically gauge how the models actually behave when confronted with messy real world host load measurements. Our study is also unique in focusing on predictability on the scale of seconds as opposed to minutes or longer time scales, and thus is perhaps most useful in the context of interactive applications such as scientific visualization tools~\cite{DV-PRELIM-REPORT-PDPTA99} running on distributed systems. Finally, we used the codebase of a real distributed resource prediction service, RPS~\cite{DINDA-RPS-TR}. RPS is used in the Remos resource measurement system~\cite{REMOS-HPDC98} and in BBN's QuO distributed object quality of service system~\cite{QUO-JOURNAL-98}. We began by choosing to measure host load by the Digital Unix five-second load average. We found that this measure, which can be easily acquired by user-level programs, is closely related to the execution time of short compute-bound tasks (Section~\ref{sec:load_running_time}.) We collected a large number of 1 Hz benchmark load traces, which capture all the dynamics of the load signal, and subjected them to a detailed statistical analysis, which we summarize in Section~\ref{sec:stats}. On the one hand, this analysis suggested that linear time series models such as those in the Box-Jenkins~\cite{BOX-JENKINS-TS} AR, MA, ARMA, and ARIMA classes might be appropriate for predicting load. On the other hand, the existence of self-similarity induced long-range dependence suggested that such models might require an impractical number of parameters or that the much more expensive ARFIMA model class~\cite{HOSKING-FRACT-DIFF,GRANGER-JOYEUX-FRACT-DIFF,BERAN-STAT-LONG-RANGE-METHODS}, which explicitly captures long-range dependence, might be more appropriate. Since it is not obvious which model is best, we empirically evaluated the predictive power of the AR, MA, ARMA, ARIMA, and ARFIMA model classes, as well as a simple ad hoc windowed-mean predictor called BM and a long-term mean predictor called MEAN. We describe these model classes and our implementations of them in Section~\ref{sec:ltsmodels}. Our evaluation methodology, which we describe in detail in Section~\ref{sec:method}, was to run randomized testcases on the benchmark load traces. The testcases (152,000 in all, or about 4000 per trace) were randomized with respect to model class, the number of model parameters and their distribution, the trace subsequence used to fit the model, and the length of the subsequent trace subsequence used to test the model. We collected a large number of metrics for each testcase, but we concentrate on the mean squared (prediction) error metric in this paper. This metric is directly comparable to the raw variance in the load signal, and, with a normality assumption, can be translated into a confidence interval for the running time of a task. Moreover, by the central limit theorem, these estimates of a model class's expected mean squared error are themselves normally distributed, and thus we can determine, to a given significance level, whether one model provides lower expected mean squared error than another. Further, we determine how much the error varies from testcase to testcase. The results of our evaluation are presented in Section~\ref{sec:results}. We found that host load is consistently predictable to a useful degree from past behavior. Except for the MA models, which performed quite badly, the predictive models were all roughly similar, although statistically significant differences were indeed found. These marginal differences do not seem to warrant the use of the more sophisticated model classes because their run-time costs are much higher. Our recommendation is to use AR models of relatively high order (16 or better) for load prediction within the regime we studied. %Section~\ref{sec:related_work} puts our study in the context of other %work in the area, and Section~\ref{sec:conclusions} summarizes and %concludes with a discussion of future directions of this research. %A typical %testcase in our study fit a randomly selected model with a randomly %selected number of up to 8 parameters to a randomly selected five %minute to three hour subsequence of a trace. This fitted model was %then evaluated on an immediately following subsequence with a length %randomly selected from five minutes to three hours. As each %subsequent value in this test sequence became available, the fitted %model was used to make 1,2,3,...,30 second ahead predictions. At the %end of the testcase, the mean squared error for each of the prediction %horizons, as well as a variety of other metrics was recorded. It is %from an aggregate of about 152,000 such testcases (about 4000 per %trace) \Section{Host load and running time} \label{sec:load_running_time} \begin{figure} \centerline{ \colfigwidth \epsfbox{load2time.eps} } \caption{Relationship between average load during execution and execution time} \label{fig:load2time} \end{figure} For CPU-bound tasks, there is an intuitive relationship between host load and execution time. Consider Figure~\ref{fig:load2time}, which plots execution time versus average load experienced during execution for tasks consisting of simple wait loops. The data was generated by running variable numbers of these tasks together simultaneously, at identical priority levels, on an otherwise unloaded Digital Unix machine. Each of these tasks sampled the Unix five-second load average at roughly one second intervals during their execution and at termination printed the average of these samples as well as their execution time. It is these pairs that make up the 42,000 points in the figure. Notice that the relationship between the measured load during execution and the execution time is almost perfectly linear ($R^2>0.99$.) If load were presented as a continuous signal, we would summarize this relationship between execution time and load as \begin{displaymath} \int_{t_{now}}^{t_{now}+t_{exec}}\frac{1}{1+z(t)}dt=t_{nominal} \end{displaymath} where $t_{now}$ is when the task begins executing, $t_{nominal}$ is the execution time of the task on a completely unloaded machine, $(1+z(t))$ is the load experienced during execution ($z(t)$ is the continuous ``background'' load), and $t_{exec}$ is the task's execution time. Notice that the integral simply computes the inverse of the average load during execution. In practice, we can only sample the load with some non-infinitesimal sample period $\Delta$ so we can only approximate the integral by summing over the values in the sample sequence. In this paper, $\Delta=1$ second, so we will write such a sequence of load samples as $\langle z_t \rangle =\ldots,z_{t-1},z_t,z_{t+1},\ldots$ We will also refer to such a sequence as a {\em load signal}. If we knew $z_{t+k}$ for $k>0$, we could directly compute the execution time. Unfortunately, the best we can do is predict these values in some way. The quality of these predictions will then determine how tightly we can bound the execution time of our task. We measure prediction quality by the {\em mean squared error}, which is the average of the square of the difference between predicted values and actual values. Note that there is a mean squared error associated with every {\em lead time} $k$. Two-step-ahead predictions ($k=2$) have a different (and probably higher) mean squared error than one-step-ahead predictions ($k=1$), and so on. During some intervals of time, the load signal may be less predictable than others. Estimates of the mean squared error for those intervals would be larger than for other intervals. Because of this variability, care must be taken in comparing different prediction techniques. However, if we aggregate the estimates for many different intervals, the central limit theorem tells us that the resulting {\em expected mean squared error} is normally distributed, and thus we use this quantity in our comparisons. Intuitively, this is the mean squared prediction error one would expect from a given prediction technique confronted with a randomly chosen interval. For the simplest prediction, the long-term mean of the signal, the mean squared error is simply the variance of the load signal. As we shall see in Section~\ref{sec:stats}, load is highly variable and exhibits other complex properties, which leads to very loose estimates of execution time. The hope is that sophisticated prediction schemes will have much lower mean squared errors. \Section{Statistical properties of load traces} \label{sec:stats} %We collected and studied long, fine grain load traces from a wide %variety of machines in order to discover common statistical properties %and their implications for load prediction. This work, which we %summarize here, is reported in detail in~\cite{DINDA-STAT-PROPS-HOST-LOAD-LCR98}. The %results suggest that classical Box-Jenkins linear models may be %appropriate for load prediction, but that large numbers of parameters %may be necessary to capture its complicated dynamics. Two interesting %new results are that load signals are self-similar and that they %display what we call epochal behavior. Self-similarity suggests the %use of expensive ARFIMA linear models, which can capture long-range %dependence, may be appropriate. Epochal behavior suggests that %reasonably frequent model refitting, piecewise linear models (such as %TAR models~\cite{TAR-TONG-BOOK}), or locally weighted regression on an %``unfolded'' load signal~\cite{ABARBANEL} may be appropriate. This %paper studies the performance of the classical Box-Jenkins models and %the ARFIMA model. In Section~\ref{eval} we evaluate these schemes on %these load traces we introduce in this section. %\SubSection{Measurement methodology} %\label{sec:load_trace_measurement} The load on a Unix system at any given instant is the number of processes that are running or are ready to run, which in turn is the length of the ready queue maintained by the scheduler. The kernel samples the length of the ready queue at some rate and exponentially averages the samples to produce a load average which can be accessed from a user program. It is important to note that while this filtering does tend to correlate the load average over the short term (shorter than the time constant of the filter), the exponential-averaging filter does expose the full spectral content of the underlying signal. The Unix we used was Digital Unix (DUX.) DUX is interesting because the time constant is a mere five seconds, which minimizes the effect of phantom correlation due to the filter. This is especially important when gauging the efficacy of prediction techniques. While the analysis of this section and the prediction results of Section~\ref{sec:results} use the filtered load signal as it is directly available to applications (and which correlates strongly with running time as shown in Section~\ref{sec:load_running_time}), we have also applied our analysis and prediction techniques to the ``unfiltered'' load signal with similar results. By subjecting various DUX systems to varying loads and sampling at progressively higher rates, we determined that DUX updates the load average value at a rate of $1/2$ Hz. We chose to sample at 1 Hz in order to capture all of the load signal's dynamics made available by the kernel. We collected load traces on 38 hosts belonging to our lab %the Computing, Media, and Communication Laboratory (CMCL) %at CMU and to the Pittsburgh Supercomputing Center (PSC) for slightly more than one week in late August, 1997. A second set of week-long traces was acquired on almost exactly the same set of machines in late February through early March, 1998. The results of the statistical analysis were similar for the two sets of traces. In this paper, we describe and use the August 1997 set. A more detailed description of our analysis and the individual results for each trace is available elsewhere~\cite{DINDA-STAT-PROP-HOST-LOAD-EXT-TR} and we will happily supply the traces themselves to interested parties. All of the hosts in the August 1997 set were DEC Alpha machines running DUX. Thirteen of the hosts are in the PSC's production cluster. Of these, there are two front-end machines, four interactive machines, and seven batch machines. In addition, traces were taken on eight hosts that comprise our lab's experimental cluster, two large memory machines used by our group as compute servers, and fifteen desktop workstations owned by members of our research group. \comment{ \begin{figure} \centerline{ \epsfxsize=3.5in \epsfbox{Aug97BoxPlot.eps} } \caption{Statistical summary of load traces.} \label{fig:trace_stats} \end{figure} } \comment{ Figure~\ref{fig:trace_stats} summarizes these traces in a Box plot variant. The central line in each box marks the median load, while the lower and upper lines mark the 25th and 75th percentiles. The lower and upper ``whiskers'' extending from the box mark the actual 2.5th and 97.5th percentiles. The circle marks the mean and the triangles mark the 2.5th and 97.5th percentiles assuming a normal distribution with the trace's mean and variance. } The following points summarize the results of our statistical analysis~\cite{DINDA-STAT-PROP-HOST-LOAD-EXT-TR} that are relevant to this study: (1) The traces exhibit low means but very high variability, measured by the variance, interquartile range, and maximum. Only four traces have mean loads near 1.0. The standard deviation (square root of the variance) is typically at least as large as the mean, while the maximums can be as much as two orders of magnitude larger. This high variability indicates that there exists ample opportunity for prediction algorithms to improve things. (2) Measures of variability, such as the variance and maximum, are positively correlated with the mean, so a machine with a high mean load will also tend to have a large variance and maximum. This correlation suggests that there is more opportunity for prediction algorithms on more heavily loaded machines. (3) The traces have relatively complex, sometimes multimodal distributions that are not well fitted by common analytic distributions. However, we note here that assuming normality (but disallowing negative values) for the purpose of computing a 95\% confidence interval is a reasonable operation. \comment{ , as can be seen by comparing the 2.5th and 97.5th percentile whiskers and the corresponding normal percentiles (triangles) computed from the mean and variance in Figure~\ref{fig:trace_stats}. } (4) Time series analysis of the traces shows that load is strongly correlated over time. The autocorrelation function typically decays very slowly while the periodogram shows a broad, almost noise-like combination of all frequency components. An important implication is that linear models may be appropriate for predicting load signals. However, the complex frequency domain behavior suggests such models may have to be of unreasonably high order. \comment{ \begin{figure} \centerline{ \colfigwidth \epsfbox{hurst_all.epsf} } \caption{Hurst parameter estimates.} \label{fig:trace_hurst} \end{figure} } % % Decided to get rid of self-similarity citation altogether - we cite % the lcr paper earlier % \comment{We estimated their Hurst parameters~\cite{HURST51} using four different nonparametric methods~\cite{TAQQU-HURST-ESTIMATORS-COMPARE} which provided the results in Figure~\ref{fig:trace_hurst}.~\footnote{Note that the estimators we used cannot provide a confidence interval. Estimators that can, such as Whittle's estimators, make assumptions about the system underlying the data (ie, they are parametric) that are unjustified for our data.} } (5) The traces exhibit self-similarity with Hurst Parameters~\cite{FRACTAL-STRUCTURES-CHAOS-SCI-MED, BERAN-STAT-LONG-RANGE-METHODS} ranging from 0.63 to 0.97, with a strong bias toward the top of that range. Hurst parameters in the range of 0.5 to 1.0 indicate self-similarity with positive near-neighbor correlations. This result tells us that load varies in complex ways on all time scales and has long-range dependence. Long-range dependence suggests that using the fractional ARIMA (ARFIMA) modeling approach~\cite{HOSKING-FRACT-DIFF,GRANGER-JOYEUX-FRACT-DIFF,BERAN-STAT-LONG-RANGE-METHODS} may be appropriate. %\begin{figure} %\centerline{ %\colfigwidth %\epsfbox{sdev_mean_epoch.epsf} %} %\caption{Mean epoch lengths and +/- one standard deviation.} %\label{fig:trace_epoch} %\end{figure} (6) The traces display what we term ``epochal behavior.'' The local frequency content (measured by using a spectrogram) of the load signal remains quite stable for long periods of time (150-450 seconds mean, with high standard deviations), but changes abruptly at the boundaries of such epochs. Such abrupt transitions are likely to be due to processes being created, destroyed, or entering new execution phases. This result implies that linear models may have to be refit at these boundaries. Our evaluation (Section~\ref{sec:results}) ignores these boundaries, so we do see the effect of inadvertently crossing one during prediction. Strictly speaking, (6) means that load is not stationary. However, it is also not free to wander at will --- clearly load cannot rise to infinite levels or fall below zero. This is not incompatible with the ``borderline stationarity'' implied by (5). %Figure~\ref{fig:trace_epoch} shows the mean epoch %length and standard deviation of epoch length for each of the traces. \Section{Linear time series models} \label{sec:ltsmodels} \begin{figure} \centerline { \epsfxsize=3in \epsfbox{ltsmodel.eps} } \caption{Linear time series models.} \label{fig:ltsmodels} \end{figure} The main idea behind using a linear time series model in load prediction is to treat the sequence of periodic samples of host load, $\langle z_t \rangle$, as a realization of a stochastic process that can be modeled as a white noise source driving a linear filter. The filter coefficients can be estimated from past observations of the sequence. If most of the variability of the sequence results from the action of the filter, we can use its coefficients to estimate future sequence values with low mean squared error. Figure~\ref{fig:ltsmodels} illustrates this decomposition. In keeping with the relatively standard Box-Jenkins notation~\cite{BOX-JENKINS-TS}, we represent the input white noise sequence as $\langle a_t \rangle$ and the output load sequence as $\langle z_t \rangle$. On the right of Figure~\ref{fig:ltsmodels} we see our partially predictable sequence $\langle z_t \rangle$, which exhibits some mean $\mu$ and variance $\sigma_z^2$. On the left, we see our utterly unpredictable white noise sequence $\langle a_t \rangle$, which exhibits a zero mean and a variance $\sigma_a^2$. In the middle, we have our fixed linear filter with coefficients $\langle \psi_j \rangle$. Each output value $z_t$ is the sum of the current noise input $a_t$ and all previous noise inputs, weighted by the $\langle\psi_j\rangle$ coefficients. Given an observed load sequence $\langle z_t \rangle$, the optimum values for the coefficients $\psi_j$ are those that minimize $\sigma_a^2$, the variance of the driving white noise sequence $\langle a_t \rangle$. Notice that the one-step-ahead prediction given all the data up to and including time $t-1$ is $\hat{z}_{t-1}^1 =\sum_{j=1}^\infty \psi_ja_{t-j}$, since the expected value of $a_{t}=0$.) The noise sequence consists simply of the one-step-ahead prediction errors and the optimal coefficient values minimize the sum of squares of these prediction errors. The general form of the linear time series model is, of course, impractical, since it involves an infinite summation using an infinite number of completely independent weights. Practical linear time series models use a small number of coefficients to represent infinite summations with restrictions on the weights, as well as special casing the mean value of the sequence, $\mu$. To understand these models, it is easiest to represent the weighted summation as a ratio of polynomials in $B$, the backshift operator, where $B^dz_t =z_{t-d}$. For example, we can write $z_t = \sum_{j=1}^\infty \psi_ja_{t-j}+a_t$ as $z_t = \psi(B)a_t$ where $\psi(B)=1+\psi_1B+\psi_2B^2+\ldots$ Using this scheme, the models we examine in this paper can be represented as \begin{equation} z_t = \frac{\theta(B)}{\phi(B)(1-B)^d}a_t + \mu \label{eqn:ts_poly} \end{equation} where the different model classes we examine in this paper (AR, MA, ARMA, ARIMA, ARFIMA, BM, and MEAN) constrain $\theta(B)$, $\phi(B)$ and $d$ in different ways. In the signal processing domain, this kind of filter is known as a pole-zero filter. The roots of $\theta(B)$ are the zeros and the roots of $\phi(B)(1-B)^d$ are the poles. It is also a state-space filter, as can be seen if written as $z_t + \eta_1z_{t-1} + \eta_2z_{t-2} + \ldots + \eta_{p+d}z_{t-p-d} = a_t + \theta_1a_{t-1}+\theta_2a_{t-2} + \ldots + \theta_qz_{t-q}$ where $\eta(B)=\phi(B)(1-B)^d$ for integer $d$. In general, such a filter can be unstable in that its outputs can rapidly diverge from the input signal. This instability is extremely important from the point of view of the implementor of a load prediction system. Such a system will generally fit a model (choose the $\theta(B)$ and $\phi(B)$ coefficients, and $d$) using some $m$ previous observations. The model will then be ``fed'' the next $n$ observations and asked to make predictions in the process. If the coefficients are such that the filter is unstable, then it may explain the initial $m$ observations very well, yet fail miserably and even diverge (and crash!) when used on the $n$ observations after the fitting. {\bf AR($p$) models:} The class of AR($p$) models (purely autoregressive models) has $z_t = \frac{1}{\phi(B)}a_t + \mu$ where $\phi(B)$ has $p$ coefficients. From the point of view of a system designer, AR($p$) models are highly desirable since they can be fit to data in a deterministic amount of time. In the Yule-Walker technique that we used, the autocorrelation function is computed to a maximum lag of $p$ and then a $p$-wide Toeplitz system of linear equations is solved. Even for relatively large values of $p$, this can be done almost instantaneously. {\bf MA($q$) models:} The class of MA($q$) models (purely moving average models) has $z_t = \theta(B)a_t$ where $\theta(B)$ has $q$ coefficients. MA($q$) models are a much more difficult proposition for a system designer since fitting them takes a nondeterministic amount of time. Instead of a linear system, fitting a MA($q$) model presents us with a quadratic system. Our implementation, which is nonparametric (ie, it assumes no specific distribution for the white noise source), uses the Powell procedure to minimize the sum of squares of the $t+1$ prediction errors. The number of iterations necessary to converge is nondeterministic and data dependent. {\bf ARMA($p$,$q$) models:} The class of ARMA($p$,$q$) models (autoregressive moving average models) has $z_t=\frac{\theta(B)}{\phi(B)}a_t + \mu$ where $\phi(B)$ has $p$ coefficients and $\theta(B)$ has $q$ coefficients. By combining the AR($p$) and MA($q$) models, ARMA($p$,$q$) models hope to achieve greater parsimony --- using fewer coefficients to explain the same sequence. From a system designer's point of view, this may be important, at least in so far as it may be possible to fit a more parsimonious model more quickly. Like MA($q$) models, however, ARMA($p$,$q$) models take a nondeterministic amount of time to fit to data, and we use the same Powell minimization procedure to fit them. {\bf ARIMA($p$,$d$,$q$) models:} The class of ARIMA($p$,$d$,$q$) models (autoregressive integrated moving average models) implement Equation~\ref{eqn:ts_poly} for positive integer $d$. Intuitively, the $(1-B)^d$ component amounts to a $d$-fold integration of the output of an ARMA($p$,$q$) model. Although this makes the filter inherently unstable, it allows for modeling nonstationary sequences. Such sequences can vary over an infinite range and have no natural mean. Although load clearly cannot vary infinitely, it doesn't have a natural mean either. ARIMA($p$,$d$,$q$) models are fit by differencing the sequence $d$ times and fitting an ARMA($p$,$q$) model as above to the result. {\bf ARFIMA($p$,$d$,$q$) models:} The class of ARFIMA($p$,$d$,$q$) models (autoregressive fractionally integrated moving average models) implement Equation~\ref{eqn:ts_poly} for fractional values of $d$, $0 hatz_tcross+k^i % compare hatz_tcross+k^i to ztcorss+k+i To explore this space in a reasonably unbiased way, we ran a large number of randomized testcases on each of the traces. Figure~\ref{fig:testcasegen} illustrates the parameter space from which testcase parameters are chosen. A testcase is generated and evaluated using the following steps. \small \begin{ecompact} \item Choose a random {\em crossover point}, $t_{cross}$, from within the trace. \item Choose a random number of samples, $m$, from $600,601,\ldots,10800$ (5 minutes to three hours.) The $m$ samples preceding the crossover point, $z_{t_{cross}-m}, z_{t_{cross}-m+1}, \ldots,z_{t_{cross}-1}$ are in the {\em fit interval}. \item Choose a random number of samples, $n$ from $600,601,\ldots,10800$ (5 minutes to three hours.) The $n$ samples including and following the crossover point, $z_{t_{cross}}, z_{t_{cross}+1}, \ldots, z_{t_{cross}+n-1}$,are in the {\em test interval}. \item Choose a random AR, MA, ARMA, ARIMA, or ARFIMA {\em test model} from the table in Figure~\ref{fig:testcasegen}, fit it to the samples in the fit interval, and generate a predictor from the fitted test model. \item For $i=m$ to $1$, step the predictor with $z_{t_{cross}-i}$ (the values in the fit interval) to initialize its internal state. After this step, the predictor is ready to be tested. \item For $i=0$ to $n-1$ do the following: \begin{icompact} \item Step the predictor with $z_{t_{cross}+i}$ (the next value in the test interval.) \item For each lead time $k=1,2,\ldots,30$ seconds, produce the predictions $\hat{z}_{t_{cross}+i}^k$. $\hat{z}_{t_{cross}+i}^k$ is the prediction of $z_{t_{cross}+i+k}$ given the samples $z_{t_{cross}-m},z_{t_{cross}-m+1},\ldots,z_{t_{cross}+i}$. Compute the prediction errors $a_{t+i}^k = \hat{z}_{t_{cross}+i}^k -z_{t_{cross}+i+k}$. \end{icompact} \item For each lead time $k=1,2,\ldots,30$ seconds, analyze the $k$-step-ahead prediction errors $a_{t+i}^k$, $i=0,1,\ldots,n-1$. \item Output the testcase parameters and the analysis of the prediction errors. \end{ecompact} \normalsize For clarity in the above, we focused on the linear time series model under test. Each testcase also includes a parallel evaluation of the BM and MEAN models in order to facilitate direct comparison with the simple BM model and the raw signal variance. The lower limit we place on the length of the fit and test intervals is purely prosaic --- the ARFIMA model needs about this much data to be successfully fit. The upper limit is chosen to be greater than most epoch lengths so that we can see the effect of crossing epoch boundaries. The models are limited to eight parameters because fitting larger MA, ARMA, ARIMA, or ARFIMA models is prohibitively expensive in a real system. We did also explore larger AR models, up to order 32. The analysis of the prediction errors includes the following. For each lead time, the minimum, median, maximum, mean, mean absolute, and mean squared prediction errors are computed. The one-step-ahead prediction errors (ie, $a_{t+i}^1$, $i=0,1,\ldots,n-1$) are also subject to IID and normality tests as described by Brockwell and Davis~\cite{INTRO-TIME-SERIES-FORECASTING}, pp. 34--37. IID tests included the fraction of the autocorrelations that are significant, the Portmanteau Q statistic (the power of the autocorrelation function), the turning point test, and the sign test. Normality was tested by computing the $R^2$ value of a least-squares fit to a quantile-quantile plot of the values or errors versus a sequence of normals of the same mean and variance. We ran approximately 152,000 such testcases, which amounted to about 4000 testcases per trace, or about 1000 per model class and parameter set, or about 30 per trace, model class and parameter set. Our parallelized simulation discarded testcases in which an unstable model ``blew up,'' either detectably or due to a floating point exception. The results of the accepted testcases were committed to a SQL database to simplify the analysis discussed in the following section. %% end of part that needs to be formalized \Section{Results} \label{sec:results} The section addresses the following questions: Is load consistently predictable? If so, what are the consistent differences between the different model classes, and which class is preferable? To answer these questions we analyze the database of randomized testcases from Section~\ref{sec:method}. For the most part we will address only the mean squared error results, although we will touch on the other results as well. It is important to reiterate the comments of Section~\ref{sec:method}: we are examining the practical predictive power of the models here, not the their explanatory power, or ``fit.'' {\bf Load is consistently predictable:} For a model to provide consistent predictability of load, it must satisfy two requirements. First, for the average testcase, the model must have a considerably lower expected mean squared error than the expected raw variance of the load signal (ie, the expected mean squared error of the MEAN model.) The second requirement is that this expectation is also very likely, or that there is little variability from testcase to testcase. Intuitively, the first requirement says that the model provides good predictions on average, while the second says that most predictions are close to that average. %\begin{figure} %\centerline { %\colfigwidth %\epsfbox{all_lead1_8to8.eps} %} %\caption{All traces, 1 second lead, 8 parameters} %\label{fig:all_lead1_8to8} %\end{figure} \begin{figure} \centerline { \begin{tabular}{c} \colfigwidth \epsfbox{all_lead1_8to8.eps} \\ (a) 1 second predictions \\ \colfigwidth \epsfbox{all_lead15_8to8.eps} \\ (b) 15 second predictions \\ \colfigwidth \epsfbox{all_lead30_8to8.eps} \\ (c) 30 second predictions \\ \end{tabular} } \caption{Distributions of mean squared errors for all traces using 8 parameter models.} \label{fig:all_8to8} \end{figure} Figure~\ref{fig:all_8to8}(a) suggests that load is indeed consistently predictable in this sense. The figure is a Box plot that shows the distribution of one-step-ahead (one second) mean squared error measures (ie, the distribution of the measure $\frac{1}{n}\sum_{i=0}^{n-1}(a_{t+i}^{1})^2$, using the formalisms of Section~\ref{sec:method}) for 8 parameter models on all of the traces. In the figure, each category is a specific class of model and is annotated with the number of samples for that class. For each class, the circle indicates the expected mean squared error, while the triangles indicated the 2.5th and 97.5th percentiles assuming a normal distribution. The center line of each box shows the median while the lower and upper limits of the box show the 25th and 75th percentiles and the lower and upper whiskers show the actual 2.5th and 97.5th percentiles. Notice that the expected raw variance (MEAN) of a testcase is approximately 0.05, while the expected mean squared error for {\em all} of the model classes is nearly zero. In terms of the length of the 95\% confidence interval for the execution time of a one second task as described in Section~\ref{sec:load_running_time}, this is a reduction from $(2)(1.96)\sqrt{0.05} = 0.87$ seconds to virtually zero for all of the classes of predictive models, including the simple BM model. The figure also shows that our second requirement for consistent predictability is met. We see that the variability around the expected mean squared error is much lower for the predictive models than for MEAN. For example, the 97.5th percentile of the raw variance is almost 0.3 (2.2 second interval) while it is about 0.02 (0.6 second interval) for the predictive models. %\begin{figure} %\centerline { %\colfigwidth %\epsfbox{all_lead15_8to8.eps} %} %\caption{All traces, 15 second lead, 8 parameters} %\label{fig:all_lead15_8to8} %\end{figure} %\begin{figure} %\centerline { %\colfigwidth %\epsfbox{all_lead30_8to8.eps} %} %\caption{All traces, 30 second lead, 8 parameters} %\label{fig:all_lead30_8to8} %\end{figure} Figures~\ref{fig:all_8to8}(b) and~\ref{fig:all_8to8}(c) show the results for 15 second predictions and 30 second predictions. Notice that, except for the MA models, the predictive models are consistently better than the raw load variance, even with 30 second ahead predictions. We also see that MA models perform quite badly, especially at higher lead times. This was also the case when we considered the traces individually and broadened the number of parameters. MA models are clearly ineffective for load prediction. {\bf Successful models have similar performance:} Surprisingly, Figures~\ref{fig:all_8to8}(a)--(c) also show that the differences between the successful models are actually quite small. This is also the case if we expand to include testcases with 2 to 8 parameters instead of just 8 parameters. With longer lead times, the differences do slowly increase. \begin{figure} \centerline { \begin{tabular}{c} \colfigwidth \epsfbox{axp0_lead1_8to8.eps} \\ (a) 1 second predictions \\ \colfigwidth \epsfbox{axp0_lead15_8to8.eps} \\ (b) 15 second predictions \\ \colfigwidth \epsfbox{axp0_lead30_8to8.eps} \\ (c) 30 second predictions \\ \end{tabular} } \caption{Distributions of mean squared errors for axp0.psc trace using 8 parameter models.} \label{fig:axp0_8to8} \end{figure} %\begin{figure} %\centerline { %\colfigwidth %\epsfbox{axp0_lead1_8to8.eps} %} %\caption{Axp0, 1 second lead, 8 parameters} %\label{fig:axp0_lead1_8to8} %\end{figure} %\begin{figure} %\centerline { %\colfigwidth %\epsfbox{axp0_lead15_8to8.eps} %} %\caption{Axp0, 15 second lead, 8 parameters} %\label{fig:axp0_lead15_8to8} %\end{figure} %\begin{figure} %\centerline { %\colfigwidth %\epsfbox{axp0_lead30_8to8.eps} %} %\caption{Axp0, 30 second lead, 8 parameters} %\label{fig:axp0_lead30_8to8} %\end{figure} For more heavily loaded machines, the differences can be much more dramatic. For example, Figure~\ref{fig:axp0_8to8}(a) shows the distribution of one-step-ahead (one second) mean squared error measures for 8 parameter models on the axp0.psc trace. Here we see an expected raw variance (MEAN) of almost 0.3 (2.2 second confidence interval) reduced to about 0.02 (0.6 second interval) by nearly all of the models. Furthermore, the mean squared errors for the different model classes are tightly clustered around the expected 0.02 value, quite unlike with MEAN, where we can see a broad range of values and the 97.5th percentile is almost 0.5 (2.8 second interval.) The axp0.psc trace and others are also quite amenable to prediction with long lead times. For example, Figures~\ref{fig:axp0_8to8}(b) and (c) show 15 and 30 second ahead predictions for 8 parameter models on the axp0.psc trace, respectively. With the exception of the MA models, even 30 second ahead predictions are consistently much better than the raw signal variance. These figures remain essentially the same if we include testcases with 2 to 8 parameters instead of just 8 parameters. \begin{figure}[t] \centerline{ \begin{tabular}{c} \figfontsize \begin{tabular}{l|ccccccc} &MEAN &BM &AR &MA &ARMA &ARIMA &ARFIMA\\ \hline MEAN &$=$ &$>$ &$>$ &$>$ &$>$ &$>$ &$>$ \\ BM &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ \hline AR &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ \hline MA &$<$ &$>$ &$>$ &$=$ &$>$ &$=$ &$>$ \\ ARMA &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ ARIMA &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ ARFIMA &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ \end{tabular} \\ \normalsize (a) all traces, one step ahead \\ \figfontsize \begin{tabular}{l|ccccccc} &MEAN &BM &AR &MA &ARMA &ARIMA &ARFIMA \\ \hline MEAN &$=$ & $>$ & $>$ & $=$ & $>$ & $>$ & $>$\\ BM &$<$ &$=$ & $>$ & $<$ & $=$ & $=$ & $>$\\ \hline AR &$<$ &$<$ & $=$ & $<$ & $=$ & $<$ & $=$\\ \hline MA &$=$ &$>$ & $>$ & $=$ & $>$ & $>$ & $>$\\ ARMA &$<$ &$=$ & $=$ & $<$ & $=$ & $=$ & $>$\\ ARIMA &$<$ &$=$ & $>$ & $<$ & $=$ & $=$ & $>$\\ ARFIMA &$<$ &$<$ & $=$ & $<$ & $<$ & $<$ & $=$\\ \end{tabular} \\ \normalsize (b) axp0.psc trace, 16 steps ahead \end{tabular} } \caption{T-test comparisons: (a) all traces aggregated, lead 1, 8 parameters. Longer leads are essentially the same except MA is always worse. (b) axp0, lead 16, 8 parameters.} \label{fig:t-test-8to8} \normalsize \end{figure} %\begin{figure}[t] %\figfontsize %\centerline{ %\begin{tabular}{l|ccccccc} % &MEAN &BM &AR &MA &ARMA &ARIMA &ARFIMA\\ %\hline %MEAN &$=$ &$>$ &$>$ &$>$ &$>$ &$>$ &$>$ \\ %BM &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ %\hline %AR &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ %\hline %MA &$<$ &$>$ &$>$ &$=$ &$>$ &$=$ &$>$ \\ %ARMA &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ %ARIMA &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ \\ %ARFIMA &$<$ &$=$ &$=$ &$<$ &$=$ &$=$ &$=$ %\end{tabular} %} %\caption{T-test comparisons - all traces aggregated, lead 1, 8 %parameters. Longer leads are essentially the same except MA is always %worse.} %\label{fig:t-test-all-lead1-8to8} %\normalsize %\end{figure} Although the differences in performance between the successful models are very small, they are generally statistically significant. We can algorithmically compare the expected mean squared error of the models using the unpaired t-test~\cite{JAIN-PERF-ANALYSIS}, pp. 209--212, and do ANOVA procedures to verify that the differences we detect are significantly above the noise floor. For each pair of model classes, the t-test tells us, with 95\% confidence, whether the first model is better, the same, or worse than the second. We do the comparisons for the cross product of the models at a number of different lead times. We consider the traces both in aggregate and individually, and we use several different constraints on the number of parameters. Figure~\ref{fig:t-test-8to8}(a) shows the results of such a t-test comparison for the aggregated traces, a lead time of 1, and 8 parameters. In the figure, the row class is being compared to the column class. For example, at the intersection of the AR row and MA column, there is a '$<$', which indicates that, with 95\% confidence, the expected mean squared error of the AR models is less than that of MA models, thus the AR models are better than MA models for this set of constraints. For longer lead times, we found that the results remained essentially the same, except that MA became consistently worse. The message of Figure~\ref{fig:t-test-8to8}(a) is that, for the typical trace and a sufficient number of parameters, there are essentially no differences in the expected mean squared error of the BM, AR, ARMA, ARIMA, and ARFIMA models, even with high lead times. Further, with the exception of the MA models, all of the models are better than the raw variance of the signal (MEAN model.) %\begin{figure}[t] %\figfontsize %\centerline{ %\begin{tabular}{l|ccccccc} % &MEAN &BM &AR &MA &ARMA &ARIMA &ARFIMA \\ %\hline %MEAN &$=$ & $>$ & $>$ & $=$ & $>$ & $>$ & $>$\\ %BM &$<$ &$=$ & $>$ & $<$ & $=$ & $=$ & $>$\\ %\hline %AR &$<$ &$<$ & $=$ & $<$ & $=$ & $<$ & $=$\\ %\hline %MA &$=$ &$>$ & $>$ & $=$ & $>$ & $>$ & $>$\\ %ARMA &$<$ &$=$ & $=$ & $<$ & $=$ & $=$ & $>$\\ %ARIMA &$<$ &$=$ & $>$ & $<$ & $=$ & $=$ & $>$\\ %ARFIMA &$<$ &$<$ & $=$ & $<$ & $<$ & $<$ & %$=$\\ %\end{tabular} %} %\caption{T-test comparisons - axp0, lead 16, 8 parameters.} %\label{fig:t-test-axp0-lead16-8to8} %\normalsize %\end{figure} For machines with higher mean loads, there are more statistically significant differences between the models. For example, Figure~\ref{fig:t-test-8to8}(b) shows a t-test comparison for the axp0.psc trace at a lead time of 16 seconds for 8 parameter models. Clearly, there is more differentiation here, and we also see that the ARFIMA models do particularly well, as we might expect given the self-similarity result of Section~\ref{sec:stats}. However, notice that the AR models are doing about as well. The results of these t-test comparisons can be summarized as follows: (1) Except for the MA models, the models we tested are significantly better (in a statistical sense) than the raw variance of the trace and the difference in expected performance is significant from a systems point of view. (2) The AR, ARMA, ARIMA, and ARFIMA models are significantly better than or equal to the BM model and occasionally the difference in expected performance is also significant from a systems point of view. (3) The AR, ARMA, ARIMA, and ARFIMA models have similar expected performance. {\bf AR models are appropriate:} Clearly, using one of the model classes, other than MA, for load prediction is beneficial. Further, using an AR, ARMA, ARIMA, or ARFIMA model is occasionally preferable to the BM model. Since there are no major differences in the performance of these models, and since fitting AR models is much less expensive than the other models (and also takes deterministic amount of time), we recommend the use of AR models for load prediction. AR models of order 4 or higher seem to work fine. However, we found the knee in the performance of the AR models to be around order 16. Since fitting AR models of this order, or even considerably higher, is quite fast (2 ms for an AR($16$) fit to 2000 data points on a 500 MHz Alpha~\cite{DINDA-RPS-TR}), we recommend using AR($16$)s or better for load prediction. % % This section will probably have to go due to space constraints % {\bf Prediction errors are not IID normal:} As we described in Section~\ref{sec:method}, our evaluation of each testcase includes tests for the independence and normality of prediction errors. Intuitively, we want the errors to be independent so that they can be characterized with a probability distribution function, and we want the distribution function to be normal in order to simplify computing confidence intervals from it. For the most part, the errors we see with the different models are not independent or normal to any confidence level. However, from a practical point of view, the errors are much less correlated over time than the raw signal. For example, for AR($8$) models, the Portmanteau Q statistic tends to be reduced by an order of magnitude or more, which suggests that those autocorrelations that are large enough to be significant are only marginally so. Furthermore, assuming normality for computing confidence intervals for high confidence levels, such as 95\%, seems to be reasonable. However, since a load prediction system will probably continuously evaluate a fitted model in order to determine when to refit it, it seems reasonable for it to keep a histogram or other more detailed representation of the errors in order to more accurately compute confidence intervals. %\begin{figure} %\figfontsize %\centerline{ %\begin{tabular}{l|ccccccc} %Model & 0 & 1& 2& 3& 4& 5& 6 \\ %\hline %MEAN & 0.00\% & 0.00\% & 0.00\% & 0.00\% & 0.00\% & 0.00\% & 100.00\% \\ %BM & 13.79\% & 17.24\% & 3.45\% & 24.14\% & 27.59\% & 13.79\% & 0.00\% \\ %AR & 20.69\% & 24.14\% & 17.24\% & 13.79\% & 10.34\% & 13.79\% & 0.00\% \\ %MA & 3.45\% & 6.90\% & 3.45\% & 10.34\% & 31.03\% & 44.83\% & 0.00\% \\ %ARMA & 6.90\% & 20.69\% & 24.14\% & 34.48\% & 10.34\% & 3.45\% & 0.00\% \\ %ARIMA & 48.28\% & 6.90\% & 20.69\% & 3.45\% & 6.90\% & 13.79\% & 0.00\% \\ %ARFIMA & 6.90\% & 24.14\% & 31.03\% & 13.79\% & 13.79\% & 10.34\% & 0.00\% \\ %\end{tabular} %} %\caption{Model rankings by mean for lead 1, 8:8 params} %\label{fig:t-test-example} %\normalsize %\end{figure} %\begin{figure} %\figfontsize %\centerline{ %\begin{tabular}{l|ccccccc} %Model & 0 & 1& 2& 3& 4& 5& 6 \\ %\hline %MEAN & 2.63\% & 5.26\% & 7.89\% & 26.32\% & 7.89\% & 50.00\% & 0.00\% \\ %BM & 5.26\% & 23.68\% & 18.42\% & 15.79\% & 18.42\% & 10.53\% & 7.89\% \\ %AR & 21.05\% & 13.16\% & 13.16\% & 18.42\% & 23.68\% & 7.89\% & 2.63\% \\ %MA & 2.63\% & 7.89\% & 5.26\% & 5.26\% & 2.63\% & 7.89\% & 68.42\% \\ %ARMA & 5.26\% & 21.05\% & 15.79\% & 21.05\% & 18.42\% & 10.53\% & 7.89\% \\ %ARIMA & 31.58\% & 15.79\% & 18.42\% & 7.89\% & 10.53\% & 7.89\% & 7.89\% \\ %ARFIMA & 31.58\% & 13.16\% & 21.05\% & 5.26\% & 18.42\% & 5.26\% & 5.26\% \\ %\end{tabular} %} %\caption{Model rankings by mean for lead 16, 8:8 params} %\label{fig:t-test-example} %\normalsize %\end{figure} \Section{Related work} \label{sec:related_work} In application-centric scheduling~\cite{APPLES-HPDC96}, applications schedule themselves, adapting to the availability of resources, perhaps using a framework such as QuO~\cite{QUO-JOURNAL-98} or Dv~\cite{DV-PRELIM-REPORT-PDPTA99}. Resource monitoring systems such as Remos~\cite{REMOS-HPDC98}, the Network Weather Service~\cite{FORECASTING-NETWORK-NWS-WOLSKI-HPDC97}, or Topology-d~\cite{SERVICE-NET-AWARE-APPS-SPDT-98} provide measurements and predictions to help applications make scheduling decisions. This paper characterized one such measurement (the Unix load average) and studied how to predict it. While considerable effort has gone into characterizing workloads~\cite{LIMITED-BENEFIT-MIGRATION-EAGER,LOAD-BAL-HEURISTICS-PROCESS-BEHAVIOR-LELAND-OTT,EXPLOIT-PROC-LIFETIME-DIST-DYN-LB-SIGMETRICS}, the focus has been on load sharing systems~\cite{ADAPTIVE-LOAD-SHARING-HOMOGENEOUS-EAGER}, which schedule all of jobs in a distributed system, and not on application-centric scheduling. An important assumption in load sharing and balancing systems is that current load is a predictor of future load. This work shows that that assumption is valid in a statistically rigorous manner, and also explores the predictive power of more sophisticated predictive models. Mutka and Livny~\cite{WORKSTATION-AVAIL-STATS-CONDOR} studied workstation availability as a binary function of load average and other measures. Our previous paper~\cite{DINDA-STAT-PROP-HOST-LOAD-EXT-TR} and the summary in this paper are the first detailed study of the Unix load average we are aware of. Although linear time series methods are widely used in other areas, including networking~\cite{TIME-SERIES-MODELS-INTERNET-TRAFFIC-BASU-GT-TR-95,TIME-SERIES-MODEL-LONG-TERM-NSFNET-GROSCHWITZ-ICC94}, little work has been done in using time series methods for host load prediction. Samadani and Kalthofen's work~\cite{PRED-BASED-SCHED-DIST-COMP-SAMADANI-UNPUB-96} is the closest to ours. They found that small ARIMA models were preferable to single-point predictors and Bayesian predictors for predicting load. Their empirical study concentrated on coarse grain prediction (one minute sampling interval) of four traces that measured load as the number of non-idle processes shown by the Unix ``ps'' command. In contrast, we studied finer grain (one second sampling interval) prediction of the DUX five-second load average on a much larger set of machines using higher order models as well as the ARFIMA class of models. Additionally, our study was randomized with respect to models instead of using the Box-Jenkins heuristic identification procedure. Finally, we reached a different conclusion for our regime, namely that AR models are sufficient. The Network Weather Service (NWS) uses windowed mean, median, and AR models to predict various resource measures~\cite{FORECASTING-NETWORK-NWS-WOLSKI-HPDC97} including CPU availability~\cite{WOLSKI-PRED-CPU-AVAIL-HPDC-99}. Our study and the study described in the latter paper are quite complementary. The NWS study confirms our earlier self-similarity results, quantifies how well various measures of CPU availability, including the Unix one-minute load average, correlate with running time, proposes a new hybrid CPU availability measure, and studies the performance of NWS's predictive models for 10 second and five minute predictions of CPU availability measured at a 0.1 Hz rate. In contrast, we evaluated the performance of AR and more complex linear time series models (including one that captures the long-range dependence that self-similarity induces) applied to shorter range prediction (1 to 30 seconds) of a more dynamic load signal (the DUX five-second load average measured at 1 Hz.) In addition, our evaluation used a randomized methodology and a much larger set of traces to gauge the models' predictive power independent of any particular system. Interestingly, both studies reach the conclusion that relatively simple predictive models such as AR are adequate for host load prediction. %Hailperin~\cite{LOAD-BAL-TIME-SERIES-STAT-PERIODIC-LOADS-HAILPERIN} %used the statistical periodicity of his sensor-driven workload to %predict remote load in his load balancing system. % %The AR, MA, ARMA, and ARIMA models have long histories dating back to %the 19th century. Box, et al~\cite{BOX-JENKINS-TS} is the standard %reference and Brockwell and Davis\cite{INTRO-TIME-SERIES-FORECASTING} %provide a useful example-driven approach. ARFIMA models are a %relatively recent %innovation~\cite{HOSKING-FRACT-DIFF,GRANGER-JOYEUX-FRACT-DIFF}. %Bassingthwaighte, et al~\cite{FRACTAL-STRUCTURES-CHAOS-SCI-MED} %provide a good intuitive introduction to self-similarity and %Beran~\cite{BERAN-STAT-LONG-RANGE-METHODS} provides a good %introduction to long-range dependence. \Section{Conclusions and future work} \label{sec:conclusions} We have presented a detailed evaluation of the performance of linear time series models for predicting the Unix five-second load average on a host machine, from 1 to 30 second into the future, using 1 Hz samples of its history. Predicting this load signal is interesting because it is directly related to the running time of compute-bound tasks. We began by studying the statistical properties of week-long 1 Hz load traces collected on 38 different machines. This study suggested that Box-Jenkins (AR, MA, ARMA, and ARIMA) and ARFIMA models might be appropriate. We evaluated the performance of these classes of models and a simple windowed-mean predictor by running 152,000 randomized testcases on the load traces and then analyzed the results. The main contribution of our evaluation is to show, in a rigorous manner, that host load on real systems is predictable to a very useful degree from past behavior by using linear time series techniques. In addition, we discovered that, while there are statistically significant differences between the different classes of models we studied, the marginal benefits of the more complex models do not warrant their much higher run-time costs. We reached the conclusion that AR models of order 16 or higher are sufficient for predicting 1 Hz data up to 30 seconds in the future. These models work very well and are very inexpensive to fit to data and to use. Host load signals are not generated by linear systems. Other techniques, such as threshold autoregressive models or a methodology based on chaotic dynamics may be more appropriate for generating or understanding load signals. However, for prediction, we are satisfied with simple linear methods in the regime we studied in this paper. We are currently working on using predictions produced by RPS in scheduling the soft real-time tasks of a distributed visualization application~\cite{DINDA-CASE-BERT-WPDRTS-99,DV-PRELIM-REPORT-PDPTA99}. \small \tinyspacing \bibliography{texbib} \normalsize \normspacing \end{document}