Modeling Multi-Attribute Demand for Sustainable Cloud Computing with Copulae

1 Introduction

In 2013, US data centers used about 91 billion kWh [NRDC, 2013], showing a 50% increase since 2006. This number is 2.2% of the 4 trillion kWh consumed in total, and is equal to the amount of energy used by 5.6 million average US households [EIA, 2015]. This consumption growth is expected to continue, with another 50% increase in the next five years [NRDC, 2013]. This will require construction of an additional 16 500 MW power plants, and result in 100 million metric tons of CO2 pollution per year. While the need for such facilities will not abate, much can be done to reduce the energy required. Recent work has shown that considerable energy can be saved by using heterogeneous processor types in data centers [Guevara et al., 2014], and then assigning processes to be run on the most energy-efficient resources for handling them by sophisticated resource allocation methods

Benjamin Lubin

School of Management Boston University blubin@bu.edu

[Lubin et al., 2009; Beloglazov and Buyya, 2010; Beloglazov et al., 2012].

However, such methods are only as good as the demand models that drive them, and the models in the literature handle only single resource types, not the full joint distribution of demand. Having high fidelity demand models is critical if the enormous rise in data center energy usage is to be curtailed.

Analyzing the characteristics of workloads is a longstudied problem with application to load balancing, system availability, failure, and process behavior. The literature covering this topic originates in the 1960s. In some of the earliest such studies, Rosin [1965] and Walter and Wallace [1967] study data for more than 10, 000 jobs run at the University of Michigan Computing Center. They find that the shape of CPU process time is skewed and used a simple exponential distribution to model it. Such work continued into the 1970s, for example Agrawala et al. [1976] presented three simple models of system workload. A large literature in the high performance computing literature has also addressed demand modeling [Li et al., 2004].

Most of this work has focused on one specific attribute of demand (e.g., CPU time, memory usage, or I/O operations). To the best of our knowledge there is no existing work on the joint distribution of resource usage. We believe that there is important structure in the correlation between different attributes of demand, and that therefore an understanding the nature of the joint distribution will significantly improve allocation, scheduling and pricing algorithms among other uses for such models, which in turn will yield a significant positive environmental impact due to reduced energy usage.

Accordingly, in this paper, we develop a multivariate model for computational resource usage. Moreover, as part of our model we propose to use a more flexible distribution, the Burr XII, as our model of univariate demand. We will show that this model has significantly more power to capture real-world usage patterns than the models used to date in the literature. Having modeled the individual attributes, we then combine them into a multi-attribute model through the machinery of copulae [Sklar, 1959].

Copulae are an elegant way to build and fit complex joint distributions that have been used to great effect in finance, but which have not been used in computational demand modeling. By employing this method we have been able to accurately handle an otherwise-intractable modeling problem. In this work, we focus on the bivariate case, specifically capturing the joint demand for CPU and memory, typically the two most important resource types. We then validate our approach based on a public data set from Google [Reiss et al., 2011] that contains 29 days of workload data from a large production cluster. Using this data, we examine a large number of models from the family we propose, and demonstrate the effectiveness of a particular model that combines Burr XII marginals with the Frank copulae. The proposed model can be directly generalized to more resource types, though we are limited by our data in doing so here. In this paper we use machine learning and computational statistics to model data center usage, offering the following contributions:

• A novel application of the Burr XII distribution to modeling the distributions of individual workload attributes (Section 2.2).
• The first demand models that capture the joint distribution over multiple attribute types, built by applying copula methods to our marginal models (Section 2.3).
• An application of these new tools to a large Google data set (Section 3).
• Validation of this approach by employing a parametric bootstrapped Kolmogorov-Smirnov test on both the univariate and multi-variate data (Section 5)

2 Theory andModels

Here we provide the machinery needed to create multiattribute models and apply them to our setting.

2.1 Preliminaries

Given a set of size d of correlated random variables Xi, ∀i ∈ D = {1, ... , d}, we define the joint cumulative distribution function (CDF) as F(⃗x), as the probability that Xi ≤ xi ∀i ∈ D where ⃗x ∈ Rd . We then define the joint probability density function (PDF) in the usual way as:

