Hyperparameter optimization for recommender systems through Bayesian optimization

Recommender systems represent one of the most successful applications of machine learning in B2C online services, to help the users in their choices in many web services. Recommender system aims to predict the user preferences from a huge amount of data, basically the past behaviour of the user, using an efficient prediction algorithm. One of the most used is the matrix-factorization algorithm. Like many machine learning algorithms, its effectiveness goes through the tuning of its hyper-parameters, and the associated optimization problem also called hyper-parameter optimization. This represents a noisy time-consuming black-box optimization problem. The related objective function maps any possible hyper-parameter configuration to a numeric score quantifying the algorithm performance. In this work, we show how Bayesian optimization can help the tuning of three hyper-parameters: the number of latent factors, the regularization parameter, and the learning rate. Numerical results are obtained on a benchmark problem and show that Bayesian optimization obtains a better result than the default setting of the hyper-parameters and the random search.


Fig. 1
Two different types of ML algorithm for RS: a Content-based approach, and b Collaborative filtering approach. Source: https://towardsdatascience.com/brief-on-recommender-systems-b86a1068a4dd in selecting items from a large set of choices. Example of applications can be found in many fields, among which movies (Koren et al. 2009), music (Lee et al. 2010), books (Crespo et al. 2011), e-commerce (McNally et al. 2011) and active stock selection (De Rossi et al. 2019). The idea behind an RS is that providing personalized suggestions significantly increasing the likelihood of a customer making a purchase compared to un-personalized ones. Personalized recommendations have huge importance where the number of possible items is large such as in e-commerce related to art (books, movies, music), fashion, food, etc. Some of the major participants in e-commerce (Amazon), movie streaming (Netflix), and music streaming (Spotify) successfully apply recommender systems to deliver automatically generated personalized recommendations to their customers.
Machine learning (ML) algorithms in Recommender systems are typically classified into two categories (Aggarwal 2016) (see Fig. 1): • Content-based approaches profile users and items by identifying their characteristic features, such as demographic data for user profiling, and product information/descriptions for item profiling; • Collaborative filtering approaches (CF) identify relationships between users and items and make associations using the past user activities information to predict user preferences on new items.
The first approach is laborious because it is necessary to collect information about users/items, and it often tricky because the users must share their personal data for the creation of a database for profiling. The CF approach requires few data, basically a list of tuples containing the user ID, the item ID, and the rating done by the user to that item. Moreover, the CF algorithms are more flexible in that they can be applied to RS independently of the domain of application.
In this paper, the authors focus on CF, in which the basic data structure is the rating matrix, formed by any possible user-item combination. The problem is also called the matrix completion problem, because it is based on an incompletely specified matrix of values (the users' past rated items), and the aim is to predict the remaining values (the predicted user' rates on new items) using some learning algorithm. The main challenge in designing CF methods is that real-world databases are mostly sparse (many unknown entries), but the unknown ratings are predictable because the known ratings are often highly correlated across various users or items.
Two types of methods are commonly used in the CF framework (Cacheda et al. 2011): the memory-based methods and model-based methods. The memory-based methods, or neighborhood-based CF algorithms, were among the earliest collaborative filtering algorithms, in which the ratings of user-item combinations are predicted based on their neighborhoods (users similar to a target user or items similar to a target item). They are based on the fact that similar users display similar patterns of rating behavior (user-based) or similar items receive similar ratings (item-based). An example of neighbourhood method is the k-Nearest neighbours (Yeehuda 2010). These methods are simple to implement, and the resulting recommendations are often easy to explain. On the other hand, memory-based algorithms do not work very well with sparse rating matrices: they scale poorly with the number of dimensions, and their predictions are not accurate for user/item matrix with few ratings.
The model-based methods are an alternative approach that try to predict the ratings by characterizing both items and users using a certain number of parameters inferred from the rating patterns. More in details, they are based on the assumptions that the preferences of a user can be inferred from a small number of hidden or latent factors. The most successful realizations of latent factor models are based on matrix factorization (Salakhutdinov and Mnih 2008;Koren et al. 2009). These approaches are considered superior to classic nearest-neighbour techniques for producing recommendations and learn the latent factors within an optimization framework. This corresponds to a low-rank approximation of the rating matrix, with the assumption of correlations between rows (or columns) to guarantee the dimensionality reduction of the matrix itself.
The RS becomes a minimization problem, in which we want to determine two lowrank matrices whose product is as close as possible to the rating matrix. The associated error function depends on the number of latent factors. Moreover, a regularization term is added to reduce both the non-linearity effect and the overfitting problem. Finally, if we solve this problem by means of stochastic gradient descent (Bottou 2010), we must set also the learning rate related to the descent direction.
These three hyper-parameters (i.e., number of latent factors, regularization term and learning rate) must be tuned through an optimization procedure, usually called hyper-parameter optimization.
The objective function maps any possible hyper-parameter configuration to a numeric score quantifying the quality of the matrix factorization. This score is computed by dividing the known entries of the rating matrix in two sets: a training test and a test set. For each hyper-parameter configuration, the stochastic gradient method is applied on the training set and the score is obtained as validation on the test set. This validation schema represents a randomize time-consuming black-box optimization problem. Therefore, an efficient global optimization strategy is mandatory to obtain the best possible configuration using a small number of function evaluations.
An example of the application of BO for RS is, e.g., in Dewancker et al. (2016) and Cano (2019). In the first case, the authors use a Bayesian optimization web-service, called SigOpt, and show that BO outperforms random search in tuning three hyperparameters of the Alternating Least Squares algorithm (Udell et al. 2016) to estimate the entries of the rating matrix. In the second case, the authors show the advantage in tuning the two hyper-parameters related to the KNN method to find the top five items recommended. Finally, alternative use of BO for RS can be found in Vanchinathan et al. (2014), where the challenge of ranking recommendation lists based on click feedback by efficiently encoding similarities among users and among items is considered. In this case, the Gaussian Process is used to model the elements of the rating matrix directly.
In our paper, we want to use BO for the hyper-parameter optimization related to the parameters of the stochastic gradient descent to find the best possible configuration, in terms of the learning rate, number of latent factors, and regularization parameter. We describe the mathematical formalization of the problem and we show an example of application of BO on a benchmark dataset, the MovieLens-100k, to find the best possible configuration of the hyper-parameters.
The rest of the paper is organized as follows. Section 2 introduces the problem definition and Sect. 3 describes the Bayesian optimization algorithm. A benchmark application is presented in Sect. 4, and Sect. 5 we present the conclusions.

