{'fraction_non_alphanumeric': 0.03148378228343335, 'fraction_numerical': 0.012662243252531126, 'mean_word_length': 4.6277893484082115, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interpretable Machine Learning with an Ensemble of Gradient Boosting Machines

A method for the local and global interpretation of a black-box model on the basis of the well-known generalized additive models is proposed. It can be viewed as an extension or a modification of the algorithm using the neural additive model. The method is based on using an ensemble of gradient boosting machines (GBMs) such that each GBM is learned on a single feature and produces a shape function of the feature. The ensemble is composed as a weighted sum of separate GBMs resulting a weighted sum of shape functions which form the generalized additive model. GBMs are built in parallel using randomized decision trees of depth 1, which provide a very simple architecture. Weights of GBMs as well as features are computed in each iteration of boosting by using the Lasso method and then updated by means of a specific smoothing procedure. In contrast to the neural additive model, the method provides weights of features in the explicit form, and it is simply trained. A lot of numerical experiments with an algorithm implementing the proposed method on synthetic and real datasets demonstrate its efficiency and properties for local and global interpretation.


Introduction
Machine learning models play an important role for making decision and inferring predictions in various applications. However, a lot of machine learning models are regarded as black boxes, and they cannot easily explain their predictions or decisions in a way that humans could understand. This fact contradicts with a requirement of understanding results provided by the models, for example, a doctor has to have an explanation of a stated diagnosis and has to understand why a machine learning model provides the diagnosis in order to choose a preferable treatment [26], the doctor cannot effectively use predictions provided by the models without their interpretation and explanation. To cope with the problem of interpretation of black-box models, a lot of interpretable models have been developed to explain predictions of the deep classification and regression algorithms, for example, deep neural network predictions [4,24,38,52].
Interpretation methods can be divided into two groups: local and global methods. Local methods try to explain a black-box model locally around a test example whereas global methods derive interpretations on the whole dataset or its part. In spite of importance and interest of the global interpretation methods, most applications aim to understand decisions concerning with an object to answer the question what features of the analyzed object are responsible for a black-box model prediction. Therefore, we focus on both the groups of interpretation methods, including local as well as global interpretations.
One of the most popular post-hoc approaches to interpretation is to approximate a black-box model by the linear model. The well-known methods like SHAP (SHapley Additive exPlanations) [36,49] and LIME (Local Interpretable Model-Agnostic Explanation) [43] are based on building a linear model around the instance to be explained. Coefficients of the linear model are interpreted as the feature's importance. The linear regression for solving the regression problem or logistic regression for solving the classification problem allow us to construct the corresponding linear models.
However, the simple linear model cannot adequately approximate a black-box model in some cases. Therefore, its generalization in the form of Generalized Additive Models (GAMs) [25] is used. GAMs are a class of linear models where the outcome is a linear combination of some functions of features. They aim to provide a better flexibility for the approximation of the black-box model and to determine the feature importance by analyzing how the feature affects the predicted output through its corresponding function [34,55]. GAMs can be written as follows: y(x) = w 0 + w 1 g 1 (x 1 ) + ... + w m g m (x m ). (1) where x = (x 1 , ..., x m ) is the feature vector; y is the target variable; g i is a univariate shape function with E(g i ) = 0; w = (w 1 , ..., w m ) is the vector of coefficients.
The feature contribution to the black-box model output can be understood by looking at the shape functions g i [40,10].
One of the interpretation methods using GAMs is Explainable Boosting Machine (EBM) proposed by Nori et al. [40]. According to EBM, shape functions are gradient-boosted ensembles of bagged trees, each tree operating on a single variable. Another interesting class of models called Neural Additive Models (NAMs) was proposed by Agarwal et al. [2]. NAMs learn a linear combination of neural networks such that a single feature is fed to each network which performs the function g i (x i ). The networks are trained jointly. The impact of every feature on the prediction is determined by its corresponding shape function g i . Similar approaches using neural networks for constructing GAMs and performing shape functions called GAMI-Net and the Adaptive Explainable Neural Networks (AxNNs) were proposed by Yang et al. [54] and Chen et al. [11], respectively.
One of the popular models to explain the black-box model by means of GAMs is the well-known gradient boosting machine (GBM) [19,20]. GBMs have illustrated their efficiency for solving regression problems. An idea behind their using as interpretation models is that all features are sequentially considered in each iteration of boosting to learn shape function for all features. For example, it is shown by Lou et al. [34] that all shape functions are set to zero to initialize the algorithm. For each feature, residuals are calculated, the onedimensional function is learned on residuals, and it is added to the shape function. This procedure is performed over a predefined number of iterations. It should be noted that regression trees are often used to predict new residuals. They are built in the GBM such that each successive tree predicts the residuals of the preceding trees given an arbitrary differentiable loss function [46].
Another modification of the GBM for interpretation is based on the boosting package XGBoost [12]. Chang et al. [10] propose to limit the tree depth to 1 to avoid interactions of features and to convert XGBoost to the GAM.
In order to extend the set of interpretation methods using GAMs, we propose another method based on applying GBMs with partially randomized decision trees [32]. In contrast to extremely randomized trees proposed by Geurts et al. [23], the cut-point for partitioning each feature in these trees is determined randomly from the uniform distribution, but the best feature is selected such that it maximizes the score. The proposed method is based on a weighted sum of separate GBMs such that each GBM depends on a single feature.
In contrast to the NAM [2], the proposed model calculates coefficients w 1 , ..., w m of the GAM in the explicit form because GBMs are built in parallel, but not sequentially. Moreover, the problem of a long training of neural networks is partially solved because the proposed algorithm implementing the model optimally fits coefficients w 1 , ..., w m in each iteration of the gradient boosting due to using the Lasso method and a cross-validation procedure. When the weights stop changing, the training process can be stopped, unlike EBM, where the stopping criterion is an accuracy measure. In addition, if there exist correlations between features, the proposed model detects them. In this case, the weights continue to grow and will not reach a plateau.
In contrast to models proposed in [8], the proposed method builds many decision trees for every feature. Moreover, all GBMs are trained simultaneously and dependently. First, coefficients of the GAM play the role of an adaptive learning rate [32], which allows us to reduce the overfitting problem. Second, if the prediction is a non-linear function, for example, in the case of classification, then derivatives of the loss function for each GBM depend on predictions of other GBMs. It is important to note that the GBMs in the proposed method can be used jointly with neural networks when a part of features is processed by trees, another part is used by networks.
A lot of numerical experiments with an algorithm implementing the proposed method on synthetic and real datasets demonstrate its efficiency and properties for local and global interpretation.
The code of the proposed algorithmn can be found in https://github.com/andruekonst/egbm-gam. The paper is organized as follows. Related work is in Section 2. An introduction to the GBM, brief introductions to LIME and the NAM are provided in Section 3 (Background). A detailed description of the proposed interpretation method and an algorithm implementing it are provided in Section 4. Numerical experiments with synthetic data and real data using the global interpretation are given in Section 5. Numerical experiments with the local interpretation are given in Section 6. Concluding remarks can be found in Section 7.

Related work
Local and global interpretation methods. Due to importance of the machine learning model interpretation in many applications, a lot of methods have been proposed to explain black-box models locally. One of the first local interpretation methods is the Local Interpretable Model-agnostic Explanations (LIME) [43], which uses simple and easily understandable linear models to approximate the predictions of black-box models locally. Following the original LIME [43], a lot of its modifications have been developed due to a nice simple idea underlying the method to construct a linear approximating model in a local area around a test example. Examples of these modifications are ALIME [47], Anchor LIME [44], LIME-Aleph [42], GraphLIME [27], SurvLIME [33]. Another explanation method, which is based on the linear approximation, is the SHAP [36,49], which takes a game-theoretic approach for optimizing a regression loss function based on Shapley values. A comprehensive analysis of LIME, including the study of its applicability to different data types, for example, text and image data, is provided by Garreau and Luxburg [21]. The same analysis for tabular data is proposed by Garreau and Luxburg [22]. An interesting information-theoretic justification of interpretation methods on the basis of the concept of explainable empirical risk minimization is proposed by Jung [29].
An important group of interpretation methods is based on perturbation techniques [17,18,41,51]. The basic idea behind the perturbation techniques is that contribution of a feature can be determined by measuring how a prediction score changes when the feature is altered [14]. Perturbation techniques can be applied to a black-box model without any need to access the internal structure of the model. However, the corresponding methods are computationally complex when samples are of the high dimensionality.
Global interpretations are more difficult tasks, and there are a few papers devoted to solving the global interpretation tasks [6,28,30,35,37,53]. We also have to point out models which can be used as local and global interpretations [2,8,10,11,34,40,54].
Gradient boosting machines with randomized decision trees. GBMs are efficient tool for solving regression and classification problems. Moreover, it is pointed out by Lundberg et al. [35] that GBMs with decision trees as basic models can be more accurate than neural networks and more interpretable than linear models. Taking into account this advantage of decision trees in GBMs, a modification of the GBM on the basis of the deep forests [56] called the multi-layered gradient boosting decision tree model is proposed by Feng et al. [16]. Another modification is the soft GBM [15]. It turns out that randomized decision trees [23] may significantly improve GBMs with decision trees. As a results, several models are implemented by using different modifications of randomized decision trees [31]. These methods can be successfully applied to the local and global interpretation problems.

A brief introduction to the GBM for regression
The standard regression problem can be stated as follows.
Given N training examples D = {(x 1 , y 1 ), ..., (x N , y N )}, in which x i belongs to a set X ⊂ R m and represents a feature vector involving m features, and y i ∈ R represents the observed output or the target value such that y i = q(x i ) + ε. Here ε is the random noise with expectation 0 and unknown finite variance. Machine learning aims to construct a regression model or an approximation g of the function f that minimizes the expected risk or the expected loss function with respect to the function parameters. Here P (x, y) is a joint probability distribution of x and y; the loss function L(·, ·) may be represented, for example, as follows: There are many powerful machine learning methods for solving the regression problem, including regression random forests [5,7], the support vector regression [48], etc. One of the powerful methods is the GBM [20], which will be briefly considered below.
GBMs iteratively improve the predictions of y from x with respect to L by adding new weak or base learners that improve upon the previous ones, forming an additive ensemble model of size M : where i is the iteration index; h i is the i-th base model, for example, a decision tree; γ i is the coefficient or the weight of the i-th base model. The algorithm aims to minimize the loss function L, for example, the squared-error L 2 -loss, by iteratively computing the gradient in accordance with the standard gradient descent method. If to say about decision trees as the base models, then a single decision tree is constructed in each iteration to fit the negative gradients. The function h i can be defined by parameters θ i , i.e., h i (x) = h(x, θ i ). However, we will not use the parametric representation of the function.
The gradient boosting algorithm can be represented be means of Algorithm 1.

Algorithm 1 The GBM algorithm
Require: Training set D; the number of the GBM iterations T Ensure: Predicted function g(x) 1: Initialize the zero base model g 0 (x), for example, with the constant value.
Calculate the residual q (t) i as the partial derivative of the expected loss function L(y i , g t (x i )) at every point of the training set, t = 1, ..., N , (the negative gradient)

4:
Build a new base model (a regression tree Find the best gradient descent step-size γ t : Figure 1: A scheme of the NAM given in [2] The above algorithm minimizes the expected loss function by using decision trees as base models. Its parameters include depths of trees, the learning rate, the number of iterations. They are selected to provide a high generalization and accuracy depending on an specific task.
The gradient boosting algorithm is a powerful and efficient tool for solving regression problems, which can cope with complex non-linear function dependencies [39].

Neural Additive Models
The proposed interpretation method can be regarded as a modification of the NAM because this method takes some ideas behind the NAM, namely, the separate training on single features and summing the shape functions g i . Therefore, we consider it in detail. According to [2], a general scheme of the NAM is shown in Fig. 1. The separate networks with inputs x 1 , ..., x m are trained jointly using backpropagation. They can realize arbitrary shape functions because there are no restrictions on these networks. However, a difficulty of learning the neural network may be a possible small amount of training data in case of the global interpretation. In this case, the network may be overfitted. On the one hand, to cope with the overfitting problem, the separate networks can be constructed as small as possible to reduce the number of their training parameters. On the other hand, simple neural networks may have difficulties in realizing the corresponding shape functions.

LIME
It is important to point out that methods based on using GAMs and NAMs [2,8,10,11,34,40,54] as well as the proposed method can be regarded as some modifications of LIME for local interpretation to some extent. In the case of local interpretation, the methods have to use the perturbation technique and to minimize the loss function which measures how the interpretation is close to the prediction. The same elements are used in LIME [43]. Therefore, we briefly consider the original LIME. LIME aims to approximate a black-box interpretable model denoted as f with a simple linear function g in the vicinity of the point of interest x, whose prediction by means of f has to be interpreted, under condition that the approximation function g belongs to a set of linear functions G. According to LIME, a new dataset consisting of perturbed examples x k around point x is generated to construct the function g, and predictions corresponding to the perturbed examples are obtained by means of the black-box model, i.e., by computing f (x k ) for many points x k . Weights w x are assigned to new examples in accordance with their proximity to the point of interest x by using a distance metric, for example, the Euclidean distance.
An interpretation (local surrogate) model is trained on new generated samples by solving the following optimization problem: Here L is a loss function, for example, mean squared error, which measures how the model g is close to the prediction of the black-box model f ; Φ(g) is the regularization term.
The obtained local linear model g interprets the prediction by analyzing coefficients of g. It should be noted that LIME can also apply GAMs as an extension of G. However, the shape functions of features have to be known for constructing the optimization problem.

The interpretation method and the algorithm for regression
The basic idea behind the proposed interpretation method is to realize a weighted sum of GBMs such that each GBM depends on a single feature. Moreover, weights of features obtained by means of the proposed algorithm can be used for interpretation because their absolute values can be viewed as impacts of the features.
The proposed method of interpretation can play roles of the local interpretation when a test example x is interpreted and a new dataset consisting of perturbed samples at a local area around the test example x is generated, as well as the global interpretation when the interpretation model is built on the entire dataset. The algorithm implementing the proposed interpretation method for regression is based on the iterative use of a GBM ensemble consisting of m parallel GBMs such that each GBM deals with a separate feature. The corresponding model is composed as a weighted sum of separate GBMs such that each GBM depends on a single feature. Partially randomized decision trees [32] are used for implementing each GBM. The cut-point for partitioning each feature in these trees is randomly selected from the uniform distribution, but the optimal feature is selected to maximize the score. The tree depth is limited to 1. It is important to point out that arbitrary trees can be used for implementing the GBMs. However, various experiments show that the used depth 1 significantly reduces the risk of overfitting. Moreover, the training time is significantly reduced when the depth of trees is 1. It is a very interesting and unexpected observation that trees of depth 1 with the random splitting of each feature provide better results than deep trees. This observation again illustrates the strength of the GBM as a machine learning model.
Initially, each GBM computes functions g (s) (x i,k ) of the k-th feature in the s-th iteration of the whole algorithm by using residuals r k . In order to assign weights to vectors g (s) k , k = 1, ..., m, the Lasso method [50] is used. It should be noted that other methods for assigning the weights can be used, for example, the ridge regression or the elastic net regression. However, our numerical experiments show that the Lasso method significantly reduces time (and the number of the ensemble iterations) for convergence of the weights. Denote weights computed by using the Lasso method as v (s) k . In order to smooth the weight changes in each iteration, we propose to update the weights by means of the following rule: Here α ∈ (0, 1] is the smoothing parameter. In particular, if α = 1, then "old" weights w In contrast to the standard GBM, we have to get residual r (s) i,k for every example as well as every feature because every feature is separately processed by the corresponding GBM. It should be noted that there holds Hence, we get The above implies that the training set for learning the k-th GBM at the next iteration is represented by pairs (x i,k , r One of the problems is that the obtained weights may be of different scales. In order to overcome this problem and to correct the weights, coefficients are multiplied by the standard deviation of every feature computed with using all points from the dataset. The obtained corrected weights will be viewed as the feature importance. The same procedure has been performed by Chen et al. [11].
Several rules for stopping the algorithm can be proposed. The algorithm can stop iterations when the weights do not change or they insignificantly change with some relative deviation . However, this rule may lead to an enormous number of iteration when weights are not stabilized due to overfitting problems. Therefore, we use just some predefined value T of iterations, which can be tuned in accordance with the weight changes.
Algorithm 2 can be viewed as a formal scheme for computing the weights and shape functions.
In fact, we have the gradient boosting algorithms in the "global" boosting algorithm which computes weights and shape functions w k = w (T ) k ; g k (x) = g (T ) k for every feature. An outline scheme of this "global" boosting algorithm is shown in Fig. 3. It conditionally can be represented as a series of the gradient boosting machine ensembles (EGBMs), where every EGBM is the algorithm given in Fig. 2   k stabilization means that there is a strong correlation between features, which leads to overfitting of the GBMs.

. Each EGBM uses the initial
The interpretation algorithm for classification is similar to the studied above algorithm. However, in contrast to the regression problem, we have softmax function and the cross-entropy loss for computing residuals. Moreover, it also makes sense to pre-train each GBM in classification to predict the target value by means of performing several iterations. The pre-training iterations aim to separate points by a hyperplane at the beginning of the training process, otherwise all points may be equal to 0.

Numerical experiments with synthetic data
In order to study the proposed interpretation method, we consider several numerical examples for which training examples are randomly generated.

Linear regression function
Let us suppose that the studied regression function is known, and it is There are 7 features whose impacts should be estimated. We generate N = 1000 points of ε with expectation 0 and standard deviation 0.05. It is assumed that outputs of the black-box model correspond to these points. We expect to get the same relationship between weights obtained by means of the explanation model and coefficients of the above linear regression function. Fig. 4 illustrates how the weights are changed with increase of the number of iterations T . One can see that the weights tend to converge. This is a very important property which means that there is no overfitting of the explanation model. Moreover, one can see from Fig. 4 that the weights totally correspond to coefficients of the linear regression. In particular, weights of features 5, 6, 7 tend to 0, weights of features 1, 2, 3, 4 tend to 18.77, 21.77, 9.16, 11.37, respectively. Fig. 5 shows changes of every feature, where large points (the wide band) are generated points of the dataset. They can be regarded as true values. The narrow band in the middle of the wide band corresponds to predictions of a single GBM in accordance with the feature which is fed to this GBM. The x-axis corresponds to the feature values, the y-axis is predictions of a single GBM.
The obtained weights are (18.77, 21.77, 9.16, 11.37, 1.24, 1.01, 0.37). These weights are given before their correction. After the correction, we get the feature importance values (0.42, 0.86, 0.08, 0.12, 0, 0, 0). It can be seen from the numerical results that relationships between interpreted weights and between coefficients of the regression function are quite the same. This implies that the proposed method provides correct interpretation results for the linear function.

Non-linear regression function
Let us return to the previous numerical example, but the seventh term in (13) is g 7 (x) = 100 (x 7 − 0.5) 2 instead of 0x 7 now. We again generate N = 1000 points with the same parameters of the normal distribution for ε (expectation 0 and standard deviation 0.05).
The obtained weights are (17.68, 22.31, 7.49, 9.29, 0.37, 0.74, 21.07). These weights are presented before their correction. After the correction, we get the feature importance values (0.28, 0.58, 0.04, 0.07, 0, 0, 0.75). It can be seen from the numerical results that the seventh feature x 7 (its function g 7 (x)) has the largest weight. Moreover, relationships between interpreted weights and between coefficients of the regression function are quite the same. This implies that the proposed method provides correct interpretation results. It is interesting to point out that   the interpretation linear model used, for example, in LIME, could not detect the non-linearity in many cases. Fig. 6 shows how the weights are changed with increase of the number of iterations T . We again observe the convergence of weights to values given above. Fig. 7 is similar to Fig. 5, and we clearly observe how the seventh feature impacts on predictions.

The chess board
So far, we have considered interpretation of the regression models. The next numerical example illustrates the binary classification task. Training examples having two labels −1, 1 and consisting of two features are generated such that the corresponding points depict a small chess board as it is shown in Fig. 8. Black and white points belong to classes y = 1 and y = −1, respectively. This is a very interesting example which illustrates inability of the proposed algorithm to correctly classify the given dataset due to overfitting of GBMs. Indeed, Fig. 9 shows that weights of features are not converged and are not stabilized with increase of the iteration number T . The overfitting problem arises because there is a strong correlation between features. The same can be seen in Fig. 10, where outputs of the black-box model are depicted in the form of two lines with y = 1 and y = −1.
Outcomes of the proposed algorithm are irregularly located because the algorithm is subjected to overfitting. It is also interesting to note that the points are scattered such that the sum of GBM outcomes for all features provides correct values. It should also be pointed out that the standard GBM with decision trees of depth 1 cannot solve this problem. This implies that the proposed parallel architecture of the tree training is better in comparison with the serial architecture. Fig. 11 illustrates predicted values of the trained explanation model.
The computed weights of features before their correction are (29.21, 30.46). After correction, they are (0.52, 0.52). In this numerical example, weights do not have a specific sense, but one can see that the correction method itself provides accurate results because both the features are identical and equal to 0.52.

The polynomial regression with pairwise interactions
The next numerical example uses the following regression function for interpretation: It contains pairwise interactions of features (x 1 x 2 and x 3 x 4 ). Fig. 6 shows how the weights are changed with increase of the number of iterations T . The stabilization of weights is observed in Fig. 6. This is an interesting fact because many features are interacted. In spite of this interaction, the proposed algorithm copes with this problem and provides the stabilized weights. Fig. 13 shows how predictions of each GBM depend on values of       Fig. 15 in order to avoid a confused mixture of many curves.
We plot individual shape functions g k (x k ) also only for four important features: RM, LSTAT, CRIM, B. They are depicted in Fig. 16. Contributions of features (y-axis) are biased to 0 and scaled in order to have the same interval from 0 to 1 for all plots. The shape plot for the most important feature RM (the average number of rooms per dwelling) shows that contribution of the RM rises significantly with the average number of rooms. It is interesting to note that the contribution decreases when the number of rooms is smaller than 4. The shape plot for the second important feature LSTAT (% lower status of the population) shows that its contribution tends to decrease. A small increase of the contribution when LSTAT is larger than 25 may be caused by overfitting. Features CRIM and B have significantly smaller weights. This fact can be seen from plots in Fig. 16 where changes of their contributions are small in comparison with the RM and the LSTAT.

Breast Cancer dataset
The next real dataset is the Breast Cancer Wisconsin (Diagnostic). It can be found in the well-known UCI Machine Learning Repository (https://archive.ics.uci.edu). The Breast Cancer dataset contains 569 examples such that each example is described by 30 features. For classes of the breast cancer diagnosis, the malignant and the benign are assigned by classes 0 and 1, respectively. We consider the corresponding model in the framework of regression with outcomes in the form of probabilities from 0 (malignant) to 1 (benign).
The importance values of features obtained by using the proposed algorithm are depicted in Fig. 17. It can be seen from Fig. 17 that features "worst texture", "worst perimeter", "worst concave points", "worst smoothness" are of the highest importance. Fig. 18 shows how weights are changed with increase of the number of iterations T . We again consider only four important features in Fig. 18.
Individual shape functions are plotted for four important features: "worst texture", "worst perimeter", "worst concave points", "worst smoothness". Fig. 19 illustrates them. The shape plot for the "worst perimeter" shows that the probability of benign drops with increase of the worst perimeter and increases for the worst perimeter above 140. The shape plot for the second important feature "worst concave points" shows that the probability of benign decreases with decrease of the worst concave points. Features "worst texture" and "worst smoothness" have significantly smaller impact on the target probability.

Numerical experiments with synthetic data
First, we return to the chess board example and investigate a bound between two checkers. A neural network is used as a black-box regression model for interpretation. The network is trained by using the quadratic loss function. The black and white checkers are labelled by 1 and 0, respectively. The network consists of 3 layers having 100 units with the ReLU activation function and optimizer Adam. N = 1000 perturbed points are generated around point (0.35, 0.2) from the normal distribution with the standard deviation 0.025. The perturbed points are depicted in Fig. 20 (the right picture). Predictions of the neural network learned on the generated "chess board" training set are depicted in Fig. 20 (the left picture). Fig. 21 depicts the shape function of the second feature x 2 for the chess board example. It follows from Fig.  21 that the probability of the class 1 (the black checker) decreases with increase of the second feature. This is also seen from Fig. 20, where generated points "move" from the black checker to the white checker.   Fig. 14). Fig. 23 shows how the weights are changed with increase of the number of iterations T . Only four important features are shown in Fig. 23.
Shape functions are depicted for four important features NOX, AGE, DIS, RAD in Fig. 24. The shape plot for the most important feature RM shows that contribution of the RM rises significantly with the average number of rooms. It is interesting to note that the contribution decreases when the number of rooms is smaller than 4. The shape plot for the second important feature LSTAT shows that its contribution tends to decrease. A small increase of the contribution, when LSTAT is larger than 25, may be caused by overfitting. Features CRIM and B have significantly smaller weights. This fact can be seen from plots in Fig. 24 where changes of their contributions are small in comparison with the RM and LSTAT.

Breast Cancer dataset and the regression black-box model
For illustrating the local interpretation, we again use the Breast Cancer Wisconsin (Diagnostic) dataset. For classes of the breast cancer diagnosis, the malignant and the benign are assigned by labels 0 and 1, respectively. We consider the corresponding model in the framework of regression with outcomes in the form of probabilities from 0 (malignant) to 1 (benign), i.e., we use the regression black-box model based on the SVM with the RBF Figure 24: Shape functions of four important features learned on the Boston Housing dataset (the x-axis indicates values of each feature, the y-axis indicates the feature contriburion) for local interpretation kernel having parameter γ = 1/m. The penalty parameter C of the SVM is 1.0. The point of interest is determined by means of the same procedure as for the Boston Housing dataset. The perturbation procedure also coincides with the same procedure for the Boston Housing dataset. All features are standardized before training the SVM and interpreting results. Importance of features obtained by using the proposed algorithm are depicted in Fig. 25. It can be seen from Fig. 25 that features "worst symmetry", "mean concave points", "worst concavity", "worst concave points" are of the highest importance. They mainly differ from the important features obtained for the global interpretation. Fig. 26 shows how the weights are changed with increase of the number of iterations T for the above important features.
For these important features, individual shape functions are plotted in Fig. 27. The shape plot for the "worst concave points" shows that the probability of benign drops with increase of the worst concave points. The shape plot for the second important feature "worst concavity" also shows that the probability of benign significantly decreases with increase of the worst concavity. Features "worst symmetry" and "mean concave points" have significantly smaller impacts on the target probability.

Breast Cancer dataset and the classification black-box model
We consider the Breast Cancer Wisconsin (Diagnostic) dataset and corresponding model in the framework of classification with the same outcomes in the form of probabilities from 0 (malignant) to 1 (benign), i.e., we use the classification black-box model based on the SVM. Parameters of the black-box model, the point of interest and the perturbation procedure do not differ from those used for the Breast Cancer dataset with the regression black-box model. All features are standardized before training the SVM and interpreting results.
The importance values of features obtained by using the proposed algorithm are depicted in Fig. 25. It can be seen from Fig. 28 that the same features as in the example with the regression black-box, including "worst symmetry", "mean concave points", "worst concavity", "worst concave points", are of the highest importance. Fig. 29 shows how the weights are changed with increase of the number of iterations T for the above important features. Individual shape functions for these important features are plotted in Fig. 30. They are similar to those computed for the regression black-box model.

Conclusion
A new interpretation method which is based on applying an ensemble of parallel GBMs has been proposed for interpreting black-box models. The interpretation algorithm implementing the method learns a linear combination of GBMs such that a single feature is processed by each GBM. GBMs use randomized decision trees of depth 1 as base models. In spite of this simple implementation of the base models, they show correct results whose intuitive interpretation coincides with the obtained interpretation. Various numerical experiments with synthetic and real data have illustrated the advantage of the proposed method.
The method is based on the parallel and independent usage of GBMs during one iteration. Indeed, each GBM deals with a single feature. However, the architecture of the proposed algorithm leads to the idea to develop a method which could consider combinations of features. In this case, the feature correlation can be taken into account. This is a direction for further research. Another interesting direction for research is to consider a combination of the NAMs [2] with the proposed algorithm.
The proposed algorithm consists of several components, including, the Lasso method, a specific scheme for updating weights. Every component can be replaced by another implementation which could lead to better results. The choice of the best configuration for the algorithm is also an important direction for further research.
The proposed method is efficient mainly for tabular data. However, it can be also adapted to the image processing which has some inherent peculiarities. The adaptation is another interesting direction for research in future.