The marginal density function for variable Xi is then simply its distribution without regard to all the other variables:

Where 〈·; ·〉 represents vector composition. Here the integration ranges over all the dimensions but i, the marginal being specified. We will denote the corresponding CDF as Fi( ⃗X).

Complex multivariate distributions are notoriously hard to model and fit. Consequently, we here appeal to the powerful formalism of copulae. Copulae enable the construction of a model that captures the various marginal distributions of a empirical data separately from the "coupling" of these marginals. As we shall see, this separation makes both forming and fitting a complex model far more tractable.

The word copula was first introduced by Sklar [1959] in a theorem that described the necessary properties of functions that "join together" one dimensional marginal distributions to create multivariate distributions. Concretely, a copula function C is a multivariate distribution defined on the unit hypercube [0, 1]d with uniform marginals:

where Ui ∼ Uniform(0, 1), i ∈ D, along with the boundary conditions that:

Furthermore, according to Sklar's theorem [Sklar, 1959], if CDF F(⃗x) is a continuous multivariate distribution, a unique copula function C : [0, 1]d → [0, 1] exists such that

for continuous marginal distributions Fi(xi). Consequently, copula C links the marginal distributions together to create a multivariate distribution. In theory, such a copula function can be arbitrarily complicated. We seek a simple parametric copula that fits our data well. Next, we describe several such copulae from literature that we evaluate in Section 5.2

The Gaussian copula, which is based upon the multivariate Normal distribution, is one such choice for modeling dependence. The Gaussian copula is defined as:

where Φ is multivariate Normal CDF, Σ is a correlation matrix for the different variables, and the F1, F2, ... , Fd are arbitrary but well-defined marginal distributions.

Another important class of copula functions are the Archimedean copulae. Archimedean copulae have simple parametric form that allows for a variety of different dependence structures. The general form of an Archimedean copula is

where ψ is a decreasing univariate function known as generator of the copula, carefully chosen such that CA will obey the copula properties identified above. The Clayton, Gumbel and Frank copulae are each defined by a specific choice of Archimedean generator, as provided in Table 1, where λ is a free parameter that is available for fitting.

Joint dependence structure has been modeled using copulae in many fields, including reliability engineering, finance, and risk management. Yet, these powerful copula-based methods are unfamiliar to many researchers, and there is little literature comparing the efficacy of copulae in real settings.

2.2 Modeling the Marginal

Before building our full multivariate model, we must first ensure that we have accurately described the marginals in our domain, i.e. that we have appropriate univariate models for each of the demand attributes. Therefore, we next consider several distributions of demand for specific resource types including a new one not appearing elsewhere in the literature.

Process CPU usage (or runtime on a fixed processor type) is typically a key feature of computational loads. Early studies show that the general shape of the process runtime distribution is skewed [Walter and Wallace, 1967; Rosin, 1965], i.e., there are a lot of short jobs and a few long jobs. These papers use a simple exponential PDF to model process lifetime: f(x) = 1 θ ex/θ , where θ is a scale parameter. Later Leland and Ott [1986] and Harchol-Balter and Downey [1997] showed that process life distributions are actually heavytailed and better modeled by a Pareto distribution, with PDF function f(x) = ξxξm/xξ+1 where xm is location parameter and ξ is a shape parameter. There are also many studies of memory usage and its linear correlation with respect to CPU usage. For example, Parsons and Sevcik [1996], as well as Li et al. [2004] studied workload characteristics, including job size and memory usage. Chiang and Vernon [2001] found that most jobs have a very small memory usage per processor (under 32M) while 5% of jobs have a memory usage greater than 1 GB per processor, and a significant correlation between requested memory and number of processors. Work of this style describes memory usage statically; there has not been to our knowledge a distributional model for memory propsed. Moreover, none of these studies goes beyond simple linear correlation to model the full joint distribution over demand attributes as we do.

