Interpretable Architecture Neural Networks for Function Visualization
In many scientific research fields, understanding and visualizing a black-box
function in terms of the effects of all the input variables is of great
importance. Existing visualization tools do not allow one to visualize the
effects of all the input variables simultaneously. Although one can select one
or two of the input variables to visualize via a 2D or 3D plot while holding
other variables fixed, this presents an oversimplified and incomplete picture
of the model. To overcome this shortcoming, we present a new visualization
approach using an interpretable architecture neural network (IANN) to visualize
the effects of all the input variables directly and simultaneously. We propose
two interpretable structures, each of which can be conveniently represented by
a specific IANN, and we discuss a number of possible extensions. We also
provide a Python package to implement our proposed method. The supplemental
materials are available online.
Shengtong Zhang, Daniel W. Apley
# Interpretable Architecture Neural Networks for Function Visualization

###### Abstract

In many scientific research fields, understanding and visualizing a black-box function in terms of the effects of all the input variables is of great importance. Existing visualization tools do not allow one to visualize the effects of all the input variables simultaneously. Although one can select one or two of the input variables to visualize via a 2D or 3D plot while holding other variables fixed, this presents an oversimplified and incomplete picture of the model. To overcome this shortcoming, we present a new visualization approach using an interpretable architecture neural network (IANN) to visualize the effects of all the input variables directly and simultaneously. We propose two interpretable structures, each of which can be conveniently represented by a specific IANN, and we discuss a number of possible extensions. We also provide a Python package to implement our proposed method. The supplemental materials are available online.

_Keywords:_ Function Visualization; Interpretable Machine Learning; Neural Network
Introduction

Scientists and engineers often want to visualize the output of a black box function in order to understand and interpret the effects of input variables on the output (aka response) variable. For example, most surrogate models of complex computer simulations are notorious for their black box property and lack interpretability. Visualizing a black box model helps to understand how the inputs affect the predicted response variable, whether the predictive relationships seem intuitively reasonable, and how decisions can be made based on the model.

We denote the black box function by \(f(\mathbf{x})\), where \(\mathbf{x}=(x_{1},x_{2},\cdots,x_{d})\) is a \(d\)-dimensional vector of input variables, and \(f(\mathbf{x})\) is the scalar response. \(f(\mathbf{x})\), as a function of \(\mathbf{x}\), is often referred to as a response surface. In reality, many functions estimated from a large data set or a complex computer simulation model have many input variables, which makes visualization of the joint effects of all the input variables challenging.

Of the existing works that aim to interpret black box models, perhaps the most common and the most closely related works are "partial dependence" (PD) plots (Friedman, 2001) and "individual conditional expectation" (ICE) plots (Goldstein et al., 2015). Each PD plot shows the marginal effect of one or two selected input variables by integrating the response variable over the marginal distribution of the omitted input variables, whereas each ICE plot displays a collection of curves that are functions of a single selected input variable of interest, one curve for each fixed combination of the \(d-1\) omitted input variables. Accumulated local effects (ALE) plots (Apley and Zhu, 2020) improve upon PD plots by offering much faster computations and more accurate results when the input variables are highly correlated. Closely related to ICE plots, trellis plots (Becker et al., 1996, 1994) are a series of plots of the response variable as a function of a pair of selected input variables with the omitted variable(s) held fixed, with a separate plot for each fixed combination of the omitted variables. Although trellis plots typically present a clear and fairly complete picture
of \(f(\mathbf{x})\) for the case of \(d=3\), for \(d>3\) they become cumbersome, since there are too many fixed combinations of the \(d-2\) omitted variables to consider. Whereas PD, ICE, trellis and ALE plots focus on visualizing the effects of one or two variables with each plot, and do not present a clear picture of the interactions between the selected and omitted input variables for large \(d\), our approach aims to visualize the joint effects of all input variables. Moreover, Agarwal et al. (2021) combined deep neural networks with generalized additive models to increase interpretablity. The primary distinction is that to enable easy visualization, the neural network architecture in (Agarwal et al., 2021) is restricted to additive functional relationships of the form \(f(\mathbf{x})=f_{1}(x_{1})+f_{2}(x_{2})+\cdots+f_{d}(x_{d})\) and, therefore, cannot be used to represent and visualize interactions between the input variables.

Other less closely related works focus on calculating a variable importance measure for each input variable. Breiman (2001) introduced the idea of permutation-based feature importance in random forest models, and it was later extended to general black box models (Fisher et al., 2019; Konig et al., 2021). Based on game theory concepts, Lundberg and Lee (2017) used Shapley values to compute feature importance. However, the feature importance measures only provide a scalar numerical score of the importance of each input variable, without revealing how the variables affect the response variable.

Instead of singling out the effects of the original variables/features, some approaches aim to visualize the topological structure of \(f(\mathbf{x})\). These methods largely focus on identifying the number and locations of local minima and maxima of \(f\) and, subsequently, on identifying paths in the input space to traverse between the minima and maxima. Gerber et al. (2010) used the Morse-Smale complex to segment the \(d\)-dimensional input variable space into several non-overlapping regions over which \(f\) is monotonic and located local minima and maxima of \(f\). To visualize the topology of \(f\) over each segmented region, they constructed certain regression surfaces and embedded them in 2D as a simplified representation of \(f\) in each region. This approach was applied to nuclear reactor safety analysis
and visualization (Maljovec et al., 2013). Harvey et al. (2012) further used the Reeb graph to shatter the loops created by the Morse-Smale complex and provide a topologically simpler visualization. Although the local effects of input variables can be interpreted from the embedded regression curve for each segment in the Morse-Smale complex, the global effects are difficult to interpret, especially as the number of segments increases. Moreover, the method is not well suited for visualizing global interactions between a set of inputs.

Our visualization approach is based on the following approximate representation of \(f\):

\[f(x_{1},x_{2},\cdots,x_{d})\approx g(x_{j},h(\boldsymbol{x}_{\backslash j})) \quad\text{for some }j\in\{1,2,\cdots,d\}. \tag{1}\]

The notation "\(\boldsymbol{x}_{\backslash j}\)" represents all the input variables excluding \(x_{j}\). For simplicity, we focus on functions \(f\) defined on the Cartesian product \(\prod\limits_{j=1}^{d}I_{j}\), where each \(I_{j}\in\mathbb{R}\) is a closed interval. In Section 2 we show that for each \(j\in\{1,2,\cdots,d\}\) any continuous \(f\) on \([0,1]^{d}\) can be arbitrarily closely approximated by the structure (1) for some continuous functions \(g\) and \(h\) (although for some \(j\) the resulting functions \(g\) and \(h\) will be better behaved than for other \(j\)). We also describe an approach for estimating \(g\) and \(h\) and for selecting the index \(j\) that leads to a good approximation.

For visualizing \(f\) via (1) with a particular input variable \(x_{j}\) singled out, one can construct a 3D plot of \(f\) vs \(x_{j}\) and \(h\). To illustrate, consider the harmonic wave function from physics

\[f(\boldsymbol{x})=x_{1}*\sin[(2\pi/x_{2})*x_{3}+x_{4}]. \tag{2}\]

Here \(f\) is the displacement of a given point on the wave, \(x_{1}\in[0.5,2]\) is the amplitude, \(x_{2}\in[0.5,2]\) is the wavelength, \(x_{3}\in[0,1]\) is the position of that point (e.g., the distance from the source of the wave), and \(x_{4}\in[0,\pi]\) is the phase of the wave. \(f\) can be represented as (1) with \(x_{1}\) singled out, \(f(\boldsymbol{x})=x_{1}*\sin[h(\boldsymbol{x}_{\backslash 1})]\), and \(h(\boldsymbol{x}_{\backslash 1})=(2\pi/x_{2})*x_{3}+x_{4}\). The corresponding 3D visualization plot (\(f\) vs \(x_{1}\) and \(h\)) is shown in Figure 1(a). Such a plot allows one to visualize the manner in which \(f\) depends on the single variable \(x_{j}\), including how \(x_{j}\) interacts with some latent function (\(h\)) that captures the full extent of the \(x_{j}\) interactions with the other \(d-1\) input variables.
In this paper, we visualize the estimated functions mainly using static 3D plots. However, considering that static 3D plots may be easily misinterpreted due to visual misperception (Cleveland and McGill, 1984), we recommend that users plot the functions using whatever visual rendering they prefer, e.g., using software that allows one to interactively rotate the 3D plots to better avoid misinterpreting or missing features of the plot, supplementing the 3D plot with a 2D heat map with contour lines added, and/or using advanced shading and perspectives. In Figure 1, we also include a 2D heatmap of the same function and a shaded perspective plot using the software _rayshader_(Morgan-Wall, 2022).

