High Dimensional Model Representation as a Glass Box in Supervised Machine Learning
Caleb Deen Bastian Herschel Rabitz

 
        
letterpaper          





















thmTheorem
corCorollary
definition
deDefinition
reRemark
lemma
pfthmProof of Theorem
lemLemma
.tifpng.png`convert #1 `dirname #1`/`basename #1 .tif`.png
øω
δ̣
∘̧
þθ
λ
𝔼
𝕍ar
ℂov
ℙ
ℝ

ℕ

γ
R_0
S̃




























packed_enum

  
  
  





















 
 








.tiffpng.pngconvert #1 
.tiff













d


propProposition





Leb


theoremTheorem
definitionDefinition
lemmaLemma
corollaryCorollary
algAlgorithm
conjectureConjecture
probProblem
remarkRemark
costCost
assumptionAssumption
hypothesisHypothesis
goalGoal
ideaIdea
noteNot to be in the paper – just a discussion note







High Dimensional Model Representation as a Glass Box in Supervised Machine Learning
    Caleb Deen Bastian, Herschel Rabitz
    December 30, 2023
===================================================================================




Prediction and explanation are key objects in supervised machine learning, where predictive models are known as black boxes and explanatory models are known as glass boxes. Explanation provides the necessary and sufficient information to interpret the model output in terms of the model input. It includes assessments of model output dependence on important input variables and measures of input variable importance to model output. High dimensional model representation (HDMR), also known as the generalized functional ANOVA expansion, provides useful insight into the input-output behavior of supervised machine learning models. This article gives applications of HDMR in supervised machine learning. The first application is characterizing information leakage in “big-data” settings. The second application is reduced-order representation of elementary symmetric polynomials. The third application is analysis of variance with correlated variables. The last application is estimation of HDMR from kernel machine and decision tree black box representations. These results suggest HDMR to have broad utility within machine learning as a glass box representation. 







§ INTRODUCTION


In many areas of science and engineering one is typically interested in using empirical data D to understand an input-output map f from some class of systems F, often having a large number of input and output variables. We indicate the input space by (X,𝒳), the output space by (Y,𝒴), and the system (model) space as F=Y^X, and put n≡|X|∈ℕ and p≡|Y|∈ℕ respectively. Suppose D is a collection of iid random variables taking values in X× Y, denoted by D={ (x_i,y_i)}, and regard each (x_i,y_i) as a realization from the probability measure ν.

This setting can be identified to supervised machine learning <cit.>, where f is a prediction function, F is some class of candidate functions, and ν is the joint probability measure on X× Y whose independent random realizations form D. The task for supervised machine learning is solving 
    f^* = _f∈ F R(f),
 where R is a risk functional 
    R: f↦1/N∑_(x,y)∈DL(y,f(x))
 with sample size N∈ℕ and where L is a loss function, L: X× Y↦ℝ_+. The learned f^* is used to assign an element f^*(x) of Y to each x in X. In this article we consider quadratic loss function L(y,f(x))=(y-f(x))^2. 



 §.§ Motivations of representation


The input-output relation f can be represented many ways. Black box representation is solely concerned with the task of assigning an element f(x) of Y to each x in X (prediction) and includes all such functions from X into Y. In other words, the particular structure of f is irrelevant so long as it is contained in F, and notions of input variable importance and dependence of f upon the important variables are abstract. General-purpose black box learning algorithms in machine learning include neural networks, kernel machines, and decision tree models. Glass box representation restricts F to include only those functions which provide information on variable importance and dependence (explanation). Glass box representation means that F is equipped with necessary and sufficient information to interpret the model output in terms of the model input. This information is referred to as interpretative diagnostics and is defined by the variable importance and variable dependence sets I(f) and P(f) respectively, each indexed on ({1,…,n}). <cit.> discusses diagnostics for high-dimensional functions based upon projection operators, and functionals can be further defined to measure variable importance from projections. These ideas are adopted in this article to construct I(f) and P(f). 




A consideration is managing the sizes of I(f) and P(f) as the number of input variables increases, each of size 2^n. The exponential growth is known as the curse of dimensionality. For large enough n, construction cost outstrips available resources, and feasibility is achieved through simplifying assumptions on F or μ, e.g., stipulating restrictions on model form or distribution. A large body of work supports the ansatz that many high-dimensional systems especially those of a physical (or real-world) nature are accurately, if not exactly, described by representations whose sizes grow polynomially in n. 


The ansatz is elevated to theorem in a result of Kolmogorov <cit.>: his superposition theorem, an existence result, establishes that every multivariate continuous function on the unit cube X=[0,1]^n can be exactly represented using a finite number of univariate continuous functions and the binary operation of addition, 
    f(x) = ∑_i=1^2n+1ϕ_i∘(∑_j=1^nψ_ij∘ x_j),
 where ϕ_i are ψ_ij are continuous functions. These functions are highly non-smooth, limiting their utility in approximation. 

Instead of seeking a representation in terms of univariate functions, suppose a hierarchy of projections into subspaces of increasing dimensions, with the expectation that subspace contributions to output(s) rapidly diminish with increasing dimension. This is the idea behind a mathematical model called high dimensional model representation (HDMR), a finite multivariate representation that efficiently manages the curse of dimensionality and that provides structured information for input-output relationships. When applied to supervised machine learning models, HDMR is a glass box that interprets the model output in terms of the model input. 



 §.§ Related work
 HDMR is discussed as early as <cit.> in ANOVA analysis and has enjoyed extensive application to statistics. When F is the collection of symmetric functionals of iid variables HDMR is known as the Hoeffding decomposition, a fundamental object in U-statistics <cit.>. The univariate terms of HDMR are sometimes known as Hajek projections and are extensively used to establish asymptotic normality of various statistics <cit.>. Global sensitivity analysis (GSA) measures the importance of variables to F <cit.>. In GSA applications functional representation of F is avoided where instead sensitivity indices are directly estimated. <cit.> discussed HDMR as a general decomposition of F where the input space resides in ℝ^n or in the n-fold product space of arbitrary linear topological function spaces. These ideas have been further developed to F for general (non-degenerate) distributions in <cit.>, wherein HDMR is known as generalized functional ANOVA, which in turn provisions GSA in terms of structural and correlative sensitivity indices <cit.>. Sometimes HDMR is known as the Hoeffding-Sobol decomposition <cit.> or Sobol decomposition <cit.>. See <cit.> for references on HDMR's earlier history.



 §.§ Contributions
  


We discuss HDMR and provide illustrative applications to motivate its utility to supervised machine learning. First, we illustrate that HDMR can diagnose information leakage. This is demonstrated for Pearson goodness-of-fit settings for “big-data” wherein HDMR characterizes estimator efficiency loss and for popular machine learning black boxes wherein HDMR identifies biases. Second, we illustrate that HDMR characterizes the information leakage experienced by partial dependence, another interpretative diagnostic, whenever input variables are correlated. In such settings, HDMR reveals that the interpretative diagnostics functionally depend upon the distribution correlation. Third, we demonstrate that HDMR admits efficient reduced-order representations of high-dimensional models, managing the curse of dimensionality. In particular, we illustrate that the input distribution regulates the efficiency of reduced-order representation in polynomials and demonstrate this effect in a learning setting. Fourth, we demonstrate that HDMR can be applied as a wrapper method for black boxes to provision glass boxes. Fifth, we estimate HDMR of trained kernel machines or ensembles of decision trees.



 §.§ Organization
 
In section <ref> we formulate interpretative diagnostics using projection operators and functionals and related quantities. In section <ref> we discuss HDMR and provide three examples of applications: (i) to goodness-of-fit in “big-data” settings (subsection <ref>), (ii) high-dimensional modeling (subsection <ref>), and (iii) analysis of variance (subsection <ref>). In section <ref> we consider alternative interpretative diagnostics and compare to HDMR for a simple mathematical model having an analytic solution. In section <ref> we apply HDMR as a wrapper method for black boxes to convey information on variable importance and dependence, and we consider two machine learning black boxes—kernel machines (subsection <ref>) and decision trees (subsection <ref>)—on analytic and empirical (test dataset) problems. We conclude with a discussion.





§ PRELIMINARIES
 In this article we consider the model space to be square-integrable functions F=L^2(X,𝒳,μ). We formulate interpretative diagnostics for the subspaces u⊆{1,…,n} (notationally we use ℕ_n={1,…,n} throughout this article) using projection operators and functionals. For every u⊆ℕ_n let P^u be a projection operator that profiles the dependence of f on x_u and put f_u≡P^u f; f_u is said to be the variable dependence or projection of f on x_u. Let I^u: F↦ℝ_+ be a functional that measures the importance of f∈ F on x_u and put S_u≡I^u f; S_u is said to be a measure of variable importance of f on x_u.  Using {I^u} and {P^u}, we define respective variable importance and variable dependence sets I(f) ≡{I^uf} and P(f) ≡{P^uf} on (ℕ_n). 

