In machine learning, the scale of our features is a key component in building models. When we work with data from applications, we have to accept it as it is. In most cases, we cannot control where the nodes are lying. Let us, e.g., take recommendations in online shopping. We are only able to analyze the customers that actually exist and what they bought in the shop. However, the features may lie on immensely different scales. If we measure, e.g., the time a customer spent in the shop in seconds, as well as their age in years, the result will be a scale that contains values with thousands of seconds and a scale ranging from up to 90 years. Bringing those features on similar scales trough normalization may significantly improve performance of our model.

Two common methods for data normalization are min-max-normalization and Z-score normalization, see, e.g.,[1]. The former method will yield data in the interval [0, 1] and is especially useful if there is an intrinsic upper and lower bound for the values, e.g., when considering age. If we come back to our previous example, the time a customer spends in the shop would be less suitable since the values may have a wide range and we will probably have very few people with a significantly small or large time. In this case, the Z-score normalization makes much more sense. It tells us how many standard deviations our value lies away from the mean of the data resulting in a distribution with zero mean and variance one.

The explainable ANOVA approximation method introduced in [2–4] is based on the well-known multi-variate analysis of variance (ANOVA) decomposition, see, e.g,.,[5–10], and relies on the existence of a complete orthonormal system in the space which is suitable for fast matrix-vector multiplication algorithms in grouped transformations, c.f.[11]. Until now, this method was always applied with min-max-normalization since it relied on the space L2([0,1]d) of square-integrable functions over the cube with the half-period cosine basis. It is our goal to modify the approach in order to create the possibility to work with standard normal distributed data, i.e., data that has been obtained trough Z-score normalization.

We aim to achieve this by using the transformation ideas from [12] and [13] in order to construct a complete orthonormal system in the space

with the probability density of the standard normal distribution

The cumulative distribution function of the standard normal distribution is given by

(see Figure 1) for visualization, with the error function defined as

Combining this transformation with the half-period cosine basis allows for fast multiplications in the grouped transformations and makes the ANOVA approximation method applicable for Z-score normalized data.

As an example, we apply this approach to a dataset about the detection of forest fires, see [14, 15]. Constructing a model with the capability of efficiently predicting the size of the fire in this dataset may provide a way of predicting the occurrences of fires. This creates the possibility of efficiently implementing appropriate counter-measures. In our time of climate change with massive forest fires every year, e.g., in Australia or the USA, this is an extremely current topic. With the interpretation capabilities of the ANOVA method, c.f.[4], we are additionally able to explain the importance of our features and give reasonable explanation for the predictions.

In this section, it is our goal to construct a complete orthonormal system in the space L2(ℝd,ω) from Equation (1) with the product density ω(x) from Equation (2). This is the probability density function of the standard normal distribution, i.e., the normal distribution with zero mean and variance one. We have ∫ℝdω(x)dx=1, as well as supx∈ℝdω(x)=(2π)-d/2 which implies ω∈L∞(ℝd).

We aim to construct the basis using transformation ideas from [12, 13] and the half-period cosine basis on L2([0,1]d). The orthonormal basis functions on L2([0,1]d) are given by

with ||k||₀ := |supp k| for supp k := {s ∈ {1, 2, …, d}:k_s ≠ 0} and |supp k| the cardinality of supp k. We start from a given function f:ℝ^d → ℝ, f∈L2(ℝd,ω), and aim to transform it onto the cube [0, 1]^d. As transformation from the interval [0, 1] to ℝ, we apply the inverse cumulative distribution function Φ⁻¹ in each variable to obtain

with the inverse transformation being

The transformation is related to inverse transform sampling, see, e.g.,[16]. As a result, we have the commutative diagram in Figure 2. This allows us to transform the half-period cosine to a complete orthonormal system on L2(ℝd,ω) with the help of the following lemma.

Lemma 2.1. Let g,h∈L2([0,1]d), u,v∈L2(ℝd,ω) with probability density ω from Equation (2), and transformation ψ, ψ⁻¹ as in Equations (5) and (6), respectively. Then

and subsequently ||h||L2([0,1]d)=||h∘ψ-1||L2(ℝd,ω) and ||u∘ψ||L2([0,1]d)=||u||L2(ℝd,ω).

Proof. Let g,h∈L2([0,1]d) and u,v∈L2(ℝd,ω). Then we insert the definition and perform a change of variables as follows

