{'fraction_non_alphanumeric': 0.03402773714085084, 'fraction_numerical': 0.011323073322020013, 'mean_word_length': 4.734731543624161, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Forecasting of commercial sales with large scale Gaussian Processes

This paper argues that there has not been enough discussion in the field of applications of Gaussian Process for the fast moving consumer goods industry. Yet, this technique can be important as it e.g., can provide automatic feature relevance determination and the posterior mean can unlock insights on the data. Significant challenges are the large size and high dimensionality of commercial data at a point of sale. The study reviews approaches in the Gaussian Processes modeling for large data sets, evaluates their performance on commercial sales and shows value of this type of models as a decision-making tool for management.


I. INTRODUCTION
This study seeks to contribute to a better understanding in the industry to the field of forecasting with Gaussian Processes (GPs) with a focus on large data sets of commercial sales data for consumer goods in the fast moving consumer goods (FMCG) industry. Whereas there has been a wide interest in academia to explore Gaussian Processes for small data sets and on developing methods to fit large data into Gaussian Processes, there has not been enough diffusion for the commercial practitioner. Yet, GPs have found success on various applications. They are in use on a wide range of fields and use cases ranging from surrogate modeling, experiment design, mining and geo-spatial data to battery health [3], [71], [87], [11], [33], [32], [8], [10], [5], [6], [9]. On techniques for demand forecasting in the FMCG sector, there has been little academic research and not enough efforts to expose practitioners to them. Other forecasting techniques for consumer demand are widely used in industry and academia. For example, approaches using time series analysis, support vector machines (SVM), neural networks and splines [58], [31]. The industry has not been completely oblivious to employing GPs to forecast demand, for example, in electricity [7], [73], windmills [16], e-commerce [52], water [84] and tourism [85], [19]. Similarly, the consumer staples industry has made an effort to apply principled approaches to demand forecasting. One case was for tobacco products using seasonal time series decomposition and support vector machines to predict aggregations either of total industry demand [78] or of an individual producer [31]. GPs is a popular non-parametric model for regression and classification tasks leveraging the power of the underlying Bayesian inference framework [69], [12]. One shortcoming is that the computational performance of GPs deteriorates fast for data sets with more than thousands of observations [55] such as those found in commercial settings [74]. The objective of this research is to provide an overview for the practitioner of the theory and implementations of GPs in a big data setting using large commercial data sets as examples. The study argues that GP is a viable tool for sales forecasting and demand prediction in the consumer staples industry at a point of sale (POS). Recent advances with a focus on big data make them robust methods for multidimensional data [55], [88], [47], [37], [23]. Similarly, their flexibility as predictors qualify them for the particular characteristics of demand prediction for consumer products. This study evaluates and describes the characteristics of GPs for big data applied to commercial data sets of the FMCG industry and assesses related software libraries. For this purpose, it poses the following questions: (1) What is the state of the art in academic research of Gaussian Processes for big data? (2) Which implementations of Gaussian Processes are currently available? (3) How do GPs perform in practical tasks in the FMCG industry?

A. Significance & Innovation
Demand planners continue leveraging traditional forecasting techniques [14]. From an academic perspective, studies on sales forecasting either do not cover FMCG or use regression techniques considered to be special cases of GPs [12]. While valuable, these techniques do not benefit from the advantages of a probabilistic approach such as GPs. Moreover, there has been a limited number of academic studies in the area of GPs combining the disciplines of management sciences and computer sciences [29], [30]. An additional difficulty to assess the suitability of GPs for big commercial data is the popularity in the literature of small and synthetic data sets for benchmarks. This study innovates by combining a multidisciplinary approach on how GPs are evaluated from a theoretical and practical perspective and by bringing closer the GP theory and application to the industrial practitioner. Apart from a methodological innovation, this research seeks to contribute to a better understanding of Gaussian Process regression (GPR) for large commercial data in the FMCG industry.

II. RELATED WORK
There is vast literature on GPs, for example Rasmussen & Williams [12]. Summarized, a Gaussian process is a stochastic process h(x) such that any sub-collection of random variables h(x 1 ), h(x 2 ), . . . , h(x n ) for any finite subset of elements x 1 , . . . , x n ∈ X ⊂ R d follows a multivariate Gaussian distribution. The key characteristic of a GP is that it is fully specified by a mean function m(x) and a covariance function k(x, x ). Thus, The GP generalizes the Gaussian distribution going from vectors to functions. Further, it allows for efficient computation since it only needs values at a finite number of points. Moreover, the smoothness of a GP is controlled by the covariance function k(x, x ). For Gaussian Process regression (GPR), the prior is assumed to be a GP, see definition (1), with a zero mean. Suppose there are n points x i , i = 1, 2, . . . , n and corresponding labelsȳ = (y 1 , . . . , y n ). For a new point x * , y * = f (x * ) should be predicted. By assuming that the regression function f is a Gaussian process, it is obtained that the joint distribution of the training outputsȳ and the test outputs y * , both vectors, at any finite number of points is a Gaussian ȳ y * ∼ N 0, This is called the joint prior, where the mean equals zero, K is the covariance matrix for the training points x i , K * is a covariance vector between the test x * and the training points x i , and K * * is the inherent measurement noise. To get the posterior distribution over functions, it is necessary to condition this joint Gaussian prior distribution on the observations as following

A. Advantages & Disadvantages
Some of the advantages attributed to GPs are according to Do [4]: they enable the possibility to quantify uncertainty in predictions resulting not only from the intrinsic noise in the problem but also from the errors in the parameter estimation procedure. Further, many methods for model selection and hyperparameter selection in Bayesian frameworks can be applied. Moreover, GPR is non-parametric. Thus, it can model any arbitrary function based on the input points. In addition, they present a natural way to introduce kernels into the regression modeling framework. At the same time, their main disadvantages are the computation and storage costs, dominated by the inversion of K, as well as the difficulty to understand conceptually the theory behind. Moreover, they do not deal well with discontinuities such as those found in financial crises, phosphorylation, collisions, edges in images, etc. In addition, Lawrence claims that the popular squared exponential covariance (RBF) often used in the literature is too smooth for practical problems [4], [55].

B. Gaussian Processes for Big Data
The largest limitations of using GPs in practice are the storage O(n 2 ) and computational O(n 3 ) requirements for n data points. Nevertheless, the strengths of GP such as its flexibility, its conceptual simplicity and desirable properties make it highly attractive to address regression tasks. For this reason, there has been an active interest in the research community to extend GPs to large data sets. This study organizes GPs methods for big data around four themes: (1) Sparse approximations, (2) Low-rank approximations, (3) Local approximations, (4) Parallelization. Other authors such as Hoang [49] and Rasmussen & Williams [12] follow other classifications.

C. Sparse approximations
The general idea of a sparse approximation is to have a sample covariance matrix that is sparser than its original. The inversion of the sparse matrix is less computationally expensive than the inversion of a non-sparse matrix of the same size. This category has a very rich list of heterogeneous approaches. Furrer et al. [36] and Kaufman et al. [53] apply a covariance tapering technique. On the other hand, Gneiting [39] uses a compactly supported covariance function. Lindgren et al. [56] use the Gaussian Markov approximation of a GP. The work of Grigorievskiy et al. [45] for temporal GPs with a 1-dimensional input space uses the sparseness property of precision matrices in a Markovian process to scale computationally to O(b 3 n) where b is the size of the matrix block. On a similar direction, Gilboa et al. [38] also use a Markov process under the umbrella of the projection pursuit method for structured GPs. A GP can be considered structured, if its marginals contain exploitable structure enabling a reduction in computational complexity. Another family of approximations is based on approximate matrix-vector-multiplications (MVMs). Some of these methods have been reviewed by Quinonero-Candela et al. [67]. Local mixtures of GPs have been used by Urtasun & Darrell [82] for efficient modeling of human poses. Gal et al. [37] use a re-parametrization of variational inference. On a different direction, Davies & Ghahramani [22] propose focusing on the kernel and using matrix-free methods. They do not require the full Gram matrix K, only the ability to calculate Kv for any arbitrary v [79].

D. Low-rank approximations
Historically, most of the work on GPs for large datasets has been focused on this area according to Hoang et al. [48]. The low-rank approximate representation of the full-rank GP (FGP) is an alternative family for sparse GP regression. Here, the idea is to build a low rank approximation of the covariance matrix based around so-called 'inducing variables'. Examples of this are the works of Csato & Opper [20], Seeger et al. [75], Quinonero-Candela & Rasmussen [66] and Titsias [79]. Lowrank approximations introduce latent variables and assume a certain independence conditioned on the latent variables leading to a computational complexity of O(nm 2 ) and storage demands of O(nm) with m as the parameter governing the number of selected inducing variables [75], [77]. For example Hensman et al. [47] do this with variational inference. According to Low et al. [57], low-rank methods are well-suited for modeling slowly-varying functions that are largely correlated, and where it is possible to use all the data for predictions. The downside is that they require a relatively high rank to capture small correlations with high fidelity; this makes them lose attractiveness. Authors such as Low et al. [57] try to overcome this by leveraging the dual computation of complementing a low-rank approximate representation of the full-rank GP based on a support set of inputs with a Markov approximation of the resulting residual process. Hoang et al. [48] see here two large subgroups. On one side the distributed approach, where it is sought to reduce the training time with all the data by a factor close to the number of machines. Some examples here are the works of Chen et al. [15], Hoang et al. [48], [51] and Low et al. [50]. The second subgroup is the stochastic implementation. The main idea behind is twofold. First, to train with a small, randomly sampled subset of data in constant time per iteration of stochastic gradient ascent update. Second, it is to achieve asymptotic convergence to their predictive distributions. One example is Hensman et al. [47]. Das et al. [21] follow a different approach and opt for an empirical method of first using a bootstrapped data set to develop a GP and then to bag the models to produce the regression estimate. Finally, Hoang et al. [48] argue that there is a less wellexplored class exploiting sparsity in the spectral representation of a GP kernel.

E. Local approximations
In this approach the general idea is to partition the data set into separate groups. Exponents of this concept are Snelson & Ghahramani [76] and Urtasun & Darrell. Low et al. [57] mention that this family comprises ideas around localized regression, covariance tapering methods and compactly supported covariance functions. Their strength is that they address the shortcomings of low-rank approximations; they can model rapidly-varying functions with small correlations. As they only use local data for predictions, they end up performing poorly in input regions with little data. One class of local approximations partitions the input domain into a set of local regions and assume an independent GP regression model within each region. The resulting sample covariance matrix is a block diagonal matrix of local sample covariance matrices. Park & Apley [62] highlight that inverting the block diagonal matrix is much cheaper computationally. Such local approximation approaches have many advantages. By their local nature, they adapt better to local and non-stationary data features and independent local approximation models can be computed in parallel to reduce total computation time. Their major weakness is that two local models for two neighboring local regions produce different predictions at the boundary between the regions, resulting in boundary discontinuity for the regression predictive function. This boundary discontinuity implies greater degradation in prediction accuracy, particularly around the boundaries of the local regions [63]. The discontinuity issue has been addressed in different ways. The most popular approach is to smooth out some of the discontinuity by using some weighted average across the local models or across multiple sets of local models via a Dirichlet mixture (Rasmussen & Ghahramani [70]), or via a treed mixture (Gramacy & Lee [44]), or with a Bayesian model averaging (Tresp [81], Chen & Ren [18] and Deisenroth & Ng [26]), or with locally weighted projections (Nguyen-Tuong et al. [61]). Other related approaches use an additive covariance function consisting of a global covariance and a local covariance. Examples are the works of Snelson & Ghahramani [76] and Vanhatalo & Vehtari [83]. Another alternative is to construct a local model for each testing point (Gramacy & Apley [42]), or to use a local partition but constrain the local models for continuity. This last proposal is reflected in the works of Park & Huang [63] and Park et al. [64].

F. Parallelization
More than a technique to deal with data, parallelization is rather used jointly with other techniques described above. Sparse GPs can for example make use of existing paradigms in distributed computing such as Apache Spark [86] or Map-Reduce [25]. In the work of Gal et al. [37], they re-parametrize variational inference for sparse GPR, see section II-D, to reformulate the evidence lower bound in a Map-Reduce setting. Low et al. [57] propose a low-rank-cum-Markov approximation (LMA) of the full GP. This work falls also into the family of methods discussed in section II-E. Similarly, Zhang & Williamson [88] propose a parallel local method under the manteau of 'embarrassingly parallel' algorithms, where the global communication occurs only after the local computation is complete. On the same vein, Deisenroth & Ng [26] exploit the fact that a single GP can be split into q independent problems, whose parameters can be inferred independently (and parallel) of each other in a weighted product-of-experts model [27], [28]. An example of a sparse approximation, see section II-C, with a parallel component is in the work of Grigorievskiy et al. [45]. Another area of research lies on the parallelization of GPs through distributed computing or a combination of CPU and GPU as shown by Gramacy et al. [43]. They show that the combined effects of approximation and massive parallelization can be applied to GPs. However, this requires careful work.

III. ADVICE FOR THE PRACTITIONER
Inspired by the work of Davies [23], this study suggests following best practices for the practitioner. In kernel-related methods, the first and most important step is to select an appropriate kernel. It should be remembered that the kernel encodes all the domain knowledge needed. Many properties can be embedded into the kernel such as smoothness, periodicity, linearity, dependence between dimensions, etc. The literature on kernels and kernel choice is vast. The reader is advised to choose a covariance function that reflects the known properties of the underlying process as much as possible. On a second step, it is necessary to choose the ideal approximation. There is not a 'one-size fits all' approximation akin to the choice of kernel. The third step is to determine the required level of accuracy for the posterior. This is closely related to the choice of kernel and of approximation. Often, problems do not require an accuracy that justifies a computational cost of O(n 3 ). The practitioner should evaluate the problem at hand and determine the desired trade-off between computational speed and accuracy. Frequently, hyperparameters need to be optimized. Thus, the fourth step is choosing an ideal optimization. A popular choice for hyperparameter optimization is the log-likelihood. Finally, probably the most important step is data preparation. Often, not all data is necessary for a problem. It is desirable to work with a subset of the original data set.

IV. IMPLEMENTATIONS OF GAUSSIAN PROCESSES
To strengthen the notion that GPs are both popular and used in practice, table I 1 portrays a comparison between GPs and other forecasting techniques on diverse Internet websites. For example GP regression is frequently listed as part of a repository description in Github 2 . Similar results were obtained from comparing the level of popularity since 2004 using Google Trends 3 and from counting the number of threads under the respective tag at CrossValidated (CV) 4 as well as at Reddit 5 respectively. The results show that GPs, although still obscure, have a stronger presence in public forums than 'classical' regression methods such as Splines and Kernel ridge regression (KRR). Similarly, it is of relevance for the practitioner to have a library of choice for prediction with GPs. This study compares popular implementations of GP regression under a criteria inspired by Golge [40] and De et al. [24] considering soft factors such as an active community and hard ones such as computational performance.

V. GAUSSIAN PROCESS LIBRARIES
To assess different implementations of GPR, this research assigns up to three points, (1) acceptable, (2) good and (3) very good, to each of the following aspects. h) Object Oriented: Libraries with object-oriented APIs are desired for deployment in production environment.
i) Test coverage: Often, organizations set guidelines on a minimum test coverage to be fulfilled. For the practitioner seeking to introduce GPs, test coverage is essential. For the purpose of this analysis, a set of GP regression libraries was selected in table II. Sci-Kit Learn [65], GPy [41], GPflow [24], Edward [80], George [35], [2], pyGP [60], pyMC3 [72]. Implementations in other programming languages were also assessed. In R, GPfit [59] and gptk 6 are popular. However, they have not been updated in the last years. Julia, another language used by statisticians, does not appear to have an established GP library. This study found juliaGP 7 , MLKernels 8 and two libraries named GaussianProcesses.jl 910 . As of this writing, none of them is either mature or active. Stan [13], another popular programming language among Bayesian statisticians has a strong community, extensive documentation and an API for Python, PyStan. However, to narrow down the scope of this study, only primarily Python libraries are considered. This is not the first study evaluating GP software empirically. Erickson et al. [34] also compared multiple GP implementations. Their conclusion was that Scikit-learn and GPy are among the best options. However, they warn that two different implementations can give significantly different results on the same data set and therefore advice to select a library on a by-case basis. The libraries listed below were chosen for further evaluation with table II summarizing the assessment.
j) Scikit-learn: It includes implementations of almost all major machine learning methods. The community is large and active. Due to its popularity, there is extensive documentation and examples. The library is stable, with clear release cycles and support for legacy versions. On the other hand, it is not optimized for GPU. Similarly, it is less-flexible than other libraries evaluated, object-oriented and has extensive test coverage.
k) GPy: The most popular library in the GP community among libraries exclusively dedicated to GP modeling. The library has a strong community and multiple examples. The documentation is good. Additionally, it is also object-oriented This library can be seen as a compromise between industrial requirements and academic needs. Downsides are limited GPU functionality and insufficient test coverage. l) GPflow: The focus is primarily on experimentation without sacrificing test coverage and speed. It uses full GPU acceleration and variational inference as the default approximation method. The downside is that the library is not as wellknown as others. The community is small but active. Similarly, although the standard documentation is comprehensive, the amount of available examples is rather limited. m) Edward: An alternative to Stan in Python. The focus is on Bayesian probabilistic modeling and on academic experimentation. Akin to GPflow, the standard documentation is optimal with limited examples outside the official material. It has a small and active community. Given that it leverages TensorFlow [1], it can speed up computations with GPU.
n) George: A library with a focus on fast GP regression. Compared to the previous libraries, it is less popular. This is reflected on the sparse official documentation and on the lack of examples. In addition, development has stalled and the community support seems limited. Moreover, it does not offer hardware acceleration and the test coverage is basic. o) pyGP: It has not been updated in many years and the community is dormant. It lacks examples and the documentation is basic. p) pyMC3: Akin to Stan and Edward, pyMC3 is a probabilistic programming library. It is well-known in both academia and industry. The community is large, the documentation comprehensive and many examples are easily found. Given its level of industrial adoption, the library is stable and has well-known development cycles. PyMC3 is an interesting option for the industrial practitioner interested in Bayesian inference on a production-ready environment. Based on the scores from this assessment, Scikit-learn, pyMC3 and GPflow are the stronger choices for GPs in Python.