Sometimes total and relative notions of variable importance are necessary. We define the total variable importance of f on x_i as T_i≡T^i f where T^i ≡∑_u⊃ iI^u. We define the relative variable importance of f on x_i as R_i≡R^i f where R^i f≡T^i f/ ∑_j T^j f. Using {T^i} and {R^i} we define total variable importance and relative variable importance sets T(f)={T^i f} and R(f)={R^i f}. 

Using I(f), P(f), T(f), and R(f), we can pose various questions for f as super level-sets:


  * which variables are important overall? 
    T_ϵ(f)={ T_i∈T(f) : T_i ≥ϵ}
 

  * what are the relative importances of the variables? 
    R_ϵ(f)={ R_i∈R(f) : R_i≥ϵ}.
 

  * which variables are important individually? 
    I_ϵ(f)={ S_u∈I(f) : S_u ≥ϵ, |u|=1}
 

  * which variables participate in interactions?  
    X_ϵ(f)={ S_u∈I(f) : S_u ≥ϵ, |u|>1}
 

  * how does the system depend upon the important variables and interactions?  
    P_ϵ(f)={ f_u∈P(f) : I(f)∋ S_u ≥ϵ}


  * does the system admit a reduced order representation? 
    ∃  T∈{ 1,…,n}:    f≃ f^T=∑_u: |u|≤ Tf_u

 







§ HIGH DIMENSIONAL MODEL REPRESENTATION
 Writing F=L^2(X,𝒳,ν), the HDMR of f(x)∈ F is the solution to the variational problem 
    min_uf(x) - u,     u∈𝒱_0⊕∑_i𝒱_i⊕∑_i_1<i_2𝒱_i_1i_2⊕…⊕∑_i_1<…<i_l𝒱_i_1… i_l,
 where the norm · is induced by the inner product as ·=⟨·,·⟩^1/2 and the {𝒱_u⊂ F: |u|≤ l} are subspaces having certain null integral properties. It is uniquely minimized by 
    u = (𝒫_0 + ∑_i𝒫_i + ∑_i_1<i_2𝒫_i_1i_2 + … + ∑_i_1<…<i_l𝒫_i_1… i_l) f(x)
 using the collection of non-orthogonal projection operators {𝒫_u}. Putting 
    ε_l(x) ≡ f(x) - u = f(x) - f_0 - ∑_i f_i(x_i) - ∑_i_1<i_2f_i_1i_2(x_i_1,x_i_2) - … - ∑_i_1<…<i_lf_i_1… i_l(x_i_1,…, x_i_l),
 we express 
    f(x)-u = ⟨ε_l,ε_l⟩^1/2,
 which for scalar valued functions is 
    f(x)-u=∫_Xε_l^2(x)ν(x).




The HDMR of f∈ F is written as 
    f(x) = f_0 + ∑_i f_i(x_i) + ∑_i_1<i_2f_i_1i_2(x_i_1,x_i_2) + … + f_1… n(x_1,…,x_n),
, a sum of 2^n component functions <cit.>. Although many additional component functions appear when compared to Kolmogorov's result, most of these are identically zero or insignificant for many F of practical interest: it is often observed 
    f(x) ⋍ f^T(x) = f_0 + ∑_i f_i(x_i) + ∑_i_1<i_2f_i_1i_2(x_i_1,x_i_2) + … + ∑_i_1<…<i_Tf_i_1… i_T(x_i_1,…,x_i_T)
 for T≪ n, i.e. f^T is a reduced-order representation of f. Whenever T<n and f(x)≃ f^T(x), we achieve a polynomial-scaling representation in n.  In most applications T≤3 is sufficient. In general, HDMR is expected to exhibit favorable convergence in T whenever the input variables are meaningfully defined in relation to the output variable(s), typically the case with physical systems where input and output variables correspond to observable physical states. Some properties of HDMR:


  * it is unique given μ

  * {𝒱_u} are hierarchically orthogonal and partition variance, (f)=∑_u,v(f_u,f_v)

  * {𝒱_u} maximize explained variance, i.e. no other hierarchically-orthogonal functions achieve higher explanatory variance than those belonging to {𝒱_u}

  * the elements of {𝒱_u} can be attained through the action of projection operators {𝒫_u} such that f_u(x_u)≡𝒫_uf(x), i.e. the component functions are  L^2 projections of f(x)

  * all HDMR's converge at the same order: if f^T(x) converges for μ with error 𝒪(ϵ), then f^T(x) converges for μ' with error 𝒪(ϵ), although the constant can be substantially different

  * for a set of functions { f^i(x)} obeying linear-superposition conservation laws, the corresponding HDMR's obey such for every T.

The notion of hierarchical-orthogonality guarantees the existence and uniqueness of the decomposition such that (f)=∑_u,v(f_u,f_v) <cit.>. Note that hierarchical-orthogonality is a generalization of mutual orthogonality where functions are orthogonal only to functions on nested subspaces. For example, for μ_123μ_1μ_2μ_3, hierarchical-orthogonality implies ⟨ f_12,f_1⟩=⟨ f_12,f_2⟩=0 but neither implies ⟨ f_1,f_2⟩=0 nor ⟨ f_12,f_3⟩=0.


Given the variance decomposition, sensitivity indices can be defined in terms of normalized variances and covariances. This is known as structural and correlative sensitivity analysis (SCSA) <cit.>, where SCSA defines structural, correlative, and overall sensitivity indices (𝕊_u^a,𝕊_u^b,𝕊_u) for each component function u⊆ℕ_n, 
    𝕊_u^a   ≡1/(f)(f_u)
    𝕊_u^b   ≡1/(f)∑_v: v u(f_u,f_v)
    𝕊_u   ≡𝕊_u^a+𝕊_u^b,

and these satisfy 
    ∑_u𝕊_u = 1.
 Note that when μ∏_iμ_i, the projections f_u=𝒫_u f are hierarchically-orthogonal with ∑_u𝕊_u^a<1 and ∑_u(𝕊_u^a+𝕊^b_u)=1. 


When the input variables are independent, μ=∏_iμ_i, the component functions are mutually orthogonal and can be recursively constructed, 

    f_i_1… i_l(x_i_1,…,x_i_l) ≡    M^i_1… i_lf(x) - ∑_j_1<…<j_l-1⊂{ i_1,…,i_l}f_j_1… j_l-1(x_j_1,…,x_j_l-1)
        - ∑_j_1<…<j_l-2⊂{ i_1,…,i_l}f_j_1… j_l-2(x_j_1,…,x_j_l-2) - …
       -∑_j⊂{ i_1,…,i_l}f_j(x_j) - f_0

where M^i_1… i_lf(x) ≡∫_X_-i_1… i_lf(x)∏_j∉{ i_1,…,i_l}μ_j(x_j) (note that (<ref>) does not hold for correlated μ). Because of the mutual orthogonality of the component functions, we have (f_u,f_v)=0 for u v such that ∑_u𝕊^a_u=1. This is known as global sensitivity analysis (GSA) <cit.>.
 



In regard to the projection operators {P^u} and functionals {I^u} of interpretative diagnostics, we put for every u 
    P^u_HDMR f(x) ≡𝒫_u f(x) = f^HDMR_u(x_u)
 and 
    I^u_HDMR f(x) ≡ (𝕊^a_u, 𝕊^b_u, 𝕊_u).


We provide three examples of HDMR in machine learning. 



 §.§ Eliminating information leakage in “big-data” settings
 Consider an independency of iid random variables D={ X_i} with distribution μ on X. Suppose (X,𝒳) is discrete and that D is represented as an infinite double-array [X_nj] of real-valued random variables (that is, for each n∈ℕ there is a k_n∈ℕ such that X_nj=0 for all j>k_n). The probability law of X_n is given by μ_n{ k}≡ℙ(X_n = k) for k∈{1,…, k_n}. We consider “big-data” to be the case of k_n→∞ as n→∞, i.e., the dimension of X_n (denoted by k_n) increases to infinity as the sample size (denoted by n) increases so. Consider the goodness-of-fit Pearson chi-square (PGOF) statistic χ_n^2. For fixed k_n=k the statistic χ_n^2 asymptotically follows the χ^2-distribution with (k-1) degrees of freedom and whenever k is large the standardized statistic (χ^2_n-(k-1))/√(2(k-1)) is approximated by the standard Gaussian distribution. However, k_n=k is not the case in “big-data” applications where k_n→∞ is observed as n→∞. As shown in <cit.> the Gaussian approximation may or may not be valid for the doubly infinite case; furthermore it turns out that the asymptotic behavior of χ_n^2 is characterized by a subset of its HDMR component functions. Employing the material in the appendices, we define a new statistic using this subset of component functions, denoted by χ_nHDMR^2,  
    χ_nHDMR^2(x_1,…,x_n)   ≡ f_0+∑_i<jf_ij(x_i,x_j) 
       ≡χ_n^2(x)-∑_if_i(x_i),
 and compare its relative efficiency to χ_n^2 for the power law distribution, 
    α∈[0,1),   μ_n{ k} = (C_α k^α)^-1,
 where 
    C_α=∑_k∈[k_n]k^-α≃ k_n^1-α/(1-α).
 The HDMR statistic dominates the PGOF statistic.