As functional determinant we obtain |detψ′(t)|=∏i=1d2πeerf-2(2ti-1) and subsequently

This proves the first equality. For the second equality, we use an analogous procedure.

Lemma 2.1 is not surprising since we have based the transformation on the cumulative distribution function Φ from Equation (3). In the following, we obtain the new orthonormal system.

Theorem 2.2. The functions (ϕktrafo)k∈ℕ0d with

form a complete orthonormal system in L2(ℝd,ω).

Proof. We have that Φ from Equation (3) is a bijection and Lemma 2.1 shows that f↦f ∘ ψ ⁻¹ is an isometric isomorphism between L2(ℝd,ω) and L2([0,1]d). Now, an isometric isomorphism between two spaces maps an orthonormal basis in one space to an orthonormal basis in the other. We obtain that (ϕktrafo)k∈ℕ0d is an orthonormal basis in L2([0,1]d).

In summary, we have constructed a complete orthonormal system (ϕktrafo)k∈ℕ0d on the weighted space L2(ℝd,ω) using transformation ideas from [13] and the well-known half-period cosine basis (ϕkcos)k∈ℕ0d on L2([0,1]d). The basis functions are visualized in Figure 3.

In this section, we briefly summarize the interpretable ANOVA (analysis of variance) approximation method and the idea of grouped transformations, see [2, 11]. The approach was considered for periodic functions, but has been expanded to non-periodic functions in [3, 4]. In this paper, we focus on functions f : ℝ^d → ℝ from L2(ℝd,ω) with probability density ω from Equation (2). Since ω is the density of the standard normal distribution, this function space is of a high relevance. It allows us, e.g., to work with data from applications that has been Z-transformed, i.e., data with zero mean and variance one, see, e.g.,[1]. Since the transformed half-period cosine (ϕktrafo)k∈ℕ0d, see Theorem 2.2, is a complete orthonormal system in the space L2(ℝd,ω), we have

and through Parseval's identity ||f||L2(ℝd,ω)2=∑k∈ℤd|ck(f)|2.

The classical ANOVA decomposition, c.f. [5–7, 9], provides us with a unique decomposition in the frequency domain as shown in [2]. We denote the set of coordinate indices with [d] = {1, 2, …, d}. The ANOVA terms are defined as

The function can then be uniquely decomposed as

into |P([d])|=2d ANOVA terms where P([d]) is the power set of [d]. Here, the exponentially growing number of terms shows an expression of the curse of dimensionality in the decomposition.

It is our goal to obtain information on how important the ANOVA terms f_u are with respect to the function f. In order to measure this, we define the variance of a function f as

Note that, we have the special case σ2(fu)=||fu||L2(ℝd,ω)2, u ⊆ [d]. The relative importance with respect to f is then measured via global sensitivity indices (GSI) or Sobol indices, see [7, 17, 18], defined as

From the GSI we get a motivation for the concept of effective dimensions, specifically the superposition dimension as one notion of effective dimension. For a given α ∈ [0, 1] it is defined as

The superposition dimension d^(sp) tells us that we can explain the α-part of the variance of f by terms f_u with u ≤ d^(sp).

Using subsets of ANOVA terms U⊆P([d]), it is our goal to find a way to circumvent the curse of dimensionality for efficient approximation. In order to achieve this, we aim to truncate the ANOVA decomposition by taking only the ANOVA terms in U into account. The truncated ANOVA decomposition is then defined as

A specific idea for the truncation comes from the superposition dimension d^(sp) in Equation (10). The idea is to take only variable interactions into account that contain d_s or less variables, i.e., the subset of ANOVA terms is

Here, we call d_s a superposition threshold. Since d_s does not necessarily have to coincide to the superposition dimension d^(sp), we call it superposition threshold. A well-known fact from learning theory is that the number of terms in U(d, d_s) grows only polynomially in d for fixed d_s<d, i.e.,

which has reduced the curse of dimensionality.

In the following, we argue why the truncation by a superposition threshold d_s works well in relevant cases. For the approximation of functions that belong to a space Hs(ℝd,ω)⊆L2(ℝd,ω) that characterizes the smoothness s>0 by the decay of the basis coefficients c_k(f), we can show upper bounds on the superposition dimension d^(sp) for α ∈ [0, 1], see, e.g.,[2]. In fact, there are types of smoothness that are proven to yield a low upper bound for the superposition dimension specifically dominating-mixed smoothness with POD (product and order-dependent) weights, c.f. [2, 19–22].