The problem definition
In the most general framework, a CF problem is based on the definition of two sets (Takács et al. 2009): • The set of users U {u 1 , u 2 , . . . , u M }, where M is the number of users; • The set of items I {i 1 , i 2 , . . . , i N }, where N is the number of items.
Each user expresses its judgement, or rating, r ∈ X , where typical rating values can be binary or integers from a given range. The set of all the ratings given by the users on the items can be represented as a partially specified matrix R ∈ R M×N , where its entries r ui express the possible ratings of user u for item i. Usually, each user rates only a small number of items, thus the matrix elements are known in a small number of positions (u, i) ∈ S, with |S| min{M, N }(see Fig. 2). Usually, the dataset S is divided in a training set S T r and a test set S T e with S T r ∩ S T e ∅ and S T r ∪ S T e S, and the aim of CF is to create a prediction of the elements of S T e using only the knowledge of S T r , minimizing some error functions: where r ui r ui (S T r ) denotes the prediction on r ui , and is obtained as a function of the training set S T r . The typical used error function is the root mean square error (RMSE),

The matrix factorization algorithm
MF algorithms are the most successful ones for RS, and represent the state-of-the-art since they have shown their superiority in terms of predictive performance and runtime in numerous publications (Koren et al. 2009;Takács et al. 2009). The idea behind MF techniques is to approximate the matrix R as the product of two matrices ( Fig. 3): where P is a M × K and Q is a K × N matrix. The first matrix is called the user-feature matrix, the second one is called the item-feature matrix, and K represent the number of latent factors in the given factorization. Typically, K min{M, N }, and both P and Q contain real numbers, even when R contains only integers.
The aim of CF approach based on matrix factorization is to minimize the error function on the training set S T r , as a function of the matrixes ( P, Q). The prediction r ui in Eq. (1) is expressed by:r ui K k 1 p uk q ki (4) where p uk and q ki denote the elements of P and Q, respectively. Therefore, the optimization problem becomes After the minimization of Eq. (5) is done, one can estimate the RMSE on the test set S T e : The total number of variables of the optimization problem of Eq. (5) depends on the number of latent factors K. Generally, increasing the number of latent factors K improves the quality of the solution of Eq. (5) but can cause an overfitting problem. Basically, with too much freedom (too many free parameters), the model starts fitting well the training data but not generalizing well the unseen test data. A common way to avoid overfitting is to apply a regularization term to Eq. (5), obtaining a new but similar optimization problem where λ ≥ 0 is the regularization factor. Besides, a simple split in the training set and test set does not guarantee robustness of the results. A more suitable procedure, widely adopted, is the cross-validation approach. More in detail, using the N-fold validation, the original dataset is divided in n subsets. Then, the model prediction is computed for n times, using the nth set as test set, and the other n − 1 folds as the training set.
The results of this procedure can be, for example, the mean of the n different values of the RMSE on the n folds. In Fig. 4 we report an example of pseudo-code for a numeric score function, using the n-fold validation.