In Figure <ref>, we exhibit empirical distributions for PGOF, HDMR chi-square, and chi-square for k_n=2500 and n∈{50,500,5000}, having corresponding λ∈{1,10,100}, for the power law distribution with α=1/2. We show statistics for the estimators in Table <ref>. For λ=100, Figure <ref> shows both χ_n^2 and χ_nHDMR^2 are asymptotically χ^2. For λ = 10, Figure <ref> shows that χ_n^2 is not asymptotically χ^2 but χ_nHDMR^2 is. For λ = 1, Figure <ref> shows that neither χ_n^2 nor χ^2_nHDMR is asymptotically χ^2. In particular, χ^2_nHDMR is expressed in terms of a Poisson law (see appendices for a precise characterization), although observe that the empirical distribution is truncated due to undersampling of the power law tails. The HDMR statistic χ_nHDMR^2 dominates the PGOF statistic χ_n^2, i.e. PGOF experiences information leakage in comparison to the HDMR. For λ = 1 we observe PGOF to have a relative efficiency of ∼10% to that of HDMR. 









 §.§ Reduced-order representation of high-dimensional behavior
 Note that in this section we use ν to indicate the measure on X, and μ and σ denote the mean and standard deviation common to the input variables.


In this application we compute HDMR for polynomial models, in particular the function ∏_i=1^n x_i and a sum of elementary symmetric polynomials in n variables. We show that the coefficient of variation ρ of the input distribution ν regulates the efficiency of low-dimensional HDMR approximations. We illustrate these properties for polynomial models. 



Consider 
    f(x) = ∏_i=1^n x_i,   iid x,   ρ≡σ/μ 0
 with 
    f =μ^2n((1+ρ^2)^n-1).

Per (<ref>) the component functions of f(x) are 

    f_0    = 𝒫_0f(x) = μ^n
    
    f_i(x_i)    = 𝒫_if(x) =μ^n-1x_i - f_0
    
    f_ij(x_i,x_j)   =𝒫_ijf(x) = μ^n-2x_ix_j - f_i(x_i) - f_j(x_j) - f_0
    ⋮   .

Employing the material in the appendices, the sensitivity indices satisfy 
    ∑_k
    i_1<…<i_k𝕊_i_1… i_k = 1
 and at each order follow

    p{ k}≡∑_i_1<…<i_k𝕊_i_1… i_k   = nkρ^2 k/(ρ^2+1)^n-1

where ∑_k p{ k} = 1. We have 𝔼  p = n ρ^2 (ρ ^2+1)^n-1/(ρ ^2+1)^n-1 and p = n ρ ^2 (ρ ^2+1)^n-2((ρ ^2+1)^n-n ρ ^2-1)/((ρ ^2+1)^n-1)^2. Observe that when ρ<1 the explained variance of low dimensional approximations increases. For context, the uniform distribution on the unit interval (maximum entropy on [a,b]) has ρ(Unif[0,1])=√(3)/3≈ 0.58. In Figure <ref>, we plot the percent of explained variance by subspace order p≡(p{ k}: k∈ℕ_n) for n=100 and ρ∈{1/4,1/2,1,2,4}. We observe that for ρ<1 the probability mass of the probability vector p is concentrated about small k. Similar results are observed for sums of elementary symmetric polynomials in n variables 
    E_n(x) = 1+∑_k∈ℕ_n
    i_1<…<i_kx_i_1… x_i_k,
 where τ≡σ/1+μ regulates efficiency, and are contained in the appendices. 





We demonstrate this effect through Monte Carlo estimation of mean-squared-error test performance using the gradient boosting regressor (GBR) machine (Monte Carlo cross-validation), i.e. the generalization ability of the learning algorithm. We take n=10 and form D={(x_i,f(x_i))}_i=1^10^4 from iid random variables from a Beta(a,b) distribution on X=[0,1]^n. For the beta distribution, we take ρ∈{1/4,1/2,1} using (a,b)∈{(15/2,15/2),(3/2,3/2),(1/2,3/2)}. We take a tree depth of six with 5×10^3 trees. We split the data into training and test data. This is performed for 50 independent estimates, where (D_i) are independently standardized using training data. Mean and standard deviation Monte Carlo estimates are shown below in Table <ref> and density truth plots in Figure <ref> (inverse-transformed): MSE markedly increases as ρ increases.












The total variable importance of x_i to f is computed as 
    T_i=∑_kρ ^2 in-1k-1/(ρ ^2+1)^n-1 = ρ^2/ρ^2+1(ρ ^2+1)^n/(ρ ^2+1)^n-1≈ρ^2/ρ^2+1.
 Observe that T_i is independent of the value of n whenever n>k for sufficiently large k. This is depicted in Figure <ref> and observed for k∼ 100. Relative variable importance is trivially R_i=n^-1 and is independent of ρ. 










 §.§ Analysis of variance for correlated input variables
 Note that in this section and subsequent sections, we use μ to indicate the probability measure on (X,𝒳).


In this application we show how variable importance I(f) and variable dependence P(f) by HDMR depend upon correlation ρ in the input variables and this is necessary to preserve model variance. This provisions an analysis of variance that is consistent for correlated or degenerate input variables: total variance is preserved for all values of ρ.


Let 
    F=L^2(ℝ^n,ℬ(ℝ^n),μ)=Π^n≡Span{ x^α: |α|≤ n,  α∈ℕ_0^n}
 where α = (α_1,…), |α|=∑_iα_i, and μ is Gaussian. In particular, take n=2 and let 
    f(x) = β_0 + β_1x_1 + β_2x_2 + β_12x_1x_2.
 Its projections { f_u𝒫_u f} decompose f, 
    f(x)=f_0+f_1(x_1)+f_2(x_2)+f_12(x_1,x_2),

where

*!


    f_0    = β_0+β_1μ_1+β_2μ_2+β_12μ_1 μ_2+β_12ρσ_1 σ_2
    
    f_1(x_1)    = (x_1-μ_1) (β_1+β_12μ_2)+(ρ/ρ^2+1)β_12σ_2/σ_1((x_1-μ_1)^2-σ_1^2)
    
    f_2(x_2)    =(x_2-μ_2) (β_2+β_12μ_1)+(ρ/ρ^2+1)β_12σ_1/σ_2((x_2-μ_2)^2-σ_2^2)
    
    f_12(x_1,x_2)   = β_12(-ρσ_2^2 ((ρ^2-1) σ_1^2+(x_1-μ_1)^2)+(ρ^2+1) σ_1 σ_2 (x_1-μ_1) (x_2-μ_2)+ρσ_1^2 (-(x_2-μ_2)^2))/(ρ^2+1) σ_1 σ_2.


Computational details are contained in the appendices and these results are special cases of those in <cit.>.  Observe that HDMR projections contain elements not observed in f and that these additional elements are functions on X that depend upon the parameters of μ. For example, f_1(x_1) contains a quadratic term in x_1, whereas f(x) does not. The sensitivity indices of f exhibit complex dependence on the model coefficients and distribution parameters,

*!


    𝕊^a_1    = σ_1^2 ((β_1+β_12μ_2)^2+2 β_12^2 ρ^2 σ_2^2/(ρ^2+1)^2)/2 ρσ_1 σ_2 (β_1+β_12μ_2) (β_2+β_12μ_1)+σ_1^2 (β_1+β_12μ_2)^2+σ_2^2 ((β_2+β_12μ_1)^2+β_12^2 (ρ^2+1) σ_1^2)
    𝕊^b_1    = ρσ_1 σ_2 ((ρ^2+1)^2 (β_1+β_12μ_2) (β_2+β_12μ_1)+2 β_12^2 ρ^3 σ_1 σ_2)/(ρ^2+1)^2 (2 ρσ_1 σ_2 (β_1+β_12μ_2) (β_2+β_12μ_1)+σ_1^2 (β_1+β_12μ_2)^2+σ_2^2 ((β_2+β_12μ_1)^2+β_12^2 (ρ^2+1) σ_1^2))
    𝕊^a_12   = β_12^2 (ρ^2-1)^2 σ_1^2 σ_2^2/(ρ^2+1) (2 ρσ_1 σ_2 (β_1+β_12μ_2) (β_2+β_12μ_1)+σ_1^2 (β_1+β_12μ_2)^2+σ_2^2 ((β_2+β_12μ_1)^2+β_12^2 (ρ^2+1) σ_1^2))
    𝕊^b_12   = 0,