In our work, we have found empirically that the distribution of CPU usage (and memory) are not always decreasing, and a skewed bell-type structure is in fact possible. When confronted with such a bell shape, the first distributions one typically considers are the Normal and Log-Normal distributions. However, these are exponential family distributions which will not fit heavy-tailed data well. Because we believe our data may be heavy tailed, we additionally consider a somewhat more esoteric distribution, the Burr XII [Burr, 1942], also known as the Singh-Maddala distribution. The Burr XII distribution can fit a wide range of empirical data, as different values of its parameters cover a broad class of skewness and kurtosis. Hence, it is used in various fields such as finance, hydrology, and reliability to model a variety of data types. We choose it as a target because it can take on both a decreasing and a bell shape; it also generalizes the Gamma, Generalized Pareto and the Log-Logistic distributions which would otherwise be candidates in their own right [Burr, 1968; Tadikamalla, 1980].

The three parameter Burr XII PDF is given by

and CDF

where ξ > 0 and k > 0 are shape parameters and β > 0 is the scale parameter. Because of its flexibility and appropriate properties, we will consider the Burr XII distribution both for CPU and memory marginals.

2.3 Modeling the Joint Distribution

Having modeled a single attribute like CPU runtime or memory usage by an appropriate marginal distribution we now turn to handling the joint distribution of multiple attributes. At present, we focus on a joint distribution of CPU and memory because this is what is in our data, although a generalization to further dimensions is immediate.

We first model the marginals as in the previous section. Then, we perform a probability integral transform on the data along each marginal. This converts a random variable A to a new random variable B, which is uniformly distributed on [0, 1] by applying the CDF of A, FA to itself: B = FA(A). Concretely, this means we apply each of our fitted marginal distributions Fi ∼ Burr XII to the corresponding dimension of our data. This procedure produces a distribution on the unit square with uniform margins but some more complex joint structure. We then seek a parametric model for this jointuniform distribution: i.e., a copula as defined in Section 2.1.

By compositing our chosen copula with the marginal models, we then obtain a model of the full joint distribution that has the appropriate marginal distributions exactly, and where the joint dependence of the model is represented solely in the copula. This enables us to both define and fit a complete joint distribution in a principled manner.

3 The Data

Large open data sets for memory and CPU load are hard to come by, with the notable exception of the public Google workload traces we employ [Reiss et al., 2011]. The traces are compromised of 11, 000 machines' logs from 29 days in May 2011. We restrict the set to the 54.7% of jobs that completed successfully, resulting in about 760, 000 processes. Reported CPU time is in microseconds, and memory usage is scaled between 0 and 107 . We interpret a process' memory usage as the maximum observed over its runtime.

Before analyzing the data set some preparation is required. First, processes with 0 usage are omitted, which is a further 5% of the data. Further, while the bulk of the data comes from a reasonable spread of job types, the data includes a set of extreme outliers that are different in kind to the balance of the data, and thus the raw empirical distribution is multimodal in the tail. We believe this is a limitation of the data size (large though it is), and that if the number of job types included in the data were larger, these modes in the tail would disappear. To avoid this confounding our results, we restrict ourselves to the bulk of the jobs outside of these extreme outliers, or specifically to CPU usage below 5 ∗ 106 and memory usage below 7.5 ∗ 104 , which is 82% of the raw data. Next, to interpolate between the job types (modes) in the data, we perform a bivariate kernel-density smoothing with a Normal kernel and a bandwidth of 0.4 and 0.3 for memory and CPU respectively.

4 Bootstrap Goodness-of-Fit Testing

Algorithm 1 The Bootstrap Kolmogorov-Smirnov goodnessof-fit test for a hypothesized family of distributions against observed data.