The optimization problem
The optimization problem associated with (Eq. 7) represents a non-convex and multiextremal global optimization problem. Global optimization methods could be used to solve it such as evolutionary algorithms or multi-star method. However, finding the global minimum results very hard to solve since also the high number of variables. Then, it is usually content to find only a local minimum through a local gradient method. The state-of-the-art approaches are the Alternating Least Squares and the Stochastic Gradient Descent (SGD), that is the one we consider for this work.
First, one need to compute the partial derivative of the objective function with respect to the variable p uk and q ki : At this point, one could update the entire vector of decision variables as indicates the entire vector of the partial derivatives and η is the step size or learning rate. However, when the data size is very large, it is preferable to use the SGD, in which the update is stochastically approximated in terms of the error in a (randomly chosen from a uniform distribution) observed entry (i, j) as follows: By this way, one can cycle through the observed entries in R one at a time (in random order) and updates only the relevant set of 2 · K entries in the factor matrices rather than all (m · K + n · K ). Indeed, for each observed rating r i, j , the error e ui r ui − K k 1 p uk • q ki is used to update the k entries in row i of U and the k entries in the row j of Q. A typical initialization for the elements P and Q is done randomly according, e.g., to a normal distribution with mean μ 0 and small variance σ . The pseudo-code for the stochastic gradient descent method is illustrated in Fig. 5.

Tuning the hyper-parameters
The optimization algorithm presented has three hyper-parameters to set: the learning rate η, the regularization factor λ, and the number of latent factors K. Different combinations of these hyper-parameters can alter the performance of the optimization algorithm significantly. If we want to know which parameter combination yields the best results, we must perform a hyper-parameter optimization.
Formally, we define the hyperparameter optimization task as follows. Let A be the target algorithm with n number of parameters to be tuned. Each hyper-parameter θ i can be a value taken from an interval [a i , b i ] (continuous or integer) in a hyper-parameter configuration space [a 1 , b 1 ] × · · · × [a n , b n ]. Defining a performance function H : → R + that maps each possible configuration θ ∈ to a numeric score, the aim of the hyper-parameter optimization is to find the best configuration θ * that minimizes H (θ): The objective function H is characterized by the following properties: 1. The function is time-consuming in the sense that each evaluation takes a substantial amount of time (minutes for a small dataset, hours for a bigger one); 2. The function is "black-box", that means we do not know any properties about its structure like linearity or concavity, and we do not have information about derivatives; 3. The evaluation of the function is characterized by noise. This fact is because the evaluation of the performance function depends on a series of random factors related to the stochastic gradient descent algorithm (e.g., the starting point of the gradient-based optimization) and on how the dataset is divided in training and validation set.
All these facts make the optimization problem hard to solve and require a specific procedure to obtain the best configuration in a few evaluations. In Fig. 6 we show the pseudo-code of hyper-parameter optimization.

Bayesian optimization for hyperparameter optimization
BO is a sequential model-based approach to solve the problem (Eq. 10) and, in general, global optimization problems based on expensive-to-evaluate black-box functions. The search space can be a compact subset of R d , but the BO framework can be applied also to search spaces involving categorical or conditional inputs. Moreover, recent studies extend BO also for more complex input space, such as combinatorial structures (Baptista and Poloczek 2018) and graph structured input space (Oh et al. 2019).