where 
    𝕊_i = 𝕊^a_i + 𝕊^b_i,   i∈{1,2,12}
 and satisfy for general parameters 
    𝕊_1 + 𝕊_2 + 𝕊_12 = 1.
 Note that when ρ=0, then 𝕊^b_1=𝕊^b_2=𝕊^b_12=0. All non-trivial projections depend upon ρ and β. In other words, the projections are functions of the measure μ and of the model f. Note that 𝕊^b_12=0 for any ρ due to the hierarchical-orthogonality of f_12, e.g., for every ρ we have ⟨ f_12,f_1⟩=⟨ f_12,f_2⟩=0 but ⟨ f_1,f_2⟩=0 only for ρ=0.
 

Consider ρ→ 1. We know from the direct form of f(β,x) that in the limit we have x_1=x_2 and f(β,x)=β_0+(β_1+β_2)x + β_12x^2. That is, f degenerates from a two-dimensional function into a one-dimensional function as the coherence indicated by ρ achieves unity. HDMR uniquely captures this behavior in a variance-preserving manner for general distributions. We illustrate this property in Figure <ref> with plots of the sensitivity indices as functions of ρ for fixed μ, σ and β. Colored regions indicate increased variance (green) or reduced variance (red) relative to ρ=0. For example, Figure <ref> illustrates that the HDMR component function subspace in X_1× X_2 experiences annihilation for |ρ|→ 1.







In the next section we make a study of (<ref>) using other interpretative diagnostics in machine learning and compare to the results of this section.



§ INTERPRETATIVE DIAGNOSTICS
 Interpretative diagnostics provide explanation for prediction <cit.>. Many common interpretative diagnostics are dependent upon the choice of F  and are said to be class-dependent. For example, variable depth in decision nodes of ensembles of decision trees and node weight pruning in neural networks are commonly employed for supervised machine learning settings to attain measures of variable importance. Other interpretative diagnostics are defined in terms of an objective function J. For suitable J finite differences or derivative-based sensitivity analysis are conducted to assess the importance of various subsets of variables. We use the risk function R as the objective function J. We compare variable importance and dependence measures using I^u and P^u defined by partial dependence <cit.>, derivative-based global sensitivity indices <cit.>, and HDMR for the model (<ref>). We illustrate that HDMR preserves explained variance independent of correlation in μ, whereas partial dependence is valid for modest correlation (ρ<0.3). 



 §.§ Partial dependence
The partial dependence of f(x) on x_u can be defined a couple ways <cit.>. One is given by 
    P^u_PD f(x) ≡M^u f(x)=∫_X_-uf(x)μ_-u(x_-u) = f^PD_u(x_u),
 which is useful whenever the dependence between x_u and x_-u is not too strong. Another formulation of partial dependence is given by 
    P^u_PD f(x) ≡N^u f(x)=∫_X_-u f(x)μ_-u|u(x_-u) = f^PD_u(x_u),
 which considers the effect of dependencies. Here, the function is averaged with respect to the conditional distribution. If f̃(x)=f(x)-f_0 and the inputs are independent, then partial dependence is related to HDMR, where P^i_PDf̃(x)=P^i_HDMRf̃(x)=f̃_i^HDMR(x_i) for i∈ℕ_n. Because of this relationship, we define the importance functional for partial dependence in a manner similar to HDMR, 
    I^i_PD f(x) ≡(S^a PD_i, S^b PD_i, S_i^PD).
 



 §.§ Derivative-based global sensitivity measures
 <cit.> introduced derivative-based global sensitivity measures. Weighted derivative-based global sensitivity measures (DGSM) are defined as 
    I^u_DGSMf(x)≡∫_X (D^α_uf(x))^2w_u(x_u)μ(x) = S_u^DGSM,
 where 
    D^α_uf(x) = ∂^|α| f(x)/∂ x_1^α_u1…∂ x_n^α_un
 and w_u(x_u) is a weight function. Putting D^α f=(D^α_uf(x): u), a normalized DGSM is given as 
    I^u_DGSMf(x)≡∫_X (D^α_uf(x))^2w_u(x_u)μ(x)/D^α f_2;F^2 = S_u^DGSM,
 where 
    f_2;F=(∑_f∈ff_L^2(Ω,ℱ,μ)^2)^1/2.
 We take w_u(x_u)=1. 




 §.§ Illustration
 We illustrate these interpretative diagnostics for equation (<ref>), 
    f(β,x) = β_0 + β_1x_1 + β_2x_2 + β_12x_1x_2.
 We take β=1, μ=1, and σ=1, and examine variable importance as a function of correlation, ρ. For DGSM, we also consider σ_f^2(β) = β_2^2+2 β_2 β_3 ρ+β_3^2+β_12^2 (ρ^2+1). Whenever ρ∈{-1,1} such that x_1=x_2 the function is one-dimensional and all variable importance resides in the univariate importance measures. The component functions are analytically given by 

    f̃_1^HDMR(x_1)    = x_1 + ρ/ρ^2+1(x_1^2 - 1)
    f̃_1^PD (marg.)(x_1)    = x_1
    f̃_1^PD (cond.)(x_1)    = x_1 + ρ(x_1^2+x_1-1).