A. Bimbo data set
There have been previous analyses on this data set, for example by Kosar et al. [54]. Important to highlight, the test set does not contain numerical variables. Thus, demand has to be predicted using categorical features. Additional files with attributes related to POS and their location were also made available. There is a relation between sales and returns. Locations with low levels of sales have high level of returns. In addition, evaluating the data showed irregularities in its collection. After pre-processing, the resulting data set contained 147 features. The five most relevant based on their

B. Fortune 500 company data set
Data generated at POS for 28 weeks in 2016 was provided containing attributes such as sales volume per SKU, type of retail outlet, location, stock levels, trade marketing activities & more. Figure 1 exemplifies the stability of the industry by decomposing and plotting the time series of the log demand for 198 days. It is possible to observe a seasonality component and a stable demand. Based on discussions with industry experts, it was decided to discard non-relevant attributes. The reduced data set contained following attributes: (1) 180000 points of sale, (2) 2085 cities, (3) 6 POS types, (4) 39 trade categories, (5) 35 products (SKU), (6) 282 depots. Out of 112 features, the top 5 features according to XGBoost were (1) Revenue by week, SKU and POS, (2) Total volume difference by POS between two consecutive weeks, (3) Total weekly volume by POS, (4) Share of volume by POS at national level, (5) Total revenue difference at depot level between two consecutive weeks. Whereas figure 1 hints at the possibility of using ARIMA or related methods, it would be difficult to obtain significant good results at POS level. The observations by POS were collected on weekly intervals at best and not all SKUs were tracked on each visit. Thus, each POS has at most 28 (7 for Bimbo) points to predict demand for each of the 35 SKUs. This was corroborated during the pre-processing phase. A POS is visited with varying frequency, every 3 to 21 days. On each visit, there could have been notable changes in inventory levels, stocked products and sales volume.