A framework for Bayesian optimization
Bayesian optimization framework has two key components (see Fig. 7). The first component is a probabilistic surrogate model, which consists of a prior distribution that models the unknown objective function. The second component is an acquisition function that is optimized for deciding where to sample next. The surrogate model provides a posterior probability distribution that describes potential values for H (θ ) at a candidate configuration θ . The basic idea is that each time we observe H at a new point θ , we update this posterior distribution. The most used surrogate model is the Gaussian Process (GP), which is completely specified by a mean function μ(θ ) : → R and a definite positive covariance function, also called kernel, k θ , θ : 2 → R, An important property of the kernel function is that the closer two points are in the input space, the larger is the value of the kernel function. Common choices are the Gaussian kernel or the Matern kernel. Regarding the mean function, the most common choice is a constant value, μ(θ ) μ 0 .
The BO algorithm starts with an initial set of k configurations {θ i } k i 1 and their associated function values{y i } k i 1 , withy i H (θ i ). At each iteration t ∈{k + 1,…,N}, we update the GP model using the Bayes rule, to obtain posterior distribution conditioned on the current training set S t {(θ, y i )} t i 1 containing the past evaluated configurations and observations. For any, potentially non-evaluated, configuration θ ∈ , the posterior mean μ t (θ ) and the posterior variance σ 2 t (θ ) of the GP, conditioned on S k , are known in closed-form: where K is the t × t matrix whose entries are K i, j k θ i , θ j , k(θ ) is the t × 1 vector of covariance terms between θ and {θ i } t i 1 , y is the t × 1 vector whose i th entry is y i , and τ 2 is the noise variance. If a new point θ t+1 is selected and evaluated to provide an observationy t+1 H (θ t+1 ), we add the new pair{(θ t+1 , y t+1 )} to the current training set S t , obtaining a new training set for the next iteration S t+1 S t ∪ {(θ t+1 , y t+1 )}.
The next candidate point to evaluate is selected by solving an auxiliary optimization problem, typically of the form: where U t is the acquisition function to maximize. The rationale is that, because the optimization run-time or cost is dominated by the evaluation of the expensive function f , time and effort should be dedicated to choosing a useful and informative (in a sense defined by the auxiliary problem) point to evaluate. In Fig. 8 we show the pseudocode of the BO. In BO, typical utility functions used to select the next candidate point include: the probability of improvement (Kushner 1964), the expected improvement (Mockus 1989), GP Confidence Bound (Auer 2003) (i.e., Lower Confidence Bound, LCB, in case of minimization and Upper Confidence Bound, UCB, in case of maximization),  (Frazier et al. 2008). Since the low cost to evaluate the acquisition function, greedy optimization strategies can be used to solve the problem (Eq. 14) such as random search, multi-star approach or genetic algorithms.
In this paper, we use the Expected Improvement (EI), that is based on the maximization of the following function: where N and are the normal distribution function and the normal cumulative distribution function, respectively. The EI has two components: the first can be increased by reducing the mean function μ t (θ ), the second can be increased by increasing the variance σ t (θ ). These two terms balance the trade-off between exploitation (evaluating at points with low mean) and exploration (evaluating at points with high uncertainty).