Notice that the HDMR component function has non-linear dependence on ρ, whereas PD exhibits linear dependence. Variable importance results are exhibited in Table <ref> and Figure <ref> (note that  the minimum and curvature for the `U'-shaped behavior changes in β and for β=1 in Figure <ref> it is not symmetric), and variable dependence results are exhibited in Figure <ref>. For HDMR and partial dependence, we breakdown the importance measure into structural and correlative terms. We compute importance measures for ρ∈{-1,0,1}. All importance measures, with the exception of I_DGSMf(x), depend upon ρ. When ρ=0, all coincide except I_DGSMf(x). Only HDMR preserves variance, S(-1)=S(1)=1. 


The variance decomposition property of HDMR, i.e. S(-1)=S(1)=1, means that the HDMR expansion correctly captures the contributions of subsets of variables to the variance in the output. Observe that partial dependence given by (<ref>) is accurate when correlation is roughly less than 0.3. As the strength of dependence grows, the univariate partial dependency interpretative diagnostics of (<ref>) become increasingly distorted. The diagnostic given by (<ref>) exhibits far less variance stability in one-dimensional subspaces than that given by (<ref>). PD is also estimated using a gradient boosting machine (PD (GBR), where GBR-estimated PD is estimated thirty times on a grid for ρ, each GBR having 500 estimators of max depth of four and 10^3 independent random samples. PD (GBR) sensitivity indices are computed from the GBR-estimated PD's using quadrature. GBR estimation of PD (<ref>) is significantly biased as correlation increases. DGSM gives different importance of f(x) than HDMR and partial dependence, even for ρ=0, although it coincides for σ_f^2(β).


Figure <ref> reveals the profile of f̃_1(x_1) in ρ for HDMR and partial dependence (marginal, conditional, and GBR-estimated). PD (<ref>) does not change in ρ, whereas PD (<ref>) does. PD (GBR) significantly deviates from both HDMR and PD (<ref>). When ρ=0, the four coincide.


These results illustrate that PD-based interpretative diagnostics experience information leakage whenever ρ>0. In a similar manner, GBR-estimated PD exhibits information leakage. It is interesting that GBR-estimated PD diverges from both analytic PD-measures. This may be attributed to biases in the underlying algorthmic implementation, as tree-based ensemble methods such as random forest are known to exhibit biases in variable importance <cit.>.











§ GLASS BOXES FROM BLACK BOXES
  General-purpose black box learning algorithms in machine learning, such as kernel machines or decision tree models, exhibit favorable predictive performances and configuration or tuning costs. Suppose the availability of a black box representation f∈ F.  We demonstrate how HDMR assesses variable importance and dependence of kernel machines or ensembles of decision trees. In this manner, HDMR is said to be a wrapper method <cit.>. We give HDMR constructions for kernel machines and ensembles of decision trees and provide illustrations for polynomial kernel machines and ensembles of gradient boosting machines. 



 §.§ Kernel machines
 HDMR provide interpretative diagnostics for the output of a kernel machine through a decomposition of the inner product of its RKHS. This is similar to the approach used in smoothing spline ANOVA RKHS but such is based on orthogonal projections <cit.>. We outline the approach below and discuss polynomial kernels. 


Assume we have a kernel 
    K(x,x')=⟨ϕ(x),ϕ(x')⟩_ℓ^2,
 where ϕ: X↦ℓ^2(α) is a feature map with index set α of size b∈ℕ and ϕ(x) = (ϕ_α(x): α∈α) is a feature vector. We assume these feature vectors are represented by the Hadamard (entrywise) product of constants and basis elements, ϕ(x) = β⊙B(x), where B(x)=(B_α(x) : α∈α) are bases and β∈ℝ^b are coefficients (note that B(x) is known as a Mercer map). The HDMR of the kernel is attained by decomposing the feature vector as 
    ϕ(x)=∑_α∈αν_αψ_α(x)=∑_α∈αν_α(η_α⊙B(x)),
 where {ψ_α(x)≡η_α⊙B(x)} are component feature vectors and {η_α = (η_αα' : α'∈α)} are coefficient vectors attained from μ. These coefficients reflect the hierarchical-orthogonality of the feature vectors and their generation is detailed in the appendices. Putting A≡(η_α)∈ℝ^b× b, we have 
    β=A ν
 where ν=(ν_α : α∈α)∈ℝ^b. Given i) β from the kernel and ii) A from the measure μ, this system is solved for ν, 
    ν= A^-1β.
 Using {ν_αψ_α(x)}, we form the HDMR of K by putting 
    ϕ(x) ≡∑_u⊆ℕ_nΨ_u(x),
 where the collection {Ψ_u(x)} is formed as 
    {Ψ_u(x):∀ u⊆ℕ_n,   α_u⊂α,   Ψ_u(x)≡∑_α∈α_uν_αψ_α(x)}.
 If a subset of component feature vectors is sought, {Ψ_u(x) : |u|≤ T<d}, then we have A∈ℝ^a× b, an underdetermined system a<b, and a least-norm solution can be efficiently attained using QR decomposition. Given a finite dataset D and a kernel K, a kernel machine attains a set coefficients {ξ_x'} such that

    f(x)   =∑_x'∈Dξ_x' K(x,x')
       =∑_x'∈Dξ_x'⟨ϕ(x),ϕ(x')⟩_ℓ^2
       =⟨∑_u⊆ℕ_nΨ_u(x),∑_x'∈Dξ_x'ϕ(x')⟩_ℓ^2
       =∑_u⊆ℕ_n⟨Ψ_u(x), f⟩_ℓ^2
    
    
    
       =∑_u⊆ℕ_nf^HDMR_u(x_u).



	


  §.§.§ Polynomial kernels
 The (real-valued) polynomial kernel is given by 
    K_poly(x,x';c,d,γ) = (γ x· x'+c)^d,
 where c,γ≥ 0, d∈ℕ, and x,x'∈ℝ^n. Defining α(n,d)={α∈ℕ_0^n+1, |α|= d}, 
    B(x;n,d) = (x^α_1…α_n : α∈α(n,d))
 and 
    β(n,d) = (β_α : α∈α(n,d)),
 
    β_α = (dαc^α_0γ^d-α_0)^1/2,
 we have (suppressing notation in n and d) 
    ϕ(x)=β⊙B(x).
 Together, 

    K_poly(x,x';c,d,γ)   =⟨ϕ(x),ϕ(x')⟩_ℓ^2
       =∑_α∈αβ_α^2 x^α_1…α_nx'^α_1…α_n.
  

Defining C(x;n,d)=(e^2π i⟨α_1… n,x⟩ : α∈α(n,d)) and ϕ(x)=β⊙C(x), we write a complex-valued polynomial kernel,  

    K_poly(x,x';c,d,γ)   =⟨ϕ(x),ϕ(x')⟩_ℓ^2
       =(γ e^2π ix· e^2π ix'+c)^d
       =∑_α∈αβ_α^2 e^2π i⟨α_1… n,x⟩e^2π i⟨α_1… n,x'⟩
       =∑_α∈αβ_α^2 e^2π i⟨α_1… n,x+x'⟩.
 




  §.§.§ Illustration: analytic function
 We consider the mathematical function (<ref>). In Figure <ref>, we compare f_1(x_1) as computed by kernel machines—f^HDMR_1(x_1), attained from the empirical HDMR (<ref>) (“empirical-polynomial HDMR”) and f^ANOVA_1(x_1), attained from the ANOVA representation of the polynomial kernel (<ref>) (“ANOVA-polynomial kernel”). The ANOVA-kernel of the polynomial kernel is defined as K_ANOVA(x,x') = ∏_i∈ℕ_n(1+k_i(x_i,x_i')) <cit.>, where { k_i(x_i,x_i')} are univariate zero-mean polynomial kernels. Both kernel machines use the same empirical measure of 10^3 data elements with ρ=1/2 and are compared to the exact HDMR. As seen in Figure <ref>, the ANOVA-polynomial kernel does not use the correlative information of the input data, whereas the empirical-polynomial HDMR better approximates the exact underlying HDMR.







 §.§ Decision trees
 A decision tree is a function that uses a partition of X=∪_R_i∈ RR_i, R={ R_i}_i, such that 
    g(x) = ∑_R_i∈ R c_i 1_R_i(x).
 The partition R_u is defined using a subset of variables x_u such that a tree is written as 
    g_u(x_u) ≡∑_R_i∈ R_u c_i 1_R_i(x_u).
 A collection of tree ensembles, indexed by ℐ, is denoted by 
    F_ℐ≡{ f^i(x)}_i∈ℐ,
  where f^i(x) is a sum of trees across the subspaces 𝒫_i⊆(ℕ_n) 
    f^i(x)≡∑_u∈𝒫_i g^i_u(x_u)=∑_u∈𝒫_i∑_j g^i_uj(x_u)
 and has depth d∈ℕ <cit.>. A key property of trained decision tree models is that tree subspaces are highly sparse, |𝒫_i|≪ |(ℕ_n)| for i∈ℕ_n. A linear combination of tree ensembles is defined as 
    f(x) ≡∑_i∈ℐβ_i f^i(x).
 Forming F_ℐ≡(f^i: i∈ℐ) = (g_u^i: i∈ℐ, u∈𝒫_i) we take 
    β(λ)≡(β_ui(λ): i∈ℐ, u∈𝒫_i)=(F_ℐ⊗F_ℐ+λ I)^-1(F_ℐ⊗ f).
 We index the tree ensemble feature vectors on α as

    α   ≡ (α_0)∥( (i,u)∈ℐ×(ℕ_n) : i∈ℐ,   u∈𝒫_i)

where 

    α_0    = ({},{})

and define  

    α_∙ u   ≡{(i,w)∈α : i∈ℐ,    w = u}
    α_i∙   ≡{(j,u)∈α : u∈𝒫_j,    j = i}.

The tree ensemble feature vectors are given by

    ϕ^i(x)   ≡β^i(λ)⊙B^i(x)
    B^i(x)    ≡(g^i_u(x_u)1_α_i∙(j,u) : (j,u)∈α)
    β^i(λ)   ≡ (β_i(λ)1_α_i∙(j,u) : (j,u)∈α).
 This is written as 
    ϕ(x)=∑_i∈ℐ(β^i(λ)⊙B^i(x))=β(λ)⊙B(x),
 where 
    B(x)=∑_i∈ℐB^i(x)
 and 
    β(λ)=∑_i∈ℐβ^i(λ).
 Note that β_0=0.
As before, we have the coefficients attained from μ as η_α and in organized form, 
    A≡ (η_α).
 Then, weights ν=(ν_α : α∈α) are attained that satisfy 
    β(λ)=A ν
 as 
    ν= A^-1β(λ).
 Noting

    ϕ(x)=∑_α∈αν_αψ_α(x)=∑_α∈αν_α(η_α⊙B(x)),
 we form the HDMR of F_ℐ as 
    {Ψ_u(x):∀ u⊆ℕ_n,   α_u⊂α,   Ψ_u(x)≡∑_α∈α_uν_αψ_α(x)}.


We summarize the computations: 



  * Equation (<ref>): forming F_ℐ from the collection of GBR models { f^i(x)}_i∈ℐ, indexed on ℐ

  * Equations (<ref>) and (<ref>): linearly combining the elements of F_ℐ as f(x) using the coefficients β(λ)≡(β_ui(λ): i∈ℐ, u∈𝒫_i)∈ℝ^b, attained from regularized least-squares using with grid-searched λ estimated using K-fold cross-validation

  * Equation (<ref>): attaining η_α from the measure μ using a collection of QR-decompositions, each identified to a component function subspace and each QR decomposition having linear cost in the dataset size 𝒪(N) and quadratic cost 𝒪(|𝒫(x_u)∩𝒫_i|^2) (typically small in most HDMR applications) and forming the square matrix A∈ℝ^b× b, where b is total size of the corresponding inner product space |∪_i𝒫_i|

  * Equations (<ref>) and (<ref>): forming the system β=A·ν, computing the inverse of the square matrix A, such as through using Gaussian elimination, and putting ν≡A^-1β

  * Equations (<ref>) and (<ref>): forming the HDMR F_ℐ.


Projections into the space of decision trees, which are systems of simple functions, produce noisy approximations to continuous functions. This can be addressed through spectral filtering of the component function subspaces, such as projection into smooth subspaces (Fourier). Another consideration is that if the GBR ensemble experience information leakage on the subspaces of F the interpretative diagnostics will be biased (as illustrated in the previous section for non-trivial correlation). The variance-preserving property of HDMR, however, enables introspection of black box learning algorithms for independent and/or correlated variables, to the extent that the black box contains information on the projections of the system. 


We compute the HDMR of two efficient general-purpose black-box supervised machine learning algorithms—the random forest (RF) and the gradient boosting regressor (GBR) machine—for a non-linear mathematical function (Ishigami) with analytic solution and a benchmark dataset (California housing). We demonstrate that RF experiences far more information leakage than GBR and that GBR well approximates HDMR.




  §.§.§ Illustration 1: analytic test function
 In this example we verify that the HDMR of F_ℐ given by (<ref>) well approximates the true HDMR for a test function. We consider the Ishigami function, a non-linear continuous function, which exists in closed-form and is sparse. It is defined as 
    f(x) = sin x_1 + a sin^2 x_2 + b  x_3^4sin x_1
 with independent uniformly distributed x=(x_1,x_2,x_3)∈[-π,π]^3. Its HDMR is analytic and shown below,

    f^HDMR_0    = a/2
    
    f^HDMR_1(x_1)    = (1+bπ^4/5)sin x1
    
    f^HDMR_2(x_2)    = -a/2cos 2x_2
    
    f^HDMR_3(x_3)    = 0
    
    f^HDMR_12(x_1,x_2)    = 0
    
    f^HDMR_13(x_1,x_3)    = b(x_3^4-π^4/5)sin x_1
    
    f^HDMR_23(x_2,x_3)    = 0
    
    f^HDMR_123(x_1,x_2,x_3)    = 0.

Its non-trivial sensitivity indices are

    𝕊_1    = 36 (π ^4 b+5)^2/5 (45 (a^2+4)+20 π ^8 b^2+72 π ^4 b)
    𝕊_2    = 45 a^2/45 (a^2+4)+20 π ^8 b^2+72 π ^4 b
    𝕊_13   = 64 π ^8 b^2/5 (45 (a^2+4)+20 π ^8 b^2+72 π ^4 b)




We utilize an ensemble of GBR's indexed on 

    ℐ_i   ={(depth_i,subspaces)}
       ={(k,(1,2,3)): k∈[i]}.

 These are linearly combined using regularized least-squares with grid-searched λ. J=10^3 trees are used in each ensemble and N=5× 10^3 data points are sampled from Uniform[-π,π]^3. 
 
In Figures <ref> and <ref> we estimate f^HDMR_1(x_1), f^HDMR_2(x_2), and f^HDMR_13(x_1,x_3) using GBR from F_ℐ_3. We also estimate these using RF (located in the appendices; Figures <ref>, <ref>, and <ref>). We project the one-dimensional component functions into the space of zero-mean band-limited Fourier series (degree four) to attain smooth low-frequency approximations. Sensitivity indices are shown below in Table <ref>. The GBR-based decision tree approximations closely approximate the true HDMR component functions, whereas RF performs poorly, experiencing significant information leakage in comparison to GBR. GBR estimation of the two-dimensional component function f^HDMR_13(x_1,x_3) experiences some information leakage but its general qualitative behavior is satisfactorily reproduced as shown in Figure <ref>. A tree-based measure of variable importance based on the fraction of samples variables contribute through trees is computed as `Tree'  and gives substantially different results than R by HDMR: x_3 is the most important variable by fraction of samples but the least important by HDMR. The defect in RF is appears to be insensitive to the number of trees, sample size, or tree depth. 















 
 








  §.§.§ Illustration 2: California housing dataset
 We consider the California housing dataset. This dataset has nine variables, eight predictors and a single response variable `median house value.' We standardize all variables to have zero mean and unit variance and utilize an ensemble of GBR's indexed on 
    ℐ={(i,(1,…,8)): i∈[3]},
 each having 5×10^3 trees, to estimate the HDMR per (<ref>). 





  
*Variable importance  We conduct a structural and correlative sensitivity analysis (SCSA) and show sensitivity indices with ϵ≥ 0.01 below in Table <ref>. Latitude and longitude each have strong negative correlative contributions to explained variance and participate in an modest-sized interaction. The relative sensitivity indices R (derived from the total T) are compared to the tree-based feature importance measure. HDMR places more emphasis on `MedInc', `Latitude', and `Longitude' than the tree-based indices. 


































































  
*One-dimensional variable dependence  

Figures <ref> and <ref> show HDMR and partial dependence for `MedInc', `Latitude', and 'Longitude.' HDMR and partial dependence are similar for `MedInc,' although HDMR generally is smoother. For 'Latitude' and 'Longitude' they are somewhat different: HDMR profiles are generally vertically shifted to more negative values.   







  
*Two-dimensional variable dependence
In Figure <ref> we show the HDMR component function of latitude and longitude. In comparison to partial dependence (Figure 10.17 of <cit.>), HDMR reveals large positive values in scattered eastern localities and reveals positive contributions in northern-most and eastern-most locations. 







§ CONCLUDING REMARKS
 HDMR provides structured information on the variable importance and dependence for square-integrable functions. As such, it provides information necessary and sufficient to interpret the model output in terms of the model input and constitutes a glass box representation. This research highlights key utilities of HDMR to supervised machine learning applications and introduces algorithmic constructions from black box models. One is that HDMR can identify sources of information leakage. This is demonstrated for goodness-of-fit settings for “big-data” and popular machine learning interpretative diagnostics and black boxes. Second, HDMR is a useful tool for managing the curse of dimensionality, as it often admits efficient reduced-order representations of dense high-dimensional systems. Third, whenever input variables are correlated, the HDMR component functions are functions of the input distribution parameters including those for correlation. Fourth, HDMR can be applied as a wrapper method for black boxes to provision glass boxes. Collectively, these results suggest HDMR to have broad utility within machine learning. 






§ ACKNOWLEDGEMENTS
 CB acknowledges support for this project from the National Institute of Dental and Craniofacial Research of the National Institutes of Health (NIH-NIDCR Grant No. K08-DE-021430). HR acknowledges support from the National Science Foundation. We thank Grzegorz A. Rempala for comments helpful to improving the article.



§ COMPUTATIONAL IMPLEMENTATIONS
 PD (GBR) is computationally implemented using the scikit learn (version 0.18) classes GradientBoostingRegressor and partial-dependence. 



tocsectionAppendices




§ APPENDICES






































 §.§ χ_n^2 decomposition
 Recall the PGOF statistic 
    χ_n^2= n∑_k∈ [k_n](μ̂_n{ k}-μ_n{ k})^2/μ_n{ k}
 with empirical probability μ̂_n defined as 
    μ̂_n{ k} = n^-1∑_j∈ℕ_n I(X_nj,k).
 Put (x_1,…,x_n)≡ (X_n1,…, X_nn). Observe that the statistic χ^2_n is a function of the random observation vector (x_1,…,x_n). In <cit.> χ_n^2 is decomposed into two uncorrelated components, 
    χ_n^2 = n^-1(S_n+U_n)-n,
 where

    S_n = ∑_iμ^-1_n{ x_i}
 
and

    U_n = ∑_i,j
    i jμ_n^-1{ x_i} I(x_i, x_j)

with 
    𝔼  S_n = n  k_n,     𝔼 U_n=n(n-1),     𝔼 χ_n^2 = k_n - 1.
 This gives 
    χ^2_n = n^-1( μ_n^-1{ X_n}+2(n-1)(k_n-1)).
 
This is a second-order HDMR, 
    χ^2_n(x_1,…,x_n) = f_0 + ∑_i f_i(x_i) + ∑_i<j f_ij(x_i,x_j),
 where 

    f_0    = k_n-1
    
    f_i(x_i)    =n^-1(μ_n^-1{ x_i}-k_n)
    
    f_ij(x_i,x_j)   =2  n^-1(μ_n^-1{ x_i}  I(x_i, x_j) -1),

with 
    f_i = n^-2  μ_n^-1{ X_n}
 and 
    f_ij= 4 n^-2(k_n-1).
 
Note that if 
    (k_nn)^-1 μ_ n^-1{ X_n}0,
 then (2k_n)^-1/2∑_i f_i(x_i)→ 0 in probability. This is satisfied for uniform and power law (α∈[0,1)) discrete random variables. Hence, for k_n→∞ as n→∞, the first-order terms converge to constants and the asymptotic influence of the second-order terms dominates, i.e. the distributional limit of χ_n^2 is determined by the second-order component functions. This is equivalently stated using the sensitivity indices, ∑_i<j𝕊_ij 1. Putting n / √(k_n)→λ, it turns out that χ_n^2Gaussian, χ_n^2Poisson, and χ_n^2Degenerate. Hence, the distributional limit of χ_n^2 is standard Gaussian for [label=(*)] 
  * n, k_n→∞ with λ=∞ or for 
  * n→∞ and k_n=k<∞ ; however, for n, k_n→∞ and λ <∞ this is not true. Consequently, an improved statistic for χ_n^2 is defined using HDMR, 
    χ_nHDMR^2(x_1,…,x_n)   ≡ f_0+∑_i<jf_ij(x_i,x_j) 
       ≡χ_n^2(x)-∑_if_i(x_i).
 



A precise characterization of the asymptotic behavior of χ_n^2 is stated below. Note that this result also holds for power-law distributions. 
 Assume μ_n{ x} = k_n^-1 for x∈ K_n and n=1,2,…, as well as 
    n/√(k_n)→λ.
 Then, 

    χ_n^2-k_n/√(2 k_n)
    0    when λ = 0
    √(2)/λ Z - λ/√(2),    Z∼Poisson(λ^2/2)    when λ∈(0,∞)
    
    N∼Gaussian(0,1)   when λ = ∞.





 §.§ Simple product function
 For the function f(x)=∏_i=1^n x_i, the partition of variance is given by 
    (f) = σ_f^2 = ∑_i σ_i^2 + ∑_i<jσ_ij^2 + … + σ_1… n^2,

where

    σ_u^2    = ∫_X_uf^2_u(x_u)ν_u(x_u)
    ν_u(x_u)   =∏_i∈ uν_i(x_i)

and

    σ_f^2     =μ^2n((1+ρ^2)^n-1)
    
    σ_i^2    =μ^2nρ^2
    
    σ_ij^2    =μ^2nρ^4
    
    ⋮   
    σ_i_1… i_k^2     =μ^2nρ^2k.


Figure <ref> shows the percent of explained variance of f^T(x_1,…,x_n) in (n,T) for ρ=1/2. The curves correspond to different values of T, where (1,…,n) is left-to-right. As exhibited, increasing T increases the explained variance.




Figure <ref> reveals that the necessary value of T to attain a given percent of explained variance depends upon ρ.







 §.§ Polynomials as linear combinations of monomials


Defining the monomial 
    x_u = ∏_i∈ ux_i
 with iid x, having common μ and σ^2, we consider a linear combination of monomials, 
    e_T(β,x) = ∑_u
    |u|≤ Tβ_u x_u.
  The variance of e_T is 
    (e_T) = 𝔼[(∑_u
    |u|≤ Tβ_u(x_u - μ^|u|))^2] = ∑_u
    |u|≤ Tβ_u^2 (x_u) + 2∑_u<v
    |u|,|v|≤ Tβ_uβ_v(x_u,x_v),
 where, putting ρ=σ/μ, 
    (x_u) = μ^2|u|((1+ρ^2)^|u|-1).
 
    (x_u,x_v) = μ^|u|+|v|((1+ρ^2)^|u∩ v| - 1).
 

e_T(β,x) is factorized as










    e_T(x)    = x_i(∑_u
    |u|≤ T-1
    |i∩ u|=0β_iu x_u) + (∑_u
    |u|≤ T
    |i∩ u|=0β_u x_u)
       = x_ij(∑_u
    |u|≤ T-2
    |ij∩ u|=0β_iju x_u) + x_i(∑_u
    |u|≤ T-1
    |ij∩ u|=0β_iu x_u) + x_j(∑_u
    |u|≤ T-1
    |ij∩ u|=0β_ju x_u) + (∑_u
    |u|≤ T
    |ij∩ u|=0β_u x_u),









For β=1, we call e_T(x) an elementary symmetric polynomial and can compute its conditional expectations as

    𝔼[e_T]    = a_00
    𝔼[e_T|ℱ_i]    = a_11x_i + a_10
    𝔼[e_T|ℱ_ij]    = a_22x_ij + a_21(x_i+x_j) + a_20.
    ⋮   
    𝔼[e_T|ℱ_u]    =∑_w⊆ ua_|u||w|x_w,

where  
    a_rs = ∑_u
    |u|≤ T-s
    |r∩ u|=0β_ruμ^|u| = ∑_k≤ T-sn-rkμ^k.
 



The component functions of e_T(x) are defined recursively,


    f_0    = 𝔼[e_T]
    
    f_i(x_i)    = 𝔼[e_T|ℱ_i]  - f_0
    
    f_ij(x_i,x_j)    = 𝔼[e_T|ℱ_ij] - f_i(x_i) - f_j(x_j) - f_0
    ⋮   
    
    f_u(x_u)    = 𝔼[e_T|ℱ_u] - ∑_w⊂ u f_w(x_w).


For T=n, we have 
 
    σ_f^2 = (μ +1)^2 n((τ^2+1)^n-1)

 
    f_u(x_u) = (1+μ)^n-|u|∏_i∈ u(x_i-μ)

 
    σ^2_u = τ^2|u| (μ +1)^2n
 
 
    q{ k} = ∑_i_1<…<i_k𝕊_i_1… i_k = nkτ^2k/(τ^2+1)^n-1
 where 
    τ = σ/μ+1= ρ/1+1/μ.
 When τ<1, q{ k}>p{ k} for k≪ n. For the uniform distribution on the unit interval, we have τ=1/3 √(3)≈ 0.19. Furthermore, 
    lim_μ,σ→∞
    σ/μ=ρτ=ρ.

 

 
  

 

 







 §.§ Feature vector construction for hierarchically-orthogonal subspaces


Suppose for the measure space (X_u,𝒳_u,μ_u) of subspace u⊆{1,…,n}, where μ_u(x_u)=∫_X^-uμ(x_u, x_-u), we have a collection of basis vectors indexed on (u) and having dimensions b=(b_v∈ℕ: v∈(u)). This is denoted by 
    B≡{B_v≡(B_v1,…,B_vb_v): v∈(u)}.
 We can construct a non-orthogonal basis vector Φ_u having dimension b_u whose elements are hierarchically-orthogonal with respect to μ (and μ_u). To do this, we order (u) in size and use the Gram-Schmidt process to generate Q=(Q_v≡(Q_v1,…,Q_vb_v): v∈(u)) from B. If the bases are given on an empirical measure as column vectors of data B∈ℝ^N×(b), then QR decomposition is conducted on B, returning Q. Then we put Φ_u=Q_u. Empirically, this is Φ_u∈ℝ^N× b_u. Repeating this process for all subspaces of interest α we form the feature vector 
    Φ=Φ_∅∥Φ_{1}∥…∥Φ_{ n}∥Φ_{1,2}∥…=_α∈αΦ_α.
 Putting 
    Φγ = f,
 the coefficients are given as 
    γ = (γ_α: α∈α)=(Φ⊗Φ)^-1(Φ⊗ f).
 These outer products can be calculated exactly, as they are computations of moments of μ. The component functions are given by 
    f_u(x_u) = ⟨γ_u, Φ_u(x_u)⟩_ℓ^2.
 















 §.§ Correlated expansion


A non-orthogonal HDMR basis {Φ_0,Φ_1,Φ_2,Φ_12} is formed as 

*!


    Φ_0    = (1)
    Φ_i    = (x_i-μ_i/σ_i,(x_i-μ_i)^2/σ_i^2-√(2)),     i=1,2
    Φ_12   = (-ρσ_2^2 ((ρ^2-1) σ_1^2+(x_1-μ_1)^2)+(ρ^2+1) σ_1 σ_2 (x_1-μ_1) (x_2-μ_2)+ρσ_1^2 (-(x_2-μ_2)^2)/(ρ^2+1) √(ρ^2+4/ρ^2+1-3)σ_1^2 σ_2^2).



The coefficients are given by

    γ_0    = (β_0+β_1μ_1+β_2μ_2+β_12μ_1 μ_2+β_12ρσ_1 σ_2)
    γ_1    = (σ_1 (β_1+β_12μ_2),√(2)β_12ρσ_1 σ_2/ρ^2+1) 
    γ_2    = (σ_2 (β_2+β_12μ_1),√(2)β_12ρσ_1 σ_2/ρ^2+1) 
    γ_12   = (β_12√(ρ^2+4/ρ^2+1-3)σ_1 σ_2)























































 §.§ Analytic test function


 

 

 





 §.§ California housing dataset
 






























0.2in
plainnat
25


urlstyle

[Arwade et al.(2010)Arwade, Moradi, and Louhghalam]Arwade20101
Sanjay R. Arwade, Mohammadreza Moradi, and Arghavan Louhghalam.
Variance decomposition and global sensitivity for structural systems.
Engineering Structures, 320 (1):0 1 – 10,
  2010.
ISSN 0141-0296.
http://dx.doi.org/10.1016/j.engstruct.2009.08.011.
URL
  <http://www.sciencedirect.com/science/article/pii/S0141029609002788>.

[Breiman(2001)]Breiman2001
Leo Breiman.
Random forests.
Machine Learning, 450 (1):0 5–32, 2001.
ISSN 1573-0565.
10.1023/A:1010933404324.
URL <http://dx.doi.org/10.1023/A:1010933404324>.

[Chastaing et al.(2012)Chastaing, Gamboa, and Prieur]hoeff
Gaelle Chastaing, Fabrice Gamboa, and Clémentine Prieur.
Generalized hoeffding-sobol decomposition for dependent variables -
  application to sensitivity analysis.
Electron. J. Statist., 6:0 2420–2448, 2012.
10.1214/12-EJS749.
URL <http://dx.doi.org/10.1214/12-EJS749>.

[Durrande et al.(2013)Durrande, Ginsbourger, Roustant, and
  Carraro]kernel
N. Durrande, D. Ginsbourger, O. Roustant, and L. Carraro.
{ANOVA} kernels and {RKHS} of zero mean functions for model-based
  sensitivity analysis.
Journal of Multivariate Analysis, 115:0 57 – 67,
  2013.
ISSN 0047-259X.
http://dx.doi.org/10.1016/j.jmva.2012.08.016.
URL
  <http://www.sciencedirect.com/science/article/pii/S0047259X1200214X>.

[Fisher(1921)]fisher
Sir Fisher, Ronald Aylmer.
On the "probable error" of a coefficient of correlation deduced from
  a small sample.
Metron, 10 (4):0 3–32, 1921.

[Friedman(2001)]friedman2001
Jerome H. Friedman.
Greedy function approximation: A gradient boosting machine.
Ann. Statist., 290 (5):0 1189–1232, 10 2001.
10.1214/aos/1013203451.
URL <http://dx.doi.org/10.1214/aos/1013203451>.

[Guyon and Elisseeff(2003)]guyon
Isabelle Guyon and André Elisseeff.
An introduction to variable and feature selection.
J. Mach. Learn. Res., 3:0 1157–1182, March 2003.
ISSN 1532-4435.
URL <http://dl.acm.org/citation.cfm?id=944919.944968>.

[Hajek(1968)]hajek1968
Jaroslav Hajek.
Asymptotic normality of simple linear rank statistics under
  alternatives.
Ann. Math. Statist., 390 (2):0 325–346, 04
  1968.
10.1214/aoms/1177698394.
URL <http://dx.doi.org/10.1214/aoms/1177698394>.

[Hastie et al.(2009)Hastie, Tibshirani, and Friedman]esl
Trevor J. Hastie, Robert John Tibshirani, and Jerome H. Friedman.
The elements of statistical learning : data mining, inference,
  and prediction.
Springer series in statistics. Springer, New York, 2009.
ISBN 978-0-387-84857-0.
URL <http://opac.inria.fr/record=b1127878>.
Autres impressions : 2011 (corr.), 2013 (7e corr.).

[Hoeffding(1948)]hoeffding
Wassily Hoeffding.
A class of statistics with asymptotically normal distribution.
Ann. Math. Statist., 190 (3):0 293–325, 09
  1948.
10.1214/aoms/1177730196.
URL <http://dx.doi.org/10.1214/aoms/1177730196>.

[Hooker(2007)]hooker
Giles Hooker.
Generalized functional anova diagnostics for high-dimensional
  functions of dependent variables.
Journal of Computational and Graphical Statistics, 160
  (3):0 709–732, 2007.
10.1198/106186007X237892.
URL <http://dx.doi.org/10.1198/106186007X237892>.

[I.M. Sobol(2004)]Sobol2004
S.S. Kucherenko I.M. Sobol.
On global sensitivity analysis of quasi-monte carlo algorithms.
Monte Carlo Methods and Applications, 11:0 83–92,
  2004.

[Li and Rabitz(2012)]li12
Genyuan Li and Herschel Rabitz.
General formulation of hdmr component functions with independent and
  correlated variables.
Journal of Mathematical Chemistry, 500 (1):0
  99–130, 2012.
ISSN 0259-9791.
10.1007/s10910-011-9898-0.
URL <http://dx.doi.org/10.1007/s10910-011-9898-0>.

[Li and Rabitz(2014)]Li2014
Genyuan Li and Herschel Rabitz.
Analytical hdmr formulas for functions expressed as quadratic
  polynomials with a multivariate normal distribution.
Journal of Mathematical Chemistry, 520 (8):0
  2052–2073, 2014.
ISSN 1572-8897.
10.1007/s10910-014-0365-6.
URL <http://dx.doi.org/10.1007/s10910-014-0365-6>.

[Lorentz et al.(1996)Lorentz, Golitschek, and Makovoz]kol
G.G. Lorentz, M.V. Golitschek, and Y. Makovoz.
Constructive Approximation.
Spring, New York, 1996.

[Rabitz and Alis(1999)]alis99
Herschel Rabitz and Omer F. Alis.
General foundations of high?dimensional model representations.
Journal of Mathematical Chemistry, 250 (2-3):0
  197–233, 1999.
ISSN 0259-9791.
10.1023/A:1019188517934.
URL <http://dx.doi.org/10.1023/A3A1019188517934>.

[Rempała and Wesołowski(2016)]rempala
Grzegorz A. Rempała and Jacek Wesołowski.
Double asymptotics for the chi-square statistic.
Statistics & Probability Letters, 119:0 317–325, 12
  2016.
http://dx.doi.org/10.1016/j.spl.2016.09.004.
URL
  <http://www.sciencedirect.com/science/article/pii/S0167715216301699>.

[Sobol(1990)]Sobol93
I. M. Sobol.
On sensitivity estimation for nonlinear mathematical models.
Matem. Mod., 20 (1):0 112–118, 1990.

[Sobol' and Kucherenko(2009)]dgsm1
I. M. Sobol' and S. Kucherenko.
Derivative based global sensitivity measures and their link with
  global sensitivity indices.
Mathematics and Computers in Simulation, 790
  (10):0 3009–3017, 6 2009.
http://dx.doi.org/10.1016/j.matcom.2009.01.023.
URL
  <http://www.sciencedirect.com/science/article/pii/S0378475409000354>.

[Sobol(2001)]Sobol01
I.M. Sobol.
Global sensitivity indices for nonlinear mathematical models and
  their monte carlo estimates.
Mathematics and Computers in Simulation, 550
  (1-3):0 271–280, 2001.

[Strobl et al.(2007)Strobl, Boulesteix, Zeileis, and
  Hothorn]Strobl2007
Carolin Strobl, Anne-Laure Boulesteix, Achim Zeileis, and Torsten Hothorn.
Bias in random forest variable importance measures: Illustrations,
  sources and a solution.
BMC Bioinformatics, 80 (1):0 25, 2007.
ISSN 1471-2105.
10.1186/1471-2105-8-25.
URL <http://dx.doi.org/10.1186/1471-2105-8-25>.

[Takemura(1983)]takemura
Akimichi Takemura.
Tensor analysis of ANOVA decomposition.
780 (384):0 894–900, ???? 1983.
ISSN 0162-1459 (print), 1537-274X (electronic).

[Vapnik(1995)]erm
Vladimir N. Vapnik.
The Nature of Statistical Learning Theory.
Springer-Verlag New York, Inc., New York, NY, USA, 1995.
ISBN 0-387-94559-8.

[Wahba(1990)]ss1
G. Wahba.
Spline Models for Observational Data.
Society for Industrial and Applied Mathematics, 1990.
10.1137/1.9781611970128.
URL <http://epubs.siam.org/doi/abs/10.1137/1.9781611970128>.

[Wahba et al.(1995)Wahba, Wang, Gu, Klein, and Klein]ss2
Grace Wahba, Yuedong Wang, Chong Gu, Ronald Klein, and Barbara Klein.
Smoothing spline anova for exponential families, with application to
  the wisconsin epidemiological study of diabetic retinopathy.
The Annals of Statistics, 230 (6):0
  1865–1895, 1995.
URL <http://www.jstor.org/stable/2242776>.




