ANOVAapprox – interpretable approximation of high-dimensional functions

ANOVAapprox.jl (Julia) Julia documentation pyANOVAapprox (Python) Python documentation

ANOVAapprox approximates high-dimensional functions from scattered data \(\mathcal X = \{\boldsymbol x_1,\dots,\boldsymbol x_M\}\) and values \(y_j \approx f(\boldsymbol x_j)\). It assumes that \(f\) is dominated by interactions of only a few variables, i.e., that it has a low superposition dimension or a sparse ANOVA decomposition. The method is used for regression as well as for binary classification. Besides the approximation, it yields an interpretable model: global sensitivity indices, an attribute ranking and Shapley values show which variables and variable interactions are important. The fast matrix-vector multiplications are based on the NFFT and the NFCT.

ANOVA decomposition

Let \(\{\phi_{\boldsymbol k}\}_{\boldsymbol k}\) be a complete orthonormal system of tensor product form in \(L_2(\mathbb D^d)\), e.g., \(\mathbb D = \mathbb T\) for periodic and \(\mathbb D = [0,1]\) for nonperiodic functions, with \(\phi_{\boldsymbol 0} = 1\). With the coordinate indices \(\mathcal D = \{1,\dots,d\}\), a function \(f = \sum_{\boldsymbol k} c_{\boldsymbol k}(f)\, \phi_{\boldsymbol k}\) is uniquely decomposed into \(2^d\) ANOVA terms \[ f = \sum_{\boldsymbol u \subseteq \mathcal D} f_{\boldsymbol u}, \qquad f_{\boldsymbol u}(\boldsymbol x_{\boldsymbol u}) = \sum_{\operatorname{supp} \boldsymbol k = \boldsymbol u} c_{\boldsymbol k}(f)\, \phi_{\boldsymbol k}(\boldsymbol x), \] where \(\operatorname{supp} \boldsymbol k = \{s \in \mathcal D : k_s \neq 0\}\). The term \(f_{\boldsymbol u}\) depends only on the variables \(\boldsymbol x_{\boldsymbol u}\) and describes their joint interaction. Its importance is measured by the global sensitivity index \[ \varrho(\boldsymbol u, f) = \frac{\sigma^2(f_{\boldsymbol u})}{\sigma^2(f)}, \qquad \sigma^2(f) = \sum_{\boldsymbol k \neq \boldsymbol 0} |c_{\boldsymbol k}(f)|^2 . \]

To circumvent the curse of dimensionality, the decomposition is truncated to a downward closed set \(U\) of ANOVA terms, typically to all interactions of at most \(d_s\) variables, \[ \mathrm T_{U} f = \sum_{\boldsymbol u \in U} f_{\boldsymbol u}, \qquad U_{d_s} = \{\boldsymbol u \subseteq \mathcal D : |\boldsymbol u| \le d_s\}. \] The number of terms in \(U_{d_s}\) grows only polynomially in \(d\). Functions of dominating mixed smoothness with product and order-dependent weights have a low superposition dimension, so that this truncation causes only a small error.

Approximation from scattered data

Each term \(f_{\boldsymbol u}\) is approximated by a partial sum with frequencies \(\boldsymbol k_{\boldsymbol u}\) from a \(|\boldsymbol u|\)-dimensional index set \(I_{\boldsymbol u}\) with bandwidth \(N_{\boldsymbol u}\). The coefficients \(\hat{\boldsymbol f}\) of all terms, with index set \(I(U)\), are computed by a regularised least squares problem \[ \min_{\hat{\boldsymbol f}} \; \big\| \boldsymbol y - \boldsymbol F(\mathcal X, I(U))\, \hat{\boldsymbol f} \big\|_2^2 + \lambda \big\| \boldsymbol W \hat{\boldsymbol f} \big\|_2^2 , \qquad \boldsymbol F(\mathcal X, I(U)) = \big(\phi_{\boldsymbol k}(\boldsymbol x)\big)_{\boldsymbol x \in \mathcal X,\, \boldsymbol k \in I(U)}, \] with the LSQR method. The diagonal matrix \(\boldsymbol W\) of Sobolev-type weights incorporates a-priori smoothness information. Alternatively, a group lasso penalty \(\lambda \sum_{\boldsymbol u \in U} \|\hat{\boldsymbol f}_{\boldsymbol u}\|_{\boldsymbol W}\), solved by FISTA, promotes sparsity in the ANOVA terms and thus detects the important terms during the minimisation.