Dealing with integer parameters
The previous framework considered for BO through GP assumes that all the variables for H are continuous. However, even if the regularization factor λ and the learning rate η are continuous, the number of latent factor K takes values in a closed subset of integers. Optimization problems involving continuous and integer/discrete variables is common in many tasks of optimizing for the hyper-parameters of machine learning systems (Snoek et al. 2012). Typical examples can be found, e.g., in an ensemble of decision trees generated by the gradient boosting algorithm, in which we must adjust the learning rate and the maximum depth of the trees, or in a deep neural network, in which we must adjust the learning rate, the number of layers and the number of neurons per layer, which can only take discrete values. In this case, two possible approaches can be considered. In the first case, other surrogate models, different from GP, can be proposed, such as Random Forest (RF) (Tin Kam Ho 1995;Candelieri et al. 2018b), which should be, at least ideally, well suited for optimization problem with only integer variables or mixed space in which most of the variables are integer. In the second approach, that is the one we considered, the GP is used, but an approximation is done along the integer dimensions. More in details, the optimization of the acquisition function A can be done assuming all variables take continuous values. The values for the integer-valued variables are replaced by the closest integer (the naïve approach in "Dealing with categorical and integer-valued variables in Bayesian optimization with Gaussian processes" (Garrido-Merchán and Hernández-Lobato 2018)). Another possibility is to deal with integer-valued variables considering their discrete nature: for instance, random sampling-based optimization will sample in a finite set of integer values instead of a continuous range. This is the approach followed by the software Scikit-optimize (https://scikit-optimize.github.io/), and the one we considered. Figure 9 shows an example of optimization for a 1D integer function (red), and the approximation (mean and variance) to the original function by the GP model (green), whereas the second column shows the acquisition function (EI) values after the surrogate model is fitted. Note that only integer values are considered for the acquisition function. Fig. 10 Two images of the score function as a function of λ and K, fixing two different values of η in many different works, among which (Tsai and Hung 2012;Matuszyk et al. 2016;Katarya and Verma 2017;Bogunovic et al. 2018). One of these datasets is called MovieLens-100 k and consists of 100,000 ratings (1-5) from 943 users on 1682 movies. Each user has rated at least 20 movies. The data was collected through the MovieLens web site (https://movielens.org) during the seven-month period from September 19th, 1997 through April 22nd, 1998.
To rating matrix associated to the problem consists of 943 rows (users), 1682 column (items) and about 100,000 known entries. To apply the procedure described in Sect. 2, we use the Surprise library (https://surpriselib.com/), a Python scikit library for recommender system. This library has a set of built-in prediction algorithms, among which the matrix-factorization algorithm with a default setting of the hyperparameters, 0.02 for λ, 0.005 for η, and 100 for K , respectively. As performance function H , we consider a tenfold cross-validation procedure with 10 splits, from which we extract the mean RMSE on the different folds. The objective of this problem is to find the best configuration for hyper-parameters in as few iterations as possible. One iteration is defined as one call to the tenfold cross validation procedure. Using the default setting of the hyper-parameters, we obtain a value of 0.9296, where the rating error can be at most 4.
To search which hyper-parameter combination yields the best performance results, we define a configuration search space with three dimensions (λ, η, K ): the number of latent factors range in the integer interval {10, 100}, whereas the learning rate and the regularization parameter range in the continuous interval [0.001,0.1].
First, to analyse the complexity of the metric score function, in Fig. 10 we report two images in which we represent the score function as a function of λ and K, fixing to different value of η.
In the first case (η 0.005) the value of the score function varies between 0.925 and 0.945. The score function decrease as a function of K, and increases as a function of λ, respectively. The worst score values are for λ ≈ 0 (no regularization) and K > 30, and for λ > 0.8 and K < 30. The better score values are for 0.02 < λ < 0.04 and 30 < K < 50. In this part of the graph, we have a value of the score function less than 0.928.
In the second case (η 0.0155), the values of the score function varies between 0.905 and 1.1. For λ < 0.05, the score function decreases as a function of λ, and increases as a function of K, respectively. In this part of the graph, we have values of the score function higher than 0.93. The worst case is when λ ≈ 0 (no regularization) and K > 20, in which we have value of the score function over 1. Then for λ > 0.05, the score function appears less sensitive to the value of K and λ.
To apply BO we use the gp_minimize function from the Scikit-Optimize package (Head et al. 2019), based on the Scikit-Learn library (Pedregosa et al. 2012). This function uses GP as surrogate model with a Matern kernel (ν 2.5) and EI as acquisition function. The method to optimize the acquisition function is the Random Search, in which the function is evaluated on 10,000 points.
We perform BO several times using a different seed for the random generator. In Fig. 11, we report the results on five different experiments, in which we evaluate the objective function 30 times using BO (5 initial configurations chosen randomly, and 25 iterations of the BO algorithm). The plot shows the value of the minimum found (y-axis) as a function of the number of iterations performed so far (x-axis). The BO reaches a mean final value and standard deviation of 0.9064 and 0.0009, respectively, and the final values are all lower than the one obtained for the default setting. For the first iterations, the curves are different. After iteration five, the next point at which to evaluate the function is guided by the model, which is where an important decrease starts to appear. After iteration 25, all the curves seem to converge to a common value. To make the analysis of the BO results more robust, we perform a total of 50 different tests, and we report the results in Fig. 12 as a boxplot (minimum, first quartile, median, third quartile, and maximum) for each iteration. The performances are compared with that of Random Search (RS). As we can note, the behaviour of the two algorithms is similar in the first iterations (below 10/12 iterations). On the contrary, the performance for BO is better than RS in the last iterations: the BO reaches a mean final value and standard deviation of 0.9062 and 0.0007, respectively, whereas the RS reaches a mean final value and standard deviation of 0.9086 and 0.0026, respectively. One can test if the differences between the two distribution of BO and RS are meaningful or not at different iterations. This can be done performing the Mann-Whitney U test in which the null hypothesis is that "the two groups are sampled from populations with identical distribution". If the p-value is below 0.05, then the null hypothesis is not true, and the two groups are different. Performing a Mann-Whitney U test at iteration n°1°, 10°, 20°, and 30°, we obtain a p value of 0.380, 0.496, 0.066, and 3.29e-09, respectively.
Then the two distributions can be considered similar in the first iterations but are different in the final ones. This means that the performance of BO and RS are similar in the first iterations, but then BO outperforms RS. This is due mainly to the capacity of BO to explore the most promising regions obtained by the associated surrogate model, that results better than RS for the exploitation phase. Figure 13 represents the 2D scatter plot between η and λ (a), η and K (b), and λ and K(c), in which we display circles at the locations specified by the final configurations of the 50 tests for BO (blue circles) and RS (red circles). As a general result, we can see that η must be quite low and λ must be quite high, where K can take different values without altering so much the score value. We can see that many combinations of the hyper-parameters reach a low value of the score function, and that it is possible to consider a low value of K for the sake of interpretability. Indeed, the latent factors can refer, for MovieLens-100k dataset to the genre that the movie belongs to. In Fig. 14 we analyse the distribution of the points chosen by BO (left) and RS (right), respectively, for two different tests. The found minimum is represented by the black point. We report only the value of η and λ (x-axis) over the value of the objective function, fixing K 100. From the figure, we can note that the distribution of points selected by BO is more focused on the most promising zone of the score function, located in the upper-left part. On the contrary, a few points are in the low part of the graph, where we have a high value of the score function.
This general behaviour is testified also in Fig. 15, in which we represent the boxplot of η (left) and λ (right) component for the group formed by all the best configurations obtained by the 50 tests for BO (blue) and RS (red). Note that BO group is more focused on the most promising zone of the score function, located in the upper-left part, with η ≈ 0.02 and λ ≥ 0.8.

Conclusions and remarks
In these years, many approaches (content or collaborative filtering) and types of ML algorithms (memory-based methods and model-based methods) have been studied to build efficient RS. Like many applications of ML, the aim of a RS is to predict the preferences of the users on new incoming data according to a huge amount of known data (the past histories of users and/or items). As a learning algorithm, we used the CF matrix-factorization method, which leads to an optimization problem that can be solved through the stochastic gradient descent algorithm. The effectiveness of this procedure for RS goes through the tuning of its hyper-parameters. Two of these, the number of latent factors and the regularization factor, are related to objective function, whereas the learning rate is related to the optimization procedure.
A default setting of these values cannot be done without any prior information about the problem but requires a hyper-parameter optimization procedure. The associated objective function, usually a score measure based on cross-validation, is noisy time-consuming and black-box. The Grid Search method cannot assure a good precision if we use a large grid size or results too long to compute if we use a small grid size.
In this work, we have showed how the optimal hyper-parameter configuration can be obtained optimizing through BO a performance measure based on cross-validation. The success of BO strategy can be motivated by its exploration-exploitation balance: the exploration probes a larger portion of the search space to find the most promising regions, that are not yet refined, the exploitation allows of probing these promising regions in order to improve already promising solutions.
The numerical results have been obtained on a benchmark example, called Movielens-100k. The results showed that BO obtains a value of the performance metric slightly better than the RS. However, the better result obtained by BO could revealed more significant in case of complex objective or increasing the number of variables associated to different hyper-parameter optimization problem.
Funding Open access funding provided by Università degli Studi di Milano -Bicocca within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.