In terms of real data from applications, the situation is much different. Here, we cannot make the assumption that in complete generality we have a low superposition dimension. However, there are many application scenarios where numerical experiments successfully showed that this is indeed the case, see, e.g., [5]. Since we generally do not have a-priori information, we work with low superposition thresholds d_s for truncation and validate on our test data.

We now apply the previously described method to the dataset [14] from the UC Irvine machine learning repository. The dataset contains information about forest fires in the Montesinho national park in the Trás-os-Montes northeast region of Portugal. The data was collect from 2002 to 2003. Specifically, we have d = 12 attributes about the fires and the target variable is the area of the forest that was destroyed by it. If we obtain an efficient model, it can be possible to predict the risk for a future forest fire using parameters that can be easily measured. This information can then be used to prepare appropriate countermeasures. The dataset has been thoroughly considered in [15] and we compare to the results they obtained. The authors considered several different methods: the naive average predictor (Naive), multiple regression (MR), decision tree (DT), random forest (RF), neural network (NN), and support vector machine (SVM).

We group the 12 attributes into 4 categories as in [15], i.e., spatial, temporal, FWI system, and meteorological data, (see Table 1). The spatial attributes describe the spatial location of the fire in a 9 by 9 grid of our considered region. The temporal attributes are the month of the year and the day of the week when the fire occurred. The forest fire weather index (FWI), c.f.[28], is the Canadian system for rating fire danger and the dataset collects several components of it, see also Table 2. Moreover, four meteorological attributes which are used by the FWI index were selected. The target variable describes the area that was burned by the fire.

In terms of pre-processing, we apply a Z-score transformation to the variables and the logarithmic transformation log(1+·) to the burned area. The Z-score transformation achieves that our data has zero mean and unit variance. The logarithmic transformation on the target is opportune since it shows a positive skew with a large number of fires that have a small size. We denote the data (X,y) with X={x1,x2,…,xM}⊆ℝ12, M = 517, and y ∈ ℝ^M. In the following, we do not use all of the variables, but build models based only on some groups as denoted in Table 1, e.g., STM says that we use spatial, temporal, and meteorological attributes without the FWI.

Table 3 shows the overall results of our experiment (ANOVA) combined with the benchmark data from [15]. Each value, our ANOVA results, as well as the others, were obtained by averaging over executing a 10-fold cross-validation 30 times. This results in a total of 300 experiments. We used a superposition threshold of d_s = 2, c.f. Equation (11), and, therefore, needed to detect optimal choices for the parameters N₁ and N₂ from Equation (12), (see Table 4). Every experiment utilized 90% of the data as training set (Xtrain,ytrain) and 10% of the data as test set (Xtest,ytest). The best performing model was selected based on the mean absolute deviation

with y~ the predictions of our model for the data points in the test set Xtest. As a second error measure, we use the root mean square error.

We are able to outperform the previously applied method for every subset of attributes in both MAD and RMSE error. Notably, the difference in the RMSE that penalizes larger deviations in the burned area stronger than the MAD is much more significant.

While we replicated the setting of [15] for benchmark purposes, it remains our goal to identify the most important attributes for the detection of forest fires. Therefore, we now use all 12 attributes of the dataset in obtaining our approximation and subsequently interpret the results. Figure 4 shows the attribute ranking r(i), i = 1, 2, …, 12, and Figure 5 the global sensitivity indices ϱ(u,S(Xtrain,I(U(12,2)))f), u ∈ U(12, 2), after computing for an approximation with N₁ = N₂ = 2 and λ = 1.0.

The attributes 3, 7, and 9 are clearly the most important. They represent the month of the year (3), the DC code of the FWI (7), and the outside temperature (9). Using only these three attributes and superposition threshold d_s = 2, we computed an approximation with N₁ = 2, N₂ = 10, and λ = e⁸. The resulting model yielded an MAD of 12.64 and an RMSE of 45.57 with 30 times of 10-fold cross validation as before. In summary, we know that the most important information of our problem is contained in only three attributes and we also obtained a better performing model using only those three attributes.

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Both authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

DP acknowledges funding by Deutsche Forschungsgemeinschaft (German Research Foundation)—Project–ID 416228727—SFB 1410. MS was supported by the German Federal Ministry of Education and Research grant 01|S20053A. The publication of this article was funded by Chemnitz University of Technology.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.