VII. EXPERIMENTS
XGBoost was used on as a baseline for comparisons and to obtain the most relevant features. Additionally, for both data sets a correlation matrix was generated. The observations were fitted onto variations of GPs for big data sets such as Parametric GP (PGP) [68] and Variational Fourier features for Gaussian processes (VFF) [46]. These two methods were chosen for multiple reasons. First, they represent distinct approaches to GP for big data (sparse vs low rank). Second, they are recent approaches in the literature, published in 2017 and 2016 respectively. Third, they can be implemented using GPflow benefiting from GPU acceleration. Fourth, they have shown good results with large data sets. a) Hardware: The models were trained on a server with the operating system Ubuntu Yakkety, a processor Intel Xeon E5-1650 v3 Hexa-Core Haswell, 256 GB DDR4 ECC RAM, a SSD hard drive and a graphic card NVidia GeForce GTX 1080.
For VFF, an additive Matern− 3 2 covariance function was chosen, whereas the squared exponential kernel was selected for PGP. This decision was made due to their flexibility and recommendations on the respective literature. Additionally, an automated relevance determination (ARD) kernel was used to identify the best features according to the GP.

VIII. EVALUATION & DISCUSSION
For the Bimbo data set, in table IV, XGBoost performed better achieving a RMSLE among the top 1% in the Kaggle competition. On the other hand, PGP was among the top 65% and VFF was placed among the top 20% best performers. This can be observed visually in figure 2, where the Y-axis represents the ground truth and the X-axis the prediction. Ideally, the predictive points should be as close as possible to the slope. It is important to highlight that VFF struggled evaluating the full data set. Thus, only the top 10 best features according to XGBoost and to the ARD kernel were used on separate instances. This highlighted the importance of feature selection. The lists of most relevant features according to the correlation matrix, XGBoost and the ARD kernel were very different and none of them can be considered the best one. For the Bimbo data set, the features from XGBoost showed the best results. On a similar vein, once XGBoost was put under the same constraints as VFF and fed only with the top 10 features, performance dropped significantly and fared worse than PGP and VFF. One valuable characteristic of GPs is the possibility to evaluate the posterior mean. Decision-makers are often interested in identifying factors that if optimized lead to a certain outcome. For this, the posterior mean can be used. Figure 3 depicts the log demand against the Top 5 features according to XGBoost. It can be noticed that the selected features do not have a large impact on log demand. Yet, it helps understand that an increase in returns by POS type impacts log demand, first by decreasing it and later by increasing it. Figure 5 draws similar conclusions depicting average returns by SKU and average returns by SKU and POS. It can be interpreted that returns are a sign of activity. For example, supermarkets have significant sales volume but also large amount of returns. However, once the returns become too large, it signals decreased activity at the POS and thus lower sales. Another insight drawn from evaluating the posterior mean in figure 3 is the relation between types of POS (e.g., supermarket, cafeteria, etc) and demand. For example, due to their characteristics, some types of POS generate more demand than others. For the Fortune 500 company data set, in table V, XGBoost achieved a RMSLE of 0.036. The best PGP has a RMSLE of 0.082 obtained in only 15 minutes versus more than 90 minutes for XGBoost. In comparison, VFF showed a RMSLE of 0.44. Figure 4 plots the ground truth against the predictions for the best models. It is possible to appreciate that PGP has a very low error, whereas VFF struggled. The differences between XGBoost and GPs can be attributed to various reasons. On one side, the experiments with VFF used only top features, whereas XGBoost was trained on the full data set. Once XGBoost was trained with the top 10 features only, performance dropped. Similarly, observations with null values had to be discarded to be fitted the GP models. XGBoost, on the other hand, allowed for the presence of null values. Thus, the algorithm had both more and more diverse data at hand. VFF seemed to work better with mediumsized subsets in the range of the hundreds of thousands of observations. Nevertheless, for both GPs, it was evident that the choice of covariance functions has a larger impact in the performance over other considerations. This was corroborated when the number of features for VFF was increased from 20 to 50; the results were similar but at a significantly higher time cost.
a) Other regression techniques: Other regression methods such as multivariate adaptive regression splines, e.g., Py-earth 14 , KRR and SVM regression (SVMR) could had been included. However, the main limitations were the lack of production-ready implementations or absence of native libraries. On the other hand, implementations of SVMR and KRR are widely common in Python. However, they proved unsuited for the type of data sets used by this study achieving long training periods and discouraging results. 14 https://github.com/scikit-learn-contrib/py-earth Fig. 3: Posterior mean of log demand for Top5 XGBoost features for Bimbo data set b) Other considerations: A company benefits from a reduction in costs and from added productivity. This can be achieved with the usage of GP methods for big data. Although notable for being traditionally slow, new approaches to GP turned out to be competitive and even faster than state of the art frameworks such as XGBoost, with the added advantage that GPs did not require a significant effort in parameter tuning. This is replaced with efforts in developing covariance functions. It can be argued that the covariance function is a reflection of the problem task. Thus, the analyst is investing time on understanding the business problem and honing domain knowledge. Similarly, the speed of methods such as PGP allow for quick hypothesis testing and so-called rapid prototyping. In summary, GPs can achieve very similar results to XGBoost under the same conditions, with the added benefit of providing insights to the business on the factors driving demand by analyzing the posterior mean and thus serving as a support tool for decision-makers. With two large data sets from the consumer goods sector, it was possible to show that GPs can be used as a tool for modeling and prediction for very large data sets. The results are close to state of the art methods such as XGBoost. However, recent frameworks such as PGP proved to be significantly faster and scaled better than XGBoost. An additional benefit is the possibility to analyze the posterior mean conditioned on different features. This points out an additional use case for GPs as a decision support tool for management in the FMCG sector. In the case of Bimbo, the business can focus its resources on relevant POS with an over-proportional level of returns. This study innovated by providing an overview with depth and breadth on Gaussian Process for big data. It is one of the first academic works on GPs for FMCG with an additional emphasis on the practitioner. It surveyed more than forty works on GP methods for big data. Similarly, it assessed the popularity of Gaussian process among practitioners and benchmarked different GP frameworks. Nevertheless, there are still significant areas for innovation. From a theoretical standpoint, a further area of research is to evaluate the development of new kernels considering factors such as seasonality, product life cycle and out of stock events. Similarly, it is necessary to do further work on optimal feature selection for large and high-dimensional data sets. From a practical perspective, there is significant work to assess, which GPR methods for big data are better suited for the FMCG industry, and to evaluate the suitability of using the posterior mean as a decision support tool for management. From an implementation point of view, it is interesting to have a unified software framework covering all major GPR methodologies. In conclusion, this work provided an introduction and laid the ground work for future work in GPs for FMCG. Given the importance of this industry sector in the economy and the benefits of GPs, it is to be expected that there will be further academic work in this area.