Both solvers need fast multiplications with \(\boldsymbol F(\mathcal X, I(U))\) and its adjoint. Since a term \(f_{\boldsymbol u}\) only depends on \(\boldsymbol x_{\boldsymbol u}\), the matrix splits into blocks that act in low dimensions, \[ \boldsymbol F(\mathcal X, I(U))\, \hat{\boldsymbol f} = \sum_{\boldsymbol u \in U} \boldsymbol F\big(\mathcal X_{\boldsymbol u}, I_{\boldsymbol u}\big)\, \hat{\boldsymbol f}_{\boldsymbol u}, \qquad \mathcal X_{\boldsymbol u} = \{\boldsymbol x_{\boldsymbol u} : \boldsymbol x \in \mathcal X\}, \] and every block is a \(|\boldsymbol u|\)-dimensional NFFT or NFCT. This grouped transform requires \(\mathcal O\big(\sum_{\boldsymbol u \in U} (N_{\boldsymbol u}^{|\boldsymbol u|} \log N_{\boldsymbol u} + |\log \varepsilon|^{|\boldsymbol u|} M)\big)\) arithmetic operations for the accuracy \(\varepsilon\), independently of the full dimension \(d\), and the blocks can be computed in parallel.

Bases

The following orthonormal systems are available (the name is the parameter basis of the packages). For cos, cheb and std, the univariate system consists of the constant function \(1\) and the functions listed for \(k \in \mathbb N\).

The basis std is the cosine basis composed with a transformation of \(\mathbb R\) onto \([0,1]\) and is suited for Z-score normalised data. The mixed basis allows functions that are periodic only in some variables; its fast transform (NFMT) combines the NFFT and the NFCT.

Interpretability

The global sensitivity indices \(\varrho(\boldsymbol u, S f)\) of the approximation \(S f\) follow directly from its coefficients via Parseval's identity and approximate the indices of \(f\). From them, the attribute ranking \[ r(i) = \frac{1}{C} \sum_{\boldsymbol u \in U,\; i \in \boldsymbol u} \frac{\varrho(\boldsymbol u, S f)}{n_{|\boldsymbol u|}(i)}, \qquad n_s(i) = \big|\{\boldsymbol v \in U : |\boldsymbol v| = s,\ i \in \boldsymbol v\}\big|, \] with a normalisation \(C\) such that \(\sum_{i} r(i) = 1\), measures the importance of the variable \(x_i\). Similarly, the Shapley values of the variance are \[ \mathrm{Sh}(i) = \sum_{\boldsymbol u \in U,\; i \in \boldsymbol u} \frac{\sigma^2(f_{\boldsymbol u})}{|\boldsymbol u|} . \]

This information is used to improve the model: unimportant variables are removed, and the approximation is recomputed on the active set \[ U^{(\varepsilon)} = \{\boldsymbol u \in U_{d_s} : \varrho(\boldsymbol u, S f) > \varepsilon_{|\boldsymbol u|}\} \] with thresholds \(\varepsilon_1, \dots, \varepsilon_{d_s} \in (0,1)\), where larger bandwidths \(N_{\boldsymbol u}\) can be used for the remaining terms. In this way one obtains a sparse and interpretable model and avoids overfitting.

Classification