In order to understand how \(x_{j}\) interacts with the other individual input variables and, more generally, to understand the joint effect of all \(d\) input variables on \(f\), our approach proceeds hierarchically and approximates \(h\) with a structure analogous to (1) with a second input variable singled out, and so on. We describe this approach in Section 3.1 and refer to it as the "original variable hierarchical" (OVH) structure. Moreover, we also consider a more general structure that assumes the hierarchical \(h\) functions are functions of certain linear combinations of the original inputs that can further enhance visualization. We refer to this as the "disjoint active subspace hierarchical" (DASH) structure and describe it in Section 3.2.

Figure 1: An example of a visualization per structure (1) of the harmonic wave function \(f\) with (a) a 3D plot, (b) a 2D heatmap of the same function and (c) a shaded perspective plot. See the online version for color figures.
The remainder of the paper is organized as follows. In Section 2 we present a theorem to show that any continuous function \(f\) can be arbitrarily closely approximated by a model with the more interpretable structure (1), and we describe our _interpretable architecture neural network (IANN)_ model to estimate the functions \(g\) and \(h\) in (1). In Section 3, we develop the OVH and DASH structures to visualize and interpret the joint effects of all the variables by hierarchically decomposing the latent function \(h\), both of which can be conveniently represented with a specific IANN model. We note that for visualizing complex black box simulation response surfaces that are expensive to evaluate, one should first fit a surrogate model (Sacks et al., 1989; Kennedy and O'Hagan, 2000, 2001) to the simulation data and then use the surrogate model as \(f\). In Section 4 we present algorithms for finding the appropriate ordering of the input variables in the hierarchical decompositions. In Section 5, we provide additional numerical examples to illustrate function visualization using our IANN approach. In Section 6, we discuss a number of potential extensions of the IANN approach.

## 2 IANN Structure for the First Hierarchical Level

This section describes how to use the IANN model to approximate \(f\) by structure (1) and how it facilitates the visualization of \(f\). In Section 2.1, we present Theorem 2.1 to show that the proposed structure (1) can approximate any continuous function on \([0,1]^{d}\) for some continuous functions \(g\) and \(h\). To estimate \(g\) and \(h\), we introduce the IANN architecture to approximate \(f\) for each singled out variable \(x_{j}\) in Section 2.2. In Section 2.3, we compare and draw a connection between PD plots, ICE plots, and the IANN visualization plots to show that the latter has advantages that lead to clearer interpretation of \(f\). This decomposition can be used as a stand-alone approach to visualize the effects of each input (as illustrated in Figure 1), if it is repeated in a nonhierarchical manner for each input \(x_{j}\) for \(j=1,2,\cdots,d\). The main information that is missing in this visualization is the
disentanglement of the joint effects of all inputs, which is the subject of Section 3. For this, the decomposition in this section provides the basic building block.

### IANN approximation theorem

The following theorem (see Appendix D for a proof) guarantees that the interpretable structure we proposed in (1) can approximate any continuous function \(f\) on \([0,1]^{d}\).

**Theorem 2.1**.: _Let \(f:[0,1]^{d}\to\mathbb{R}\) be a continuous function. For any \(\epsilon>0\) and \(j\in\{1,2,\cdots,d\}\), there exist continuous functions \(g\) and \(h\) satisfying:_

\[\big{|}f(\mathbf{x})-g(x_{j},h(\mathbf{x}_{\setminus j}))\big{|}<\epsilon,\quad\forall \mathbf{x}=(x_{j},\mathbf{x}_{\setminus j})\in[0,1]^{d}. \tag{3}\]

More generally, the theorem above holds if \(f\) is defined on the Cartesian product of intervals in each \(x_{j}\).

The significance of the theorem is that the structure (1) can be used to approximate arbitrarily closely a continuous \(f\), and this approximation can then be used to visualize \(f\) as in Figure 1. This visualization is related to ICE plots but provides more clarity in the sense that we discuss in Section 2.3. Moreover, since \(h\) is continuous, one can apply Theorem 2.1 again to show that \(h\) can be closely approximated by the same structure (1), and this process can be repeated to yield the hierarchical decomposition described in Section 3.

Notice that Theorem 2.1 applies for each \(j\in\{1,2,\cdots,d\}\), although there is a caveat that should be pointed out. For some \(j\), the corresponding \(h\) may be so complex that it cannot realistically be estimated and visualized in the subsequent levels. However, for all the real functions \(f\) we have considered, for at least one of the input variables \(x_{j}\) we were able to accurately approximate \(f\) with an \(h\) that was well behaved and could subsequently be approximated via the same structure (1) in the next level. Moreover, our algorithm (in Section 4) automatically finds a good hierarchical ordering of the input variables to yield an
accurate approximation of \(f\), and the accuracy of the approximation can be easily quantified as a check to verify whether the function \(f\) in question can indeed be approximated by an IANN structure.

### IANN architecture

To estimate the functions \(g\) and \(h\), we use the customized neural network architecture depicted in Figure 2, which we refer to as an IANN.

The input layer has \(d\) nodes that represent the \(d\) input variables, and the single node in the output layer is the approximation of the original function \(f(\mathbf{x})\). Due to Theorem 2.1 and the universal approximation theorem of neural networks (Hornik et al., 1989; Cybenko, 1989; Hornik, 1991), with enough layers and nodes the IANN in Figure 2 can approximate any continuous function \(f\) by the structure (1) for some functions \(g\) and \(h\) estimated in the process of fitting the IANN model.

The bottle-neck layer in the middle has only two nodes, representing \(x_{j}\) and \(\hat{h}\). \(x_{j}\) is

Figure 2: IANN architecture representing Eq. (1) at the top hierarchical level. The bottleneck layer in the middle consists of only two nodes, which represents \(x_{j}\) and \(\hat{h}(\mathbf{x}_{\setminus j})\).
directly connected to the first input node by an identity map, and \(\hat{h}\) is the estimation of the latent function \(h(\mathbf{x}_{\backslash j})\) in Eq. (1). In the layers to the left of the bottle-neck layer, we connect the \(d-1\) nodes (which represent the remaining input variables \(\mathbf{x}_{\backslash j}\)) to node \(\hat{h}\) by a fully-connected neural network denoted by \(\hat{h}(\mathbf{x}_{\backslash j})\). In the layers to the right of the bottle-neck layer, we use another fully-connected network (denoted by \(\hat{g}\), with inputs \(x_{j}\) and \(\hat{h}\)) to represent the function \(g\) as an approximation of \(f(\mathbf{x})\).

Denote the training data input/response observations by \(\{(x_{1}^{i},x_{2}^{i},\cdots,x_{d}^{i},f^{i}),\ 1\leq i\leq N\}\) and the IANN response prediction by \(\hat{f}^{i}=\hat{g}(x_{j}^{i},\hat{h}(\mathbf{x}_{\backslash j}^{i}))\), where \(N\) is the number of observations in the training set. The training data were generated using the customized Latin Hypercube sampling described in the Appendix B. Note that even when \(f\) represents a simulation response surface that is expensive to evaluate, \(N\) can be chosen quite large, because \(f\) is first replaced by a surrogate model fit to the simulation data, and generating response observations from the surrogate model is inexpensive. To fit the IANN (i.e., estimate the weights and biases) to the training data, we use standard squared error loss:

\[Loss=\frac{1}{N}\sum_{i=1}^{N}\left|f^{i}-\hat{f}^{i}\right|^{2}. \tag{4}\]

After fitting the IANN, the resulting approximation of \(f\) can be visualized by plotting \(\hat{g}\) as a function of \(x_{j}\) and \(\hat{h}\), as in Figure 1. The fitted IANN serves as an approximation of \(f\). To assess whether the approximation is adequate, we generate a large random sample in the input space to serve as a test set and compute the test \(r^{2}\) by comparing the IANN test predictions \(\hat{f}\) with the actual function values \(f\) at the test inputs.

Returning to the harmonic wave function defined in Eq. (2), we (nonhierarchically) single out each of the four input variables one-by-one to serve as \(x_{j}\) and fit the IANN architecture in Figure 2. The resulting four 3D plots (each of \(g(x_{j},h)\) as a function of \(x_{j}\) and \(h\) for the four different \(x_{j}\)) are displayed in Figure 3. The top left visualization plot has the highest test \(r^{2}\) (99.8%), which suggests that the harmonic wave function can be well approximated by the structure (1) after singling out the amplitude variable \(x_{1}\)
This is consistent with the structure of \(f\) in (2) (which in practice would typically be a more complex and less directly interpretable function, e.g., a surrogate model fit to some complex computer simulation output), since (2) can be written as \(f(\boldsymbol{x})=x_{1}*\sin\left[h(\boldsymbol{x}_{\setminus 1})\right]\) with \(h(\boldsymbol{x}_{\setminus 1})=(2\pi/x_{2})*x_{3}+x_{4}\). From the top-left plot it is clear that the effect of \(x_{1}\) is linear for any fixed \(\boldsymbol{x}_{\setminus 1}\), although the effect of \(x_{1}\) depends strongly on the value of \(\boldsymbol{x}_{\setminus 1}\). For example, for some \(h(\boldsymbol{x}_{\setminus 1})\) the effect of \(x_{1}\) is positive (positive slope for the linear trend), and for other values the effect is negative. Moreover, the slope varies periodically as a sinusoidal function of \(h(\boldsymbol{x}_{\setminus 1})\). Evidently, the IANN representation of \(h(\boldsymbol{x}_{\setminus 1})\) captures the actual function \((2\pi/x_{2})*x_{3}+x_{4}\) quite well, which is more clear from the hierarchical IANN visualization presented later in Section 5.2. The top-left plot in Figure 3 also illustrates how the approach can be used to understand interactions between inputs, as \(x_{1}\) and \(h(\boldsymbol{x}_{\setminus 1})\) clearly have a strong interaction, by which the effect of \(x_{1}\) changes from positive to negative and vice-versa as \(h(\boldsymbol{x}_{\setminus 1})\) varies.

For the other three plots, the relatively low test \(r^{2}\) values suggest that singling out those variables in the IANN cannot approximate the structure (1) well. Therefore, those plots are less reliable and should be used with caution when interpreting the effects of inputs. In this situation, we recommend using only the top-left plot of \(g(x_{1},h(\boldsymbol{x}_{\setminus 1}))\) (with test \(r^{2}=99.8\%\)) at the top level and then using the hierarchical decomposition in Section 3 to visualize the effects of \(\{x_{2},x_{3},x_{4}\}\). We do this in Section 5.2 in a continuation of this example. As pointed out after Theorem 2.1, for a particular \(x_{j}\) the latent function \(h\) may be too complex to estimate via a neural network with limited number of layers and nodes (and also too complex to easily visualize in subsequent levels), which would be revealed by a low associated test \(r^{2}\). But for all real examples that we have considered, the 3D visualization plot of \(g(x_{j},h(\boldsymbol{x}_{\setminus j}))\) for at least one \(x_{j}\) has sufficiently high test \(r^{2}\) to render the visualization reliable. If multiple plots have high test \(r^{2}\), users may rely on any or all of them to interpret the effects of inputs.
### A connection to PD and ICE plots

The most popular method for visualizing the effects of the input variables is partial dependence (PD) plots, which were first introduced in Friedman (2001). Later, Goldstein et al. (2015) proposed ICE plots, which enhance PD plots by displaying a collection of curves for each fixed combination of the omitted input variables. Here we compare and draw a connection between IANN, PD, and ICE plots using the harmonic wave function to illustrate.

Figure 4(a) shows both the PD plot and the ICE plot for the variable \(x_{1}\) using the open-source Python package, _scikit-learn_(Pedregosa et al., 2011; Buitinck et al., 2013). Notice that the PD plot, represented by the dashed line, is the average of all the curves in the ICE plots. There is no tight, direct connection between our IANN plots and ALE plots. Their only connection is via the connection between IANN plots and PD plots, and recognizing that ALE plots are intended to produce something similar to PD plots, albeit in a manner that is much more computationally efficient and less prone to problems when the inputs are highly correlated. Since the inputs have little correlation in the harmonic wave example, the ALE plot (omitted for brevity) is very similar to the PD plot in Figure 4(a). The ICE plots for this example are all straight lines with different slopes, which correctly suggests that the effect of \(x_{1}\) is linear for fixed \(\mathbf{x}_{\backslash 1}\) and that \(x_{1}\) interacts with the

Figure 3: Top-level (non-hierarchical) visualization plots of the harmonic wave function in Eq. (2) with each input singled out in turn. The test \(r^{2}\) is shown above each plot to evaluate the accuracy of IANN approximation represented by the plot.
other variables (since the slopes vary). In comparison, Figure 4(b) shows the IANN plot of \(g(x_{1},h(\mathbf{x}_{\backslash 1}))\) for the same example. Note that each individual curve in the ICE plot represents \(f\) as a function of \(x_{1}\) for a fixed \(\mathbf{x}_{\backslash 1}\). The IANN plot can be viewed as smoothly piecing together all of the ICE plot curves to create a 3D surface that serves as a more structured and interpretable visualization of the effect of \(x_{1}\) on \(f\) and of how \(x_{1}\) interacts with some function (\(h\)) of the remaining variables \(\mathbf{x}_{\backslash 1}\). Three individual curves for three fixed \(h(\mathbf{x}_{\backslash 1})\) values are shown as dashed curves in Figure 4(b), and the three corresponding ICE plot curves for the same three \(\mathbf{x}_{\backslash 1}\) values are indicated in Figure 4(a) by the connecting arrows. Compared to the 2D ICE plot visualization, the 3D IANN visualization provides a clearer picture of the function \(f\) in this case.

Figure 4: Connection between (a) existing methods (ICE and PD plots) and (b) our IANN for visualizing the effects of the amplitude variable \(x_{1}\) in the harmonic wave function. The IANN visualization can be viewed as smoothly piecing together the individual curves in the ICE plots to create a more interpretable 3D surface. The three dashed curves in (b) correspond to the three ICE plot curves indicated by connecting arrows in (a).
## 3 Hierarchical IANN for Visualizing the Joint Effects of All Inputs

In addition to being used as a stand-alone visualization of the effect of a selected input \(x_{j}\), the IANN structure described in Section 2 can be used in a hierarchical manner to similarly decompose and visualize \(h\), and so on. We present two versions of hierarchical decomposition, each of which is intended to visualize the joint effects of all \(d\) inputs. In Section 3.1, we introduce the original variable hierarchical (OVH) structure, and in Section 3.2 we introduce the disjoint active subspace hierarchical (DASH) structure, which is in some sense a generalization of the OVH structure that can simplify visualization when the various functions involved are functions of certain linear combinations of the inputs.

### Original variable hierarchical (OVH) structure

In the structure (1), \(h\) is a function of all the inputs in \(\boldsymbol{x}_{\backslash j}\). To visualize how \(h\) depends on \(\boldsymbol{x}_{\backslash j}\), we can decompose the functions hierarchically (adding subscripts to \(g\) and \(h\) to indicate the hierarchical level) via

\[f(\boldsymbol{x}) \approx g_{1}(x_{j_{1}},h_{1}(\boldsymbol{x}_{\backslash j_{1}}))\] \[h_{i-1}(\boldsymbol{x}_{\backslash(j_{1},\cdots,j_{i-1})}) \approx g_{i}(x_{j_{i}},h_{i}(\boldsymbol{x}_{\backslash(j_{1}, \cdots,j_{i})})),\quad i=2,\cdots,d-1, \tag{5}\]

with \(h_{d-1}(x_{j_{d}})\stackrel{{\text{\tiny def}}}{{=}}x_{j_{d}}\). Here, \(\{j_{1},j_{2},\cdots,j_{d}\}\) represent a permutation ordering of the input indices \(\{1,2,\cdots,d\}\) that will be determined as a preprocessing step prior to fitting the hierarchical IANN. The algorithm for determining the ordering of the inputs is described in Section 4. Similar to the notation \(\boldsymbol{x}_{\backslash j}\), the notation \(\boldsymbol{x}_{\backslash(j_{1},\cdots,j_{i})}\) represents the input variables excluding \(\{x_{j_{1}},\cdots,x_{j_{i}}\}\). Each function \(g_{i}\) is the approximation of the function \(h_{i-1}\) with two inputs: \(x_{j_{i}}\) and \(h_{i}(\boldsymbol{x}_{\backslash(j_{1},\cdots,j_{i})})\). Each function \(h_{i}\) depends on the order \(\{j_{1},j_{2},\cdots,j_{i}\}\) of the inputs selected in the prior iterations, although we omit this in the notation for
simplicity. Since the function \(h\) in Theorem 2.1 is also continuous, we can repeatedly apply Theorem 2.1 to the function \(h_{i-1}\) at each level to justify the hierarchical decomposition (5).

We refer to the first equation in (5) as the _Level 1 representation_ and the corresponding 3D plot or 2D heatmap (\(f(\mathbf{x})\) as a function of \(x_{j_{1}}\) and \(h_{1}(\mathbf{x}_{\setminus j_{1}})\)) as the _Level 1 plot_. Similarly, we refer to the \(i^{th}\) equation in (5) as the _Level \(i\) representation_ and its corresponding 3D plot or 2D heatmap (\(h_{i-1}\) as a function of \(x_{j_{i}}\) and \(h_{i}\)) as the _Level \(i\) plot_, for \(i=2,\cdots,d-1\). This approach produces a total of \(d-1\) hierarchical visualization plots that can collectively be used to understand the joint effects of all \(d\) inputs.

To estimate all the \(h_{i}\) functions simultaneously, we use an IANN architecture with multiple bottle-neck layers as illustrated in Figure 5 for the case of \(d=5\) inputs. The dashed arrows represent the identity map. To simplify the notation, we omit the hat symbols on the \(h_{i}\)'s in the IANN architecture. In Figure 5, the IANN architecture consists of 5 input variables and 4 levels in total. The fourth level is just a fully connected neural network with two input variables and an output that represents the (estimated) function \(h_{3}(x_{j_{4}},x_{j_{5}})\). Likewise, the third level has bottleneck input layer with two inputs (\(x_{j_{3}}\) and \(h_{3}\)) and \(h_{2}\) as its output, the second level has bottleneck input layer with two inputs (\(x_{j_{2}}\) and \(h_{2}\)) and \(h_{1}\) as its output, and the first level has \(x_{j_{1}}\) and \(h_{1}\) as its two bottleneck inputs and the estimated/approximated \(f\) as its output.

Similar to the loss function defined in (4), we use \(L_{2}\) loss to minimize the mean squared error between the IANN output \(\hat{f}\) and the original function \(f\). By estimating the original function \(f(\mathbf{x})\) and all the latent functions \(h_{i}\) simultaneously, the error is less prone to accumulate through the different levels when \(d\) is larger.

To illustrate the OVH IANN approach, consider the following example:

\[f(\mathbf{x})=(5x_{1}+x_{2}+x_{3}+x_{4}+x_{5}-4.5)^{2},\qquad\mathbf{x}\in[0,1]^{5}. \tag{6}\]

In this case, the input ordering \(\{j_{1},j_{2},\cdots,j_{5}\}\) that satisfies the structure (5) is not unique. In fact, any permutation of the input variables satisfies (5). For example, if we choose
\(j_{1}=2\), the response can be represented as:

\[f(\mathbf{x})=g_{1}(x_{2},h_{1}(\mathbf{x}_{\backslash 2}))=(x_{2}+h_{1}(\mathbf{x}_{ \backslash 2}))^{2}, \tag{7}\]

where \(h_{1}(\mathbf{x}_{\backslash 2})=5x_{1}+x_{3}+x_{4}+x_{5}-4.5\). Following the same logic, \(h_{1}\) can be decomposed similarly regardless of which remaining input is chosen as \(x_{j_{2}}\), and so on. Note that the above arguments would hold if any of the five inputs were selected as \(x_{j_{1}}\). In the following, we use the ordering \((j_{1},j_{2},j_{3},j_{4},j_{5})=(1,5,4,2,3)\), which was determined by the algorithm described in Section 4. The visualization 3D plots are shown in Figure 6. The test \(r^{2}\) was 99.97% indicating that the OVH structure (5) provides a good approximation of \(f\).

Figure 5: IANN architecture for the OVH structure for the case of \(d=5\) input variables.
The following are salient points taken from the IANN visualization plots in Figure 6. From the Level 1 plot the function \(f(\mathbf{x})\) appears quadratic in \(x_{1}\) when we hold \(h_{1}\) fixed, which is consistent with the true function \(f(\mathbf{x})=(5x_{1}-h_{1})^{2}\) if we take \(h_{1}=-(x_{2}+x_{3}+x_{4}+x_{5}-4.5)\). The functions \(g_{1}\) and \(h_{1}\) in (7) are not unique, of course, since we can incorporate an additive and/or multiplicative constant into \(h_{1}\) and modify \(g_{1}\) accordingly without changing the function \(g_{1}(x_{1},h_{1}(\mathbf{x}_{\backslash 1}))\). This, however, does not change the interpretation of the joint effects of the input variables. Further regarding the interpretation of the Level 1 plot, the function \(f\) reaches its minimum at a value of \(x_{1}\) that increases linearly as \(h_{1}\) increases.

In order to understand the effects of the other inputs \(\mathbf{x}_{\backslash 1}\) on \(f\), we must discern two things from the plots in Figure 6: (i) The effect of \(h_{1}\) on \(f\), and (ii) the effect of the variables in \(\mathbf{x}_{\backslash 1}\) on \(h_{1}\). Regarding the former, from the Level 1 plot, for fixed \(x_{1}\) the function \(f\) appears roughly quadratic in \(h_{1}(\mathbf{x}_{\backslash 1})\). In order to understand the effects of

Figure 6: IANN visualization plots for the example \(f\) in Eq. (6). The plots are to be read in numerical order of the levels, from the top left to the bottom right.
on \(h_{1}\), one must view the subsequent level plots. From the Level 2, 3, and 4 plots in Figure 6, each function \(h_{i}\) is approximately linear in its two arguments \(x_{j_{i+1}}\) and \(h_{i+1}\), which means that \(h_{1}\) is linear in \(\boldsymbol{x}_{\backslash 1}\).

To further illustrate how the IANN level plots can help understand the behavior of \(f\), supposed one wanted to select the values of \(\boldsymbol{x}\) that maximize \(f\). From the Level 1 plot, this occurs at the two corner points \((x_{1},h_{1})=(0,1.5)\) or \((1,-1)\). If we focus on the former, from the Level 2 plot, \(h_{1}=1.5\) occurs for \((x_{5},h_{2})=(0,2)\). In turn, from the Level 3 plot, \(h_{2}=2\) occurs for \((x_{4},h_{3})=(0,-0.5)\). Finally, from the Level 4 plot, \(h_{3}=-0.5\) occurs for \((x_{2},x_{3})=(0,0)\). Therefore, the maximum of \(f\) at \((x_{1},h_{1})=(0,1.5)\) corresponds to \(x_{1}=x_{2}=\cdots=x_{5}=0\). Following the same procedure for the other maximum at \((x_{1},h_{1})=(1,-1)\), we find that this corresponds to \(x_{1}=x_{2}=\cdots=x_{5}=1\). Both of these cases correspond to the true maxima of Eq. (6).

Moreover, since \(h_{1}(\boldsymbol{x}_{\backslash 1})\) is just a linear function of the remaining variables \(x_{2},x_{3},x_{4},x_{5}\) from the plots for Levels 2, 3, and 4, and since \(f(\boldsymbol{x})\) is quadratic in \(h_{1}(\boldsymbol{x}_{\backslash 1})\) and in \(x_{1}\), we can conclude that \(f\) is quadratic in each variable.

### Disjoint active subspace hierarchical (DASH) structure

In this section, we propose an alternative hierarchical structure that can be viewed as a generalization of the OVH structure. The DASH structure assumes the hierarchical functions \(h_{i}\)'s are functions of certain disjoint linear combinations of the input variables:

\[f(\boldsymbol{x}) \approx g_{1}(v_{1},h_{1}(\boldsymbol{v}_{\backslash 1})),\] \[h_{i-1}(\boldsymbol{v}_{\backslash(1,\cdots,i-1)}) \approx g_{i}(v_{i},h_{i}(\boldsymbol{v}_{\backslash(1,\cdots,i)})), \qquad i=2,\cdots,p-1, \tag{8}\] \[\text{where}\quad v_{i} =\boldsymbol{\beta}_{i}^{T}\boldsymbol{x}_{J_{i}},\quad i=1,2, \cdots,p\leq d.\]

Here \(\{J_{1},\cdots,J_{p}\}\) are disjoint index sets, i.e., \(J_{i}\bigcap J_{j}=\emptyset,\ \forall 1\leq i\neq j\leq p\), and \(\bigcup\limits_{i=1}^{p}J_{i}=\{1,2,\cdots,d\}\), and \(\boldsymbol{x}_{J_{i}}\) denotes the input variables with indices in \(J_{i}\). Note that the DASH
structure (8) reduces to the OVH structure when \(p=d\).

Theorem 2.1 also applies to structure (8) if we substitute the \(x\)'s with the \(v\)'s since the disjoint linear combinations of inputs (the \(v\)'s) take values in closed intervals. Similar to the proof for the OVH structure, one can repeatedly apply the theorem to \(f\) and the \(h_{i}\)'s to show that any continuous function \(f\) on \([0,1]^{d}\) can be arbitrarily closely approximated by the structure (8) for some continuous functions \(g_{i}\) and \(h_{i}\). The algorithms we describe in Appendix A automatically find the number of disjoint linear combinations \(p\), the \(v\)'s in (8), and the order of the \(p\) disjoint linear combinations \((v_{1},v_{2},\cdots,v_{p})\) that gives a good approximation of \(f\) by (8).

Similar to what was done for the OVH structure, we can make a series of 3D plots to visualize \(f\) in (8) as a function of the \(v\)'s, noting that the \(v\)'s are easy to understand as functions of the inputs, since they are disjoint linear combinations of the inputs. More specifically, in the Level 1 plot we visualize \(f\) as a function of \(v_{j_{1}}\) and \(h_{1}\), and then further visualize \(h_{i-1}\) as a function of \(v_{j_{i}}\) and \(h_{i}\) at the \(i^{th}\) level for \(i=2,3,\cdots,p-1\). For many real examples that we have considered, Eq. (8) with \(p<<d\) provides a close approximation of \(f\). The advantage of this is that we reduce the number of hierarchical 3D visualization plots from \(d-1\) to \(p-1\), which substantially simplifies the interpretation when \(p<<d\).

To represent (8) with an IANN architecture, we add one linear layer in front of the first layer in the OVH structure to learn the \(p\) disjoint linear combinations \((v_{1},v_{2},\cdots,v_{p})\) and the underlying disjoint input groups, as illustrated in Figure 7 for the case of \(p=5\). Each linear combination \(v_{i}=\boldsymbol{\beta}_{i}^{T}\boldsymbol{x}_{J_{i}}\) is represented by a single node in this layer with no bias and a linear activation function. The coefficients \(\boldsymbol{\beta}_{i}\) in (8) are the estimated weights in this linear layer. The remainder of the architecture is the same as in the IANN for the OVH structure shown in Figure 5, except that the inputs to the subsequent layers are the disjoint linear combinations. For visualization, we show the coefficients \(\boldsymbol{\beta}_{i}\) below the 3D visualization plot of \(h_{i-1}\) as a function of \(v_{j_{i}}\) and \(h_{i}\).
To illustrate, consider the function:

\[f(\mathbf{x})=\left[7\exp\left(-4(1.5x_{1}+x_{2}-2x_{3})^{2}\right)+2x_{4}-1.5x_{5}+0. 7x_{6}-1.5\right](x_{7}-1.5x_{8}+0.7x_{9}-0.3)^{2}, \tag{9}\]

for \(\mathbf{x}\in[0,1]^{9}\), which can be represented as (8) with \(p=3\) and

\[\begin{array}{l}v_{1}=x_{7}-1.5x_{8}+0.7x_{9},\\ v_{2}=2x_{4}-1.5x_{5}+0.7x_{6},\\ v_{3}=1.5x_{1}+x_{2}-2x_{3},\quad\text{and}\\ f(\mathbf{x})=(v_{1}-0.3)^{2}*h_{1}(v_{2},v_{3}),\end{array} \tag{10}\]

where \(h_{1}(v_{2},v_{3})=\left[7\exp\left(-4v_{3}^{2}\right)+v_{2}-1.5\right]\). We fit the DASH structure IANN to data from this function (the three disjoint linear combinations were automatically determined using the algorithm described in Appendix A and were not prespecified), the results of which are shown in Figure 8. The test \(r^{2}\) was 99.98%, indicating that the DASH structure (8) provides a good approximation of \(f\).

Rather than using eight 3D plots to visualize and interpret the func

Figure 7: Illustration of the IANN architecture for the DASH structure with \(p=5\).
the OVH structure, we only need two 3D plots to visualize the effects of all the input variables with the DASH structure in Figure 8. We show the coefficients \(\beta_{i,j}\) of the inputs \(x_{j}\) in each linear combination \(v_{i}\) below each IANN visualization plot. From Figure 8, we see that the estimated coefficients in each \(\boldsymbol{\beta}_{i}\) have ratios that are similar to the true ratios defined in Eq. (10). For example, the true \(v_{1}=x_{7}-1.5x_{8}+0.7x_{9}\), and the estimated \(v_{1}\) from fitting the IANN is \(v_{1}=0.667x_{7}-1.001x_{8}+0.467x_{9}\), which is virtually identical except for a constant multiplicative factor.

Regarding interpreting the plots, first, from the top plot in Figure 8 the function \(f(\boldsymbol{x})\) is non-monotonic (roughly quadratic) in \(v_{1}\) when \(h_{1}\) is fixed. Since \(v_{1}\) is a linear combination of three input variables, \(x_{7},x_{8},x_{9}\), we can conclude that each of these input variables have non-monotonic effects on \(f(\boldsymbol{x})\), at least for some fixed values of the other two variables in this linear combination. From the same plot, when \(v_{1}\) is fixed, function \(f(\boldsymbol{x})\) monotonically increases in a nearly linear manner as \(h_{1}\) increases. To visualize \(h_{1}\), the Level 2 plot suggests that \(h_{1}\) is monotonic and nearly linear in \(v_{2}\) for each fixed \(v_{3}\) but non-monotonic in \(v_{3}\) for

Figure 8: DASH structure IANN visualization of \(f\) in Eq. (9). The solid curves in the left figure are \(f\) vs \(v_{1}\) for fixed \(h_{1}=-1.5,-1.0\), and \(-0.5\).
each fixed \(v_{2}\). One conclusion from this is that \(f\) is nearly linear in \(x_{4}\), \(x_{5}\), and \(x_{6}\) (the inputs involved in \(v_{2}\)) with fixed \(v_{1}\) but nonlinear in the other input variables. Notice that the slight concave curvature in \(f\) with respect to \(h_{1}\) tends to cancel the slight convex curvature in \(h_{1}\) with respect to \(v_{2}\), so that \(f\) is nearly linear in \(v_{2}\) for \(v_{1}\) fixed.

Moreover, interactions between input variables can also be discerned from the plots. Understanding interactions (as opposed to additivity) between inputs is often important when interpreting models. We can tell whether an input variable in \(v_{1}\) interacts with the remaining variables from the Level 1 IANN plot. Since the remaining variables reside either in \(v_{1}\) or \(h_{1}\), we consider two types of interactions: (i) interaction between variables in \(v_{1}\) and in \(h_{1}\), and (ii) interaction for variables within \(v_{1}\). If neither interaction exists for that input variable, we can conclude that it has no interaction with all the remaining variables.

Regarding the interaction between \(v_{1}\) and \(h_{1}\), if the effect of \(v_{1}\) on \(f\) (e.g., the solid curves in the top left plot in Figure 8) only changes by an additive constant as \(h_{1}\) changes, then by definition, there is no interaction between \(v_{1}\) and \(h_{1}\) and thus no interaction between the inputs \(\mathbf{x}_{J_{1}}\) and the inputs \(\mathbf{x}_{\setminus J_{1}}\). Regarding the interactions within \(v_{1}\), if \(v_{1}\) has only linear effect on \(f\) for each \(h_{1}\) values, the input variables within \(v_{1}\) do not interact with each other. Conversely, if \(v_{1}\) has nonlinear effect on \(f\) for at least some \(h_{1}\) values, then there are interactions within \(v_{1}\) based on the following argument. In this case, for the \(h_{1}\) values for which the effect of \(v_{1}\) is nonlinear, \(\frac{\partial f(v_{1},h_{1})}{\partial v_{1}}=\tilde{g}(v_{1},h_{1})\) for some function \(\tilde{g}\) that varies as \(v_{1}\) varies, in which case the input variables \(\mathbf{x}_{J_{1}}\) having nonzero \(\beta\) coefficients all interact with each other.

From the IANN visualization in Figure 8, we see from the three solid curves in the Level 1 plot that there is a strong interaction between \(v_{1}\) and \(h_{1}\). Since both \(v_{2}\) and \(v_{3}\) have an effect on \(h_{1}\) from the second level plot, we conclude that the inputs in \(\mathbf{x}_{J_{1}}\) interact with those in \(\mathbf{x}_{J_{2}}\) and \(\mathbf{x}_{J_{3}}\). Moreover, since the solid curves in the Level 1 plot are nonlinear, the inputs in \(\mathbf{x}_{J_{1}}\) interact with each other.
As another example of how the IANN plots can be used, consider the so-called robust parameter design (RPD) problem (Taguchi, 1986; Robinson et al., 2004), in which some of the input variables are controllable system/product design variables and others are noise variables that vary uncontrollably during system operation or product usage. The RPD goal is to find design values for the controllable inputs such that the function \(f\) is least sensitive to variation in the uncontrollable inputs. To illustrate, suppose the inputs \(\mathbf{x}_{J_{1}}=\{x_{7},x_{8},x_{9}\}\) comprising \(v_{1}\) are noise variables, and the other inputs are controllable design variables. From the Level 1 plot in Figure 8, we see that \(f\) depends least strongly on \(v_{1}\) when \(h_{1}\approx-1.0\). Consequently, if the inputs \(\mathbf{x}_{\setminus J_{1}}\) on which \(h_{1}\) depends are selected so that \(h_{1}\approx-1.0\) (e.g., by making use of the Level 2 plot in Figure 8), then \(f\) will be least sensitive to variation in the noise variables \(\mathbf{x}_{J_{1}}\).

## 4 IANN Algorithm Details

This section describes how to determine the order of the input variables in the hierarchical IANN decomposition (i.e., which input variable appears at each level) and other details of the algorithm. Here, we mainly present algorithms for this purpose for the OVH IANN structures. The algorithms for the DASH IANN structures, the details of which can be found in Appendix A, are similar.

To best approximate \(f\) with the OVH structure and make the visualization plots more reliable for interpretation, we determine the order of input variables sequentially by selecting the variable that results in the best approximation accuracy at each level in the OVH structure (5). To see the intuition behind our approach for this, notice that if \(f(\mathbf{x})=g_{1}(x_{j_{1}},h_{1}(\mathbf{x}_{\setminus j_{1}}))\) in the first level of (5), then by the chain rule,

\[\frac{\partial f}{\partial\mathbf{x}_{\setminus j_{1}}}(\mathbf{x})=\frac{\partial g_ {1}}{\partial\mathbf{x}_{\setminus j_{1}}}(\mathbf{x})=\frac{\partial g_{1}}{\partial h _{1}}(\mathbf{x})\nabla h_{1}(\mathbf{x}_{\setminus j_{1}}), \tag{11}\]

where \(\nabla\) denotes the gradient of a function with respect to its input arguments. Since
\(\frac{\partial g_{1}}{\partial h_{1}}(\mathbf{x})\) is a scalar, if we keep \(\mathbf{x}_{\setminus j_{1}}\) fixed and consider \(\frac{\partial f}{\partial\mathbf{x}_{\setminus j_{1}}}(\mathbf{x})\in\mathbb{R}^{d-1}\) for many different \(x_{j_{1}}\) values, they are all approximately colinear and differ only in their magnitude and sign. Therefore, as candidates for \(x_{j_{1}}\), we consider the inputs whose gradients are the most colinear using principal component analysis (PCA) to measure the extent of colinearity.

We note that although the functions \(g\) and \(h\) from the IANN structure appear in Eq. (11) (and in Eq. (16) below), they are not used in the algorithm for determining the ordering of the inputs in the IANN, which is a pre-processing step that uses the gradient of \(f\) directly and requires no IANN fitting. Eqs. (11) and (16) are only used to justify the rationale behind our input ordering approach. We also note that finding the ordering of inputs \(\mathbf{x}\) requires the function \(f\) to be differentiable since we need to compute the gradient of \(f(\mathbf{x})\) with respect to \(\mathbf{x}\).

The following describes the details of this gradient projection algorithm for Level 1. **Gradient projection algorithm:** (repeat steps 1-3 for each input variable \(x_{j},\ j=1,2,\cdots,d\))

1. For input variable \(x_{j}\), draw \(N_{\setminus j}\) latin hypercube design (LHD) samples in the \(\mathbf{x}_{\setminus j}\) space, denoted by \(\{\mathbf{x}_{\setminus j}^{l};\ l=1,\cdots,N_{\setminus j}\}\). Then draw \(N_{j}\) evenly spaced points spanning the range of \(x_{j}\) in the \(x_{j}\) space, denoted by \(\{x_{j}^{m};\ m=1,\cdots N_{j}\}\).
2. For \(l=1,2,\cdots,N_{\setminus j}\), calculate the gradient vectors \(\mathbf{a}_{j}^{m,l}=\nabla_{\mathbf{x}_{\setminus j}}f(x_{j}^{m},\mathbf{x}_{\setminus j }^{l})\in\mathbb{R}^{d-1},\ m=1,2,\cdots,N_{j}\), and stack them into the \(N_{j}\times(d-1)\) matrix: \[A_{j}^{l}=\begin{bmatrix}\left(\mathbf{a}_{j}^{1,l}\right)^{T}\\ \vdots\\ \left(\mathbf{a}_{j}^{N_{j},l}\right)^{T}\end{bmatrix},\quad\text{for }l=1,2, \cdots,N_{\setminus j}.\] (12)
3. Use PCA to find the eigenvector \((\mathbf{z}_{j}^{l})\) corresponding to the largest eigenvalue of \(\left(A_{j}^{l}\right)^{T}A_{j}^{l}\) for each \(l=1,2,\cdots,N_{\setminus j}\), and normalize the eigenvector such that \(\left|\left|\mathbf{z}_{j}^{l}\right|\right|=1\). Compute the error between the gradient vectors \(\mathbf{a}_{j}^{m,l}\) and their projections onto
\(\mathbf{z}_{j}^{l}\), defined as \[E_{j}^{l}=\sum_{m=1}^{N_{j}}\left|\left|\mathbf{a}_{j}^{m,l}-\left[\mathbf{z}_{j}^{l} \left(\mathbf{z}_{j}^{l}\right)^{T}\right]\mathbf{a}_{j}^{m,l}\right|\right|^{2},\quad \text{for }l=1,2,\cdots,N_{\setminus j}.\] (13) Then, compute the gradient projection error \((E_{j})\) for \(x_{j}\) via \[E_{j}=\frac{1}{N_{\setminus j}}\sum_{l=1}^{N_{\setminus j}}\frac{E_{j}^{l}}{ \Omega_{j}^{l}},\] (14) where \[\Omega_{j}^{l}=\frac{1}{N_{j}}\sum_{m=1}^{N_{j}}\left|\left|\mathbf{a}_{j}^{m,l} \right|\right|^{2}\qquad\text{for }l=1,2,\cdots,N_{\setminus j}.\] (15) Note that we use the normalization factor \(\Omega_{j}^{l}\) in (14) to express the projection error \(E_{j}^{l}\) relative to the average squared length of the gradient vectors for that \(l\).
4. Finally, we choose the input variable \(x_{j_{1}}\) as the one having the smallest projection error, i.e., \(j_{1}=\underset{j}{\arg\min}\{E_{j}:j=1,2,\cdots,d\}\).

We use the maximin criterion (Johnson et al., 1990; McKay et al., 2000) to draw the LHD samples throughout the paper, which aims to maximize the minimum distance between any two samples in the LHD to enhance the space-filling property. Also note that the entire IANN modeling procedure uses LHDs for two different purposes: First, in the gradient projection algorithm above for determining the input ordering, and then again to generate the training data for the IANN model fitting (see Appendix B). For the former, we have found that a standard LHD works fine. For the latter, we have found that the modified LHD in Appendix B consistently works better than a standard LHD.

The algorithm for selecting the input for each subsequent level is similar: Suppose we have determined the order of input variables in the previous levels, \(\{x_{j_{1}},\cdots,x_{j_{i-1}}\}\) for some \(i>1\). Similar to the procedure for finding \(x_{j_{1}}\), for each candidate input index \(j_{i}\in\{1,2,\cdots,d\}\setminus\{j_{1},\cdots,j_{i-1}\}\), we consider the gradient of \(f(\mathbf{x})\) with respect to the remaining inputs \(\mathbf{x}_{\setminus(j_{1},\cdots,j_{i})}\), which by the chain rule is (if structure (5) holds with an equality),
\[\frac{\partial f}{\partial\mathbf{x}_{\setminus(j_{1},\cdots,j_{i})}}(\mathbf{x})=\frac{ \partial g_{1}}{\partial h_{1}}(\mathbf{x})\frac{\partial h_{1}}{\partial h_{2}}( \mathbf{x}_{\setminus j_{1}})\cdots\frac{\partial h_{i-1}}{\partial h_{i}}(\mathbf{x} _{\setminus(j_{1},\cdots,j_{i-1})})\nabla h_{i}(\mathbf{x}_{\setminus(j_{1}, \cdots,j_{i})}). \tag{16}\]

Then, if we fix \(\mathbf{x}_{\setminus(j_{1},\cdots,j_{i})}\) and vary \((x_{j_{1}},\cdots,x_{j_{i}})\) by taking some LHD samples in the \(\mathbf{x}_{(j_{1},\cdots,j_{i})}\) space, the vector \(\nabla h_{i}(\mathbf{x}_{\setminus(j_{1},\cdots,j_{i})})\) in (16) is a constant vector, so that the gradient vectors \(\frac{\partial f}{\partial\mathbf{x}_{\setminus(j_{1},\cdots,j_{i})}}(\mathbf{x})\) are all colinear and differ only in their magnitude and sign. Consequently, we apply the same gradient projection algorithm described above to the subsequent levels by substituting \(\mathbf{x}_{\setminus j_{1}}\) with \(\mathbf{x}_{\setminus(j_{1},\cdots,j_{i})}\) and selecting \(x_{j_{i}}\) to be the remaining input whose gradient vectors in Eq. (16) are the most colinear according to the gradient projection error measure analogous to (14). To avoid generating LHD samples repeatedly in each level, we only take LHD samples once in the \(d\)-dimensional input space and then project them into the lower dimensional spaces where needed.

## 5 Numerical Experiments

This section provides two examples of visualizing functions for which we have a closed-form "ground truth" expression to compare to our IANN visualizations. The Supplementary Materials section provides an example in which the function is an actual surrogate model of a complex numerical simulation of the potential energy of a strain-loaded material sample.

### IANN for the borehole function

To illustrate how IANN plots help with the visualization and interpretation of function, we use the borehole function

\[f(\mathbf{x})=2\pi x_{1}\left(x_{4}-x_{6}\right)\left[\log\left(\frac{x_{2}}{x_{3 }}\right)\left(1+2\frac{x_{7}x_{1}}{\log(\frac{x_{2}}{x_{3}})x_{3}^{2}x_{8}} +\frac{x_{1}}{x_{5}}\right)\right]^{-1}, \tag{17}\]

which models the water flow between two aquifers. Although the functional form in (17) is known, it is commonly treated as a black-box function for evaluating surrogate modeling methods (An and Owen, 2001; Harper and Gupta, 1983; Morris et al., 1993; Surjanovic and
Bingham, 2013). The response \(f\) is the water flow rate between the aquifers. Here, \(d=8\), \(x_{1}\in[0.05,0.15]\) and \(x_{2}\in[100,50\,000]\) represent the radius of a borehole and its influence respectively, \(x_{3}\in[63\,070,115\,600]\) and \(x_{5}\in[63.1,116]\) represent transmissivity of the upper and lower aquifer, \(x_{4}\in[990,1110]\) and \(x_{6}\in[700,820]\) represent the potentiometric head of the upper and lower aquifer, \(x_{7}\in[1120,1680]\) represents the length of the borehole, and \(x_{8}\in[9855,12\,045]\) represents the hydraulic conductivity of the borehole.

Since the borehole function has eight inputs, each with different scales, we first use the min-max normalization to rescale the range of inputs to \([0,1]\) when fitting the IANN model. The OVH structure will hierarchically generate seven 3D visualization plots. Instead, we use the DASH structure to interpret the effects of all the input variables through four linear combinations and their hierarchical ordering, which were automatically determined using the algorithm described in Appendix A, which is similar to the one in Section 4.

Regarding computational expense, the algorithm took 96 seconds to determine the input ordering for the borehole example. The main computational expense is in fitting the single IANN after determining the input ordering, which took 245 seconds for the borehole example. The computational expense of fitting the single IANN is comparable to that of fitting standard neural networks with comparable numbers of parameters. One can apply cross-validation to choose the optimal hyperparameters in our IANN model, the computational expense of which is proportional to the time for fitting a single IANN.

In the event that users want to understand and interpret the effect of one input in particular, they can manually choose the ordering of the linear combination groups so that the first group contains the input variable of interest. Table 1 shows the ordered groups produced by the algorithms for the DASH IANN structure with the constraint that the first group is the one containing the input of particular interest. The eight rows show the resulting orderings with each of the eight inputs specified as being of particular interest. This way, one can more directly visualize its effect on \(f\) from the Level 1 plot. The order
of the remaining linear combinations is shown in the second column in Table 1. For each of the orderings in Table 1, we also show the resulting test \(r^{2}\). Figure 9 shows the IANN visualization plot with the highest test \(r^{2}\).

For comparison, we use global sensitivity analysis (Herman and Usher, 2017) to compute the total sensitivity indices for \(\{x_{1},x_{2},\cdots,x_{8}\}\), which are \(\{0.0,0.0,0.174,0.26,0.0,0.261,\\ 0.258,0.063\}\). This suggests that \(x_{1},x_{2}\) and \(x_{5}\) have little effect on the borehole function, which is also reflected in the IANN plot in Figure 9.

We used 34,000 training samples for fitting our IANN and \(10^{6}\) randomly generated test input samples to compute the test \(r^{2}\). The test \(r^{2}=99.99\%\) for the Borehole function. The coefficients of the (normalized) inputs in each linear combination are listed below the plots in Figure 9. Since the inputs have been normalized to \([0,1]\), comparing the magnitudes of

\begin{table}
\begin{tabular}{|c|c|c|} \hline index j & input groups and ordering & test \(r^{2}\) \\ \hline
1 & [[1, 2, 3, 5], [4, 6], [7], [8]] & 99.97\% \\ \hline
2 & [[1, 2, 3, 5], [4, 6], [7], [8]] & 99.97\% \\ \hline
3 & [[1, 2, 3, 5], [4, 6], [7], [8]] & 99.97\% \\ \hline
4 & [[4, 6],[[1, 2, 3, 5], [7], [8]]] & 99.99\% \\ \hline
5 & [[1, 2, 3, 5], [4, 6], [7], [8]]] & 99.97\% \\ \hline
6 & [[4, 6],[[1, 2, 3, 5], [7], [8]]] & 99.99\% \\ \hline
7 & [[7],[4, 6],[[1, 2, 3, 5], [8]]] & 99.94\% \\ \hline
8 & [[8],[4, 6],[[1, 2, 3, 5], [7]]] & 99.99\% \\ \hline \end{tabular}
\end{table}
Table 1: For the borehole function with the DASH structure assumed, a listing of the \(p=4\) disjoint linear combination groups and their orderings, which were automatically determined by algorithms for the DASH IANN structure under the constraint that the first group contains the specified input variable of particular interest, which is shown in Column 1.
the coefficients suggests the relative importance of the inputs involved in each linear combination. For example, from the Level 1 plot in Figure 9, \(v_{1}\) consists of two input variables, \(x_{4}\) and \(x_{6}\), and the corresponding coefficients have the same magnitude with different signs. This agrees with the formula for the borehole function in Eq. (17), which can be written as \(f(\mathbf{x})=(x_{4}-x_{6})h_{1}(\mathbf{x}_{\setminus(4,6)})\) for \(h_{1}=2\pi x_{1}\left[\log\left(\frac{x_{2}}{x_{3}}\right)\left(1+2\frac{x_{ 7}x_{1}}{\log(\frac{x_{2}}{x_{3}})x_{3}^{2}x_{8}^{3}}+\frac{x_{1}}{x_{5}} \right)\right]^{-1}.\) From the Level 2 plot, the coefficients are almost zero except for the input variable \(x_{3}\). This also agrees with the global sensitivity analysis, which shows that \(x_{1},x_{2},x_{5}\) have little impact on the response \(f\).

As a reference for comparison, the test \(r^{2}\) for a linear model fitted to the borehole function is \(94.68\%\), which suggests that the borehole function is almost linear over the specified input domain. The examples in the subsequent sections demonstrate the interpretability of the IANN visualization for functions having a higher level of nonlinearity.

Figure 9: DASH structure IANN visualization for the borehole function \(f\) in Eq. (17). The axes of the plots and the \(\mathbf{\beta}\) coefficients listed below the plots correspond to the normalized inputs.
### IANN for the harmonic wave function

Reconsider the harmonic wave function \(f(\mathbf{x})\) in Eq. (2), for which the IANN visualizations for only Level 1 were shown in Figures 1 and 3. We first considered the DASH structure to approximate \(f\), but it found no relevant disjoint linear combinations other than the trivial one with \(p=d\) and each \(v_{i}\) a single input variable. Consequently, we use the OVH structure. More generally, we recommend that users first try the DASH structure, and if it does not produce any relevant linear combinations that result in an acceptably high test \(r^{2}\), then users should use the OVH structure.

From the Level 1 plot in Figure 10 (which is similar to Figure 1), we see that \(x_{1}\) has linear effect on \(f(\mathbf{x})\) when the other variables are fixed, and \(f(\mathbf{x})\) is sinusoidal for each fixed value of \(x_{1}\), which agrees with Eq. (2). From the Level 2 plot, \(h_{1}(\mathbf{x}_{\setminus 1})\) is linear in \(x_{4}\) and \(h_{2}\), and \(x_{4}\) and \(h_{2}\) have little interaction, which suggests the additive structure \(h_{1}(\mathbf{x}_{\setminus 1})=\alpha_{1}h_{2}(\mathbf{x}_{\setminus(1,4)})+\alpha_{2 }x_{4}\) for some constants \(\alpha_{1}\) and \(\alpha_{2}\). Therefore, \(x_{4}\) reflects the phase angle in the harmonic function. From the Level 3 plot, we can see that \(h_{2}\) is linear in \(x_{3}\) for each fixed \(x_{2}\) and nonlinear in \(x_{2}\) for some fixed \(x_{3}\), especially when \(x_{3}\) is higher. Moreover, it is clear from the Level 3 plot that \(x_{2}\) has zero effect on \(h_{2}\) when \(x_{3}=0\), versus a large and nonlinear effect on \(h_{2}\) when \(x_{3}=1\), which also agrees with the "true"

Figure 10: OVH structure IANN visualization of harmonic wave function in Eq. (2). The test \(r^{2}=99.95\%\).
expression for \(h_{2}\) (\(=2\pi\frac{x_{3}}{x_{2}}\) up to some constants) from Eq. (2). Likewise from the Level 3 plot, as we decrease \(x_{2}\), \(x_{3}\) has a far larger effect on \(h_{2}\), and thus a far larger effect on \(h_{1}\) (from the Level 2 plot) and on \(f\) (from the Level 1 plot). From this, for smaller \(x_{2}\), \(f\) undergoes more sinusoidal cycles as \(x_{3}\) varies over its full range, which agrees with the observation from Eq. (2) that \(x_{2}\) represents the wavelength of the sinusoidal function. In this sense, the IANN plot truly reflects the physical interpretation of the input variables.

One can find more examples in the Appendix C that use our IANN approach to visualize actual black box functions that represent a surrogate model fit to the output of a black-box computer simulation of a physical system.

## 6 Conclusions and Potential Extensions

This paper introduced the IANN approach to visualize and interpret black box functions. We have shown that, theoretically, any continuous function on \([0,1]^{d}\) can be arbitrarily closely approximated by our more interpretable structure, which can be conveniently represented with the proposed IANN architecture. To visualize the effects of all the input variables, we developed two hierarchical structures (OVH and DASH), each of which can be represented with a particular IANN architecture. We have also developed algorithms to automatically determine the ordering of input variables with the goal of providing the best approximation to the original function for each hierarchical structure. We have used a number of examples to demonstrate the interpretability advantages of our IANN method.

We envision several potential extensions of our IANN approach to either enhance the interpretation of the original function \(f\) and/or to expand the class of \(f\) to which our IANN provides a good approximation. One way to enhance interpretability with the IANN visualization plots is to develop customized graphical user interfaces (GUIs), with which users can change the value of each input variable \(\{x_{1},x_{2},\cdots,x_{d}\}\) via slide-bar controls and interactively visualize the corresponding points \(\{f,h_{1},\cdots,h_{d-2}\}\) on the IANN visualization
surfaces plotted in each level. Moreover, to facilitate the creation of 3D trellis plots of \(f\) vs \(\{x_{j},x_{l}\}\) for some collection of fixed values of \(\mathbf{x}_{\setminus(j,l)}\) (which normally requires careful selection of many fixed \(\mathbf{x}_{\setminus(j,l)}\) values when \(d\) is larger), one could modify the IANN architecture to fit a model of the form

\[f(\mathbf{x})=g(x_{j},x_{l},h(\mathbf{x}_{\setminus(j,l)})), \tag{18}\]

in which case users need only to select a few fixed values of the _scalar_\(h(\mathbf{x}_{\setminus(j,l)})\), instead of selecting many fixed values of the higher-dimensional \(\mathbf{x}_{\setminus(j,l)}\).

One might also consider removing the disjoint restriction on the linear combinations in the DASH structure, which would allow each input to appear in multiple combinations. The main reason we did not pursue this is that we are prioritizing interpretability of the effects of individual inputs over generality of the IANN approach. Having overlapping linear combinations would mean that the same input variable is present in multiple linear combinations, which makes the interpretation of that individual input variable less clear. Considering that the disjoint set restriction in the DASH structure becomes less restrictive and more general as the number of linear combinations increases (in the limiting case, each linear combination is a single input, in which case the DASH structure reduces to the OVH structure), we thought the disjoint restriction constitutes a reasonable tradeoff between generality and interpretability. As future work, we plan to explore this issue further, but as yet we are unsure how to preserve interpretability with overlapping linear combinations.

A related extension is to use the modified architecture (18) and simply plot \(f\) vs \(\{x_{j},x_{l}\}\) with \(h(\mathbf{x}_{\setminus(j,l)})\) controlled by a slide bar, which would be interpreted similarly to 3D trellis plots but completely avoids having to select any fixed values related to the omitted variables \(\mathbf{x}_{\setminus(j,l)}\). Analogous to the OVH structure in Section 3.1, one could further decompose \(h\) in Eq. (18) as a function of two additional inputs and a second latent function that is likewise controlled by a slide-bar in a Level 2 plot of \(h\) versus the two additional inputs, and so on in subsequent levels. This extension would reduce the number of required levels by roughly
half in the IANN visualizations and would also expand the applicability of the approach from functions that can be represented by (1) to functions that can be represented by (18).

Another potential extension is to use a dichotomous tree IANN structure, for which the Level 1 plot is of \(f(\mathbf{x})=g_{1}(h_{1,1}(\mathbf{x}_{J_{1,1}}),h_{1,2}(\mathbf{x}_{J_{1,2}}))\) vs \(h_{1,1}(\mathbf{x}_{J_{1,1}})\) and \(h_{1,2}(\mathbf{x}_{J_{1,2}})\), where \(J_{1,1}\) and \(J_{1,2}\) are a disjoint partition of \(\{1,2,\cdots,d\}\). Each latent \(h\) function can then be decomposed similarly in the subsequent levels to produce a tree-like nested set of 3D plots. To find the indices group partition at each level, we anticipate that something similar to the IANN algorithm in Section 4 can be developed. The main benefit of this extension is that it expands the class of \(f\) to which our IANN provides a good approximation (the IANN structure in (5) is a special case of the tree IANN structure with each group of inputs partitioned into a single input and the remaining inputs) and potentially reduces the depth of the tree in the IANN visualization.

Another potential extension is from the current setting of visualizing black box simulation functions to visualizing general supervised learning models fit to observational data, in order to interpret the effects of the input variables. The main challenge is that, unlike black box simulation functions for which the entire (rectangular) input space is meaningful, the input training data for supervised learning models are often highly correlated. Visualizing black box supervised learning models using our current IANN approach would require extrapolation to regions of the input space where data are scarce, which would render the interpretations unreliable.

## Supplementary Materials

**Supplemental Materials**: Several topics will be covered here, including the algorithms for the DASH IANN structure, customized LHD sampling techniques, additional numerical examples, and the proof of Theorem 2.1. (IANN-supplementary.pdf)
**Python-package for IANN:**: Python-package "IANN" containing code to perform the
IANN method described in the article and the related datasets. (iann-codes.zip, zipped codes for IANN)

## Disclosure Statement

The authors report there are no competing interests to declare.

## Acknowledgements

This work was funded in part by the Air Force Office of Scientific Research Grant # FA9550-18-1-0381, which we gratefully acknowledge.