Given an observed data set and a family of distributions that we hypothesize the data to be drawn from, we would like to know whether our hypothesis is a reasonable one, given the data. A standard approach to answering this kind of question is to use a goodness-of-fit test, which quantifies the plausibility of the hypothesis. In standard statistical hypothesis testing, one calculates a test statistic p, that captures the probability of the data being at least as "extreme" as that observed, given that the hypothesis is true. When this probability is less than a specified significance level, α, we reject the hypothesis. In typical usage, the hypothesis is a so called null (e.g. the data is explained by randomness), and so a rejection of this hypothesis in favor of an alternative and mutually exclusive one held by the researcher is the desired outcome, in that it provides evidence in favor of the alternative hypothesis. Thus, concretely, for disproving a null hypothesis you want p < α. Here though, in goodness-of-fit testing, we wish to find evidence that our hypothesis is true (i.e., the data is consistent with our fitted distribution). Consequently, a small pvalue (i.e., below α) is evidence against our proposed model, and we want p > α, ideally significantly so [D'Agostino, 1986].

To perform such a test, it remains to choose the notion of "extreme" that defines p. Many such test statistics have been proposed. In this work, we employ the commonly used Kolmogorov-Smirnov (KS) statistic: DKS(O, F) = supx{FO(x) - F (x)} where FO(x) is the empirical CDF of the observed data calculated as: FO(x) = 1n ∑ni I(Xi ≤ x) where I is the indicator function which is 1 when the given condition is true and 0 otherwise. DKS measures the maximum point-wise difference in the CDF of the hypothesized distribution and the observed data.

When the hypothesis is that the data is drawn from a concrete distribution, not a parametrized family, then a simple goodness-of-fit test is available: One finds the quantile of √nDKS in the well-known Kolmogorov distribution, and this quantile is taken as the value of p.

However, when testing the goodness-of-fit of a distribution fit from a parametrized family of distributions, the Kolmogorov distribution is no longer the appropriate target of comparison for DKS. Consequently, we employee a method known as parametric bootstrapping [Clauset et al., 2009]. In parametric bootstrapping, the p value is based not only on the test statistic (e.g., in our case DKS) computed over the observed data and the fitted hypothesized distribution, but also over synthetic data drawn from this fitted distribution and refits of the hypothesis model to these synthetic data.

Specially, we first fit our model to the observed data by estimating the parameters using Maximum Likelihood Estimation (MLE). MLE chooses the parameters which maximize the likelihood of observing the data under the model. Then we calculate the test statistic over the observed data and our fitted model. Next, we draw b synthetic data sets from the fitted model using the parameters we have just estimated. For each, we run MLE again to estimate the parameters of a new model from the hypothesized family of distributions, and then calculate the test statistic between the synthetic sample and its specially estimated private model. The p-value is then defined to be the fraction of the synthetic data whose test value is larger than the test value obtained from observed data. Algorithm 1 shows the steps of this parametric bootstrapping algorithm formally. The effect is a Kolmogorov-Smirnov (KS) test that is unbiased when used on hypotheses involving fitted families of distributions. The larger the number of synthetic data sets generated, b, the more reliable the p-value obtained. In our analysis we let b = 1000.

Our data was processed in Matlab, and analyzed in R. We used the copula [Yan, 2007] and CDVine [Brechmann and Schepsmeier, 2013] packages to perform our copula-based fits and as the basis for implementing Algorithm 1.

5 Results

Using this framework, we next fit models of CPU and memory usage both separately, and subsequently together as a joint model of demand.

5.1 Fitting the Marginals

In both the CPU and memory usage data, there are generally sufficiently large numbers of big jobs that we expect heavytailed distributions to perform well. We fit our data using Burr XII, Log-Normal, Generalized Pareto and Log-Logistic distributions by MLE. Table 2 lists the parameters as computed. 1

To evaluate the model fit, we begin by visually observing the correspondence between the empirical CDF and the CDF of each distribution in Table 2. Figures 1 and 2 illustrate this for CPU and memory usage respectively. We observe that the truncated Burr XII distribution is the closest to the empirical distribution for both CPU runtime and memory usage.

To back up the intuition we obtain from the graphs, we formally evaluate the goodness-of-fit of each distribution using the bootstrap KS method from Section 4, evaluating the CPU and memory data separately. Accordingly, we sample 5000 data points and run our KS test procedure using the distributions and parameters from Table 2 as our hypothesis. The resulting p-value is also shown in Table 2. From the table it is clear that Log-Normal, Log-Logistic and Pareto distribution perform poorly, with the differences between their fitted distributions and the data being strongly statistically significant. The calculated p-value of the Burr XII is very large, indicating no evidence to distinguish it from the observed data for CPU and memory. We therefore adopt it as our marginal model of both resources (with different parameters for each).

5.2 Fitting the Copula

Having created a model for each of the marginal distributions, we next turn to modeling the copula needed to relate them. To get a sense of the joint data (and to use as a baseline comparison for the fidelity of our eventual model) we show a scatter/ histogram plot of a 2500 point sample of the data in Figure 3.

Next, in Figure 4 we apply the probability integral transform described in Section 2.3 to the data in Figure 3. Specifically, we use the Burr XII marginal distributions for both CPU and memory usage with the parameters provided in Table 2. The result is the empirical copula, i.e. the "coupling" structure manifest in the data once the marginal distributions have been made uniform. We can see that our truncated Burr XII marginal distribution is doing an effective job at capturing most of the marginal structure, and thus the probability integral transform yields data that is nearly uniform in each dimension.

1 We fit truncated distributions, at 5 ∗ 106 for CPU and 7.5 ∗ 104 for memory.

Having applied our marginal models, we next seek a parametric copula that captures the structure observed in Figure 4. Accordingly, we fit five of the most commonly used families of copula including Gaussian, Frank, Gumbel, T, and Clayton. Table 3 shows the estimated parameters for each family. To get a sense of the copulae produced, we plot a random sample from the fitted Frank copula in Figure 5. We can see that this plot is a close match to the transformed empirical data shown in Figure 4. The Frank copula, as observed in Figures 5 has weak correlation in the tails compared to other copula [Trivedi and Zimmer, 2007]. However, there is still a tendency for there to be stronger correlation between CPU and Memory among either extremely large or extremely small jobs, which is reflected in our data set.

Next, to evaluate goodness-of-fit, we apply our bootstrap KS test algorithm, but now on the bivariate copulae instead of univariate marginals. The resulting p-values are listed in Table 3. From the table it is clear that the p-value for all the copulae except Frank and Clayton copula are small, indicating evidence to reject them. These large values for the Frank and Clayton indicate that we have no evidence to reject them as reasonable models of the coupling between the CPU and memory properties. Although both the Frank and the Clayton are reasonable models, the evidence for the Frank is stronger than the Clayton.

The last step of validating our model is to evaluate the complete model including both the fitted Burr XII marginals and the Frank copula together. Figure 6 shows a random sample of 2500 pints from this model, and appears nearly identical to the original data in Figure 3. Once again, we can formalize this by using a KS test to compare the observed data and data drawn from our fitted model. Table 3 shows the resulting p-values, for each copula we considered. The large p-value shown for the Frank copula indicates that our aggregate Burr XII-copula model can not be statistically differentiated from the data, and thus forms a valid model of joint demand.

Figure 6: A draw of 2500 points from our full end-to-end multivariate model, using the Frank copula parametrized in Table 3 and illustrated in Figure 5. Clearly the model closely matches the original data from Figure 3.

6 Conclusion

In this paper, we have analyzed the CPU and memory usage of jobs separately and as a joint distribution using real data from Google data-centers. Our marginal models of individual process attributes (like CPU and memory) have better fidelity than those used previously. We obtained this improvement by using a more flexible distribution, the Burr XII, than previously considered. We then leveraged this accuracy in modeling marginal distributions to create the first effective multi-attribute model of data-center demand. Concretely, we show that an Archimedean copula, specifically the Frank, can be used to great effect in modeling the "coupling" structure between the attributes of computational processes in distribution. Then, by putting both our Burr XII marginal models and our Frank copula model together, we obtain a highly accurate model of the full joint distribution across attributes.

Understanding multivariate resource usage distributions is critical to creating better resource allocation and scheduling algorithms for large computing centers, which in turn offer the promise of significant energy savings and CO2 reduction.