For binary classification with labels \(y_j \in \{-1, 1\}\), ANOVAapprox computes a classifying function \(f\) with the same ANOVA structure and predicts the class \(\operatorname{sign} f(\boldsymbol x)\). This is a support vector machine in primal form, where the grouped basis functions \(\phi_{\boldsymbol k}\), \(\boldsymbol k \in I(U)\), serve as feature map. The coefficients solve \[ \min_{\hat{\boldsymbol f}} \; \frac{1}{M} \sum_{j=1}^{M} \max\big\{0,\, 1 - y_j \big(\boldsymbol F(\mathcal X, I(U))\, \hat{\boldsymbol f}\big)_j\big\}^2 + \lambda \big\|\hat{\boldsymbol f}\big\|_1 \] with the smooth squared hinge loss. The \(\ell_1\) penalty yields sparse coefficients and thus a clearer selection of the important ANOVA terms. The problem is solved by FISTA, where every iteration needs one multiplication with \(\boldsymbol F(\mathcal X, I(U))\) and its transpose, i.e., one grouped transform. Since \(f\) is again given by its ANOVA terms, the global sensitivity indices and the attribute ranking show which features and feature couplings determine the decision. In both packages, classification is selected by the parameter classification of the constructor.

The method and numerical experiments with the cosine basis and Chui–Wang wavelets, including a variant with \(\ell_2\) regularisation, are presented in the paper Fast and interpretable support vector classification based on the truncated ANOVA decomposition. Already the least squares approach above can be used as a classifier; an example for the adult census data set is given in the paper Grouped transformations and regularization in high-dimensional explainable ANOVA approximation.

Implementation

ANOVAapprox is available as Julia package ANOVAapprox.jl, based on GroupedTransforms.jl and NFFT3.jl, and as Python package pyANOVAapprox, based on pyGroupedTransforms and pyNFFT3. Both are distributed under the GPL 3.0 and are installed with

using Pkg
Pkg.add("ANOVAapprox")

and

pip install pyANOVAapprox

ANOVAapprox.jl was developed by Michael Schmischke and Felix Bartel. Laura Lippert contributed the wavelet bases. Pascal Schröter contributed the mixed basis, maintains ANOVAapprox.jl and developed the Python version pyANOVAapprox.

References

• Potts, D. and Schmischke, M.
  Approximation of high-dimensional periodic functions with Fourier-based methods.
  SIAM J. Numer. Anal. 59, 2393–2429, 2021. doi arXiv pdf
• Potts, D. and Schmischke, M.
  Interpretable approximation of high-dimensional data.
  SIAM J. Math. Data Sci. 3, 1301–1323, 2021. doi arXiv pdf
• Potts, D. and Schmischke, M.
  Learning multivariate functions with low-dimensional structures using polynomial bases.
  J. Comput. Appl. Math. 403, 113821, 2022. doi arXiv pdf
• Bartel, F., Potts, D., and Schmischke, M.
  Grouped transformations and regularization in high-dimensional explainable ANOVA approximation.
  SIAM J. Sci. Comput. 44, A1606–A1631, 2022. doi arXiv pdf
• Potts, D. and Schmischke, M.
  Interpretable transformed ANOVA approximation on the example of the prevention of forest fires.
  Front. Appl. Math. Stat. 8, 795250, 2022. doi arXiv pdf
• Lippert, L., Potts, D., and Ullrich, T.
  Fast hyperbolic wavelet regression meets ANOVA.
  Numer. Math. 154, 155–207, 2023. doi arXiv pdf
• Potts, D. and Schröter, P.
  ANOVA approximation with mixed tensor product basis on scattered points.
  Linear Algebra Appl. 697, 528–560, 2024. doi arXiv pdf
• Akhalaya, K., Nestler, F., and Potts, D.
  Fast and interpretable support vector classification based on the truncated ANOVA decomposition.
  GAMM-Mitt. 48, e202470007, 2025. doi arXiv pdf