# Shapley Values of Structured Additive Regression Models and Application to RKHS Weightings of Functions

Gabriel Dubé, Mario Marchand

<figure><img src="image_1.png" /></figure>

## To cite this version:

Gabriel Dubé, Mario Marchand. Shapley Values of Structured Additive Regression Models and Application to RKHS Weightings of Functions. Transactions on Machine Learning Research Journal, 2025. hal-04753914v3

HAL Id: hal-04753914

https://hal.science/hal-04753914v3

Submitted on 27 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Shapley Values of Structured Additive Regression Models and Application to RKHS Weightings of Functions

**Gabriel Dubé**

Département d'informatique et de génie logiciel
Université Laval

gadub44@ulaval.ca

**Mario Marchand**

Département d'informatique et de génie logiciel
Université Laval

mario.marchand@ift.ulaval.ca

Reviewed on OpenReview: https://openreview.net/forum?id=aWRMuXTvPf

## Abstract

Shapley values are widely used in machine learning to interpret model predictions. However, they have an important drawback in their computational time, which is exponential in the number of variables in the data. Recent work has yielded algorithms that can efficiently and exactly calculate the Shapley values of specific model families, such as Decision Trees and Generalized Additive Models (GAMs). Unfortunately, these model families are fairly restricted. Consequently, we present STAR-SHAP, an algorithm for efficiently calculating the Shapley values of Structured Additive Regression (STAR) models, a generalization of GAMs which allow any number of variable interactions. While the computational cost of STAR-SHAP scales exponentially in the size of these interactions, it is independent of the total number of variables. This allows the interpretation of more complex and flexible models. As long as the variable interactions are moderately-sized, the computation of the Shapley values will be fast, even on high-dimensional datasets. Since STAR models with more than pairwise interactions (e.g. GA2Ms) are seldom used in practice, we also present a new class of STAR models built on the RKHS Weightings of Functions paradigm. More precisely, we introduce a new RKHS Weighting instantiation, and show how to transform it and other RKHS Weightings into STAR models. We therefore introduce a new family of STAR models, as well as the means to interpret their outputs in a timely manner.

## 1 Introduction

The ability to interpret a machine learning model is an increasingly important consideration. Highly complex models can achieve great accuracy, but understanding their inner workings is challenging. Why did the model offer this particular prediction for this input? The decision was based on which features of the data? On the other hand, simpler models, like Generalized Additive Models (GAMs<fnref n="1" />, Hastie (2017)) or Decision Trees (Kotsiantis, 2013), are inherently easier to understand, but do not always reach the same levels of performance (though recent work appears to show otherwise (Rudin et al., 2024)). Work to improve explanation methods, especially the interpretation of more complex models, can help bridge this gap.

Informally, we mean by *interpreting a model* the task of explaining its output. In sensitive contexts such as healthcare (Caruana et al., 2015), energy systems (Zhang & Chen, 2024) or machine fault diagnosis (Brusa et al., 2023), to give only a few examples, mistakes can be costly, and so the ability to understand, and thus trust, an AI model is crucial. In this paper, we specifically consider the **attribution problem**, which consists in distributing the prediction score between the data features (Sundararajan & Najmi, 2020). On one

<footnote n="1">¹Table 1 lists all the acronyms used in this paper.</footnote>
end of the spectrum of interpretability lie linear models, where the contribution of each variable is obvious and independent of all other variables. On the other end, we have large deep neural networks, the incredibly complex inner workings of which are singularly difficult to elucidate (Montavon et al., 2018).

Post-hoc interpretability methods attempt to explain even complex models through various means. Just some examples are LIME (Ribeiro et al., 2016), which explains the prediction of a complex model by approximating the model locally around the instance using a simple predictor, and Partial Dependence Plots (Friedman, 2001), which visualize the effect of a variable by holding others constant.

Another widely used tool for interpreting models, and the subject of this paper, is the Shapley value (see e.g. Winter (2002)), originally devised in the context of game theory (Shapley et al., 1953). Adapted to the machine or statistical learning context, it can be used to assess the contribution of each variable to the output of a model, for a given example. Additionally, it can be averaged over a dataset to see the relative importance of each variable.

With all of this in mind, the purpose of this paper is twofold. First, we present an algorithm for efficiently calculating the Shapley values of a large family of prediction models, Structured Additive Regression (STAR) models. STAR models are a further generalization of GAMs that allow any number of variable interactions in each term, making them potentially much more expressive models. The algorithm that we present does have an exponential complexity in the size of these interactions, but has no explicit dependency on the total number of variables. This allows interpreting models with many more variable interactions, helping to make complex models more viable in sensitive domains. We show that our algorithm is especially useful in the context of high-dimensional datasets, such as Wood et al. (2015), because of this lack of an explicit dependency on the total number of variables.

Second, we improve the understanding and usability of a newly introduced class of models, the RKHS Weightings of Functions (Dubé & Marchand, 2024), mainly in the context of STAR models. To do so, we present a new RKHS Weighting instantiation (Dubé & Marchand, 2024) which is better suited to regression. Then, we show how to transform it and other RKHS Weightings into STAR models. This yields an entirely new family of STAR models. We hope that these developments will foster interest and further research into RKHS Weightings, which we believe have significant untapped potential.

The rest of the paper is as follows. Section 2 covers Shapley values and RKHS Weightings of Functions. Section 3 presents the formula and corresponding algorithm for calculating the Shapley values of a STAR model. Section 4 presents a new RKHS Weighting instantiation, and shows how to transform an RKHS Weighting model into a STAR model. Section 5 contains various experiments, showing the computation time of Shapley values using our new algorithm, the regression and time series performance of the models introduced in this paper, and an example of interpreting a model using our new algorithm. Finally, in Section 6, we go over some limitations of this paper as well as future work to be done, and conclude in Section 7.

# 2 Review of Shapley values and RKHS Weightings

The two main subjects of this paper are Shapley values and RKHS Weightings, two tools we use within the context of machine learning. We start by defining the basic machine learning notions that we use in this paper. Then, we present our definition of the Shapley value, and end on a summary of RKHS Weightings.

## 2.1 Machine learning

The regression setting in machine learning consists of learning to predict real values from examples. Given an instance space $\mathcal{X}$, usually $\mathcal{X} = \mathbb{R}^n$, the label space $\mathcal{Y} = \mathbb{R}$, and a dataset of examples $S = \{(x_i, y_i)\}_{i=1}^m \subset \mathcal{X} \times \mathcal{Y}$, our task is to find a predictor (also called model) $h: \mathcal{X} \to \mathcal{Y}$ that optimizes some objective function of the sample data. To measure the quality of a model, we use a loss function $\ell: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}$, which in this paper is the squared loss:

$$ \ell(h(x), y) := \frac{1}{2}(h(x) - y)^2. \qquad (1) $$
<table><caption>Table 1: Acronyms used in this paper.</caption><thead><tr><th>Acronym</th><th>Meaning</th><th>Relevant link</th></tr></thead><tbody><tr><td>GAM</td><td>Generalized Additive Model</td><td>Hastie (2017)</td></tr><tr><td>GA2M</td><td>GAM with pairwise interactions</td><td>Lou et al. (2013)</td></tr><tr><td>STAR (model)</td><td>Structured Additive Regression (model)</td><td>Equation (16)</td></tr><tr><td>k-STAR model</td><td>STAR model using only k-way interactions</td><td></td></tr><tr><td>RKHS</td><td>Reproducing Kernel Hilbert Space</td><td>Berlinet & Thomas-Agnan (2011)</td></tr><tr><td>MSE</td><td>Mean Square Error</td><td>Equation (2)</td></tr><tr><td>R<sup>2</sup></td><td>Coefficient of determination</td><td>Equation (3)</td></tr><tr><td>SHAP</td><td>SHapley Additive exPlanations</td><td>Lundberg & Lee (2017)</td></tr><tr><td>SHAP-IQ</td><td>SHAPley Interaction Quantification</td><td>Fumagalli et al. (2023)</td></tr><tr><td>STAR-SHAP</td><td>Structured Additive Regression SHAP</td><td>Algorithm 1</td></tr><tr><td>erf</td><td>Error function, erf(x) := 2/√π ∫<sup>x</sup><sub>0</sub> e<sup>-t<sup>2</sup></sup>dt</td><td></td></tr><tr><td>ReLU</td><td>Rectified Linear Unit, ReLU(z) := max(0, z)</td><td></td></tr><tr><td>RWSign</td><td>RKHS Weighting using sign function as base predictor</td><td>I1 in Dubé & Marchand (2024)</td></tr><tr><td>RWStumps</td><td>RKHS Weighting using decisions stumps as base predictor</td><td>I2 in Dubé & Marchand (2024)</td></tr><tr><td>RWReLU</td><td>RKHS Weighting using ReLU as base predictor</td><td>Section 4.1</td></tr><tr><td>EBM</td><td>Explainable Boosting Machine</td><td>Lou et al. (2013)</td></tr><tr><td>SVR</td><td>Support Vector Regression</td><td>Pedregosa et al. (2011)</td></tr></tbody></table>

The empirical risk is the sample average loss. When using the squared loss, this is called the Mean Square Error (MSE):

$$ \text{MSE}_{\mathcal{S}}(h) := \frac{1}{2m} \sum_{i=1}^{m} (h(x_i) - y_i)^2. \qquad (2) $$

The MSE has units of the labels squared. To successfully compare experimental results on different datasets, we instead report the coefficient of determination:

$$ R_{\mathcal{S}}^{2}(h) := 1 - \frac{\sum_{i=1}^{m} (y_i - h(x_i))^2}{\sum_{i=1}^{m} (y_i - \bar{y})^2}, \qquad (3) $$

where $\bar{y} := \frac{1}{m} \sum_{i=1}^{m} y_i$ is the average label.

In practice, simply minimizing the sample MSE is insufficient. We seek a model that will successfully generalize to new, yet-to-be-seen instances. Indeed, we can assume that the data comes from some unknown distribution $\mathcal{D}$ over $\mathcal{X} \times \mathcal{Y}$. The true risk, or true MSE, is the expected squared loss over this distribution:

$$ \text{MSE}_{\mathcal{D}}(h) := \frac{1}{2} \mathbb{E}_{(x,y) \sim \mathcal{D}} [(h(x) - y)^2]. \qquad (4) $$

We estimate this value by calculating the MSE on a set of examples distinct from the sample.

A learning algorithm $\mathcal{A}$ is a function that takes in a sample $\mathcal{S}$, and returns a predictor, i.e. $\mathcal{A}(\mathcal{S}) = h(x)$. It is this algorithm that seeks to minimize the true risk by optimizing a function of the best available proxy, the empirical MSE. Most learning algorithms have hyperparameters; the number of trees in a random forest; the learning rate of a neural network; the regularization parameter of the kernel ridge regression. Each combination of parameters yields an algorithm. Suppose we have $C$ combinations to choose from. This implies we have the candidate algorithms $A_1$ to $A_C$. To select the best one, we use cross-validation, more specifically $k$-fold cross-validation. The sample $\mathcal{S}$ is first divided into $k$ equal disjoint subsets $S_1$ to $S_k$. Each algorithm $A_c$ is applied $k$ times, giving us the $k$ models $A_c(\mathcal{S} \setminus S_1), \dots, A_c(\mathcal{S} \setminus S_k)$. We can estimate the true risk of each of these models through its MSE on the holdout set, e.g. $\text{MSE}_{\mathcal{D}}(A_c(\mathcal{S} \setminus S_1)) \approx \text{MSE}_{S_1}(A_c(\mathcal{S} \setminus S_1))$. In turn, we can estimate the true risk of $A_c(\mathcal{S})$ (the model trained on the entire dataset) as the average MSE of these $k$ models:

$$ \text{MSE}_{\mathcal{D}}(A_c(\mathcal{S})) \approx \frac{1}{k} \sum_{i=1}^{k} \text{MSE}_{S_i}(A_c(\mathcal{S} \setminus S_i)). \qquad (5) $$
The algorithm with the lowest approximated MSE is chosen, and trained again on the entire sample.

## 2.2 Shapley values

The Shapley value was originally derived by Shapley et al. (1953) as a uniquely interesting attribution method in the context of a game with $n$ players. Denoting $[n] := \{1, \dots, n\}$ the set of natural numbers up to $n$, and $\mathcal{P}([n])$ its power set, the Shapley value revolves around a **value function** $\nu: \mathcal{P}([n]) \to \mathbb{R}$. Given a coalition of players $I \subseteq [n]$, $\nu$ is the payoff that this coalition can achieve in the game by playing together. The Shapley value of player $i$ is defined as:

$$ \phi_i(\nu) = \sum_{I \subseteq [n] \setminus \{i\}} \frac{|I|!(n - |I| - 1)!}{n!} (\nu(I \cup \{i\}) - \nu(I)). \quad (6) $$

It is the average value added by player $i$ when it joins a coalition where $i$ is initially absent. The Shapley value is the unique attribution method that satisfies four key properties (Fryer et al., 2021):

1. **Efficiency.** The Shapley values of all players sum to the value of all players together:

$$ \sum_{i=1}^{n} \phi_i(\nu) = \nu([n]). \quad (7) $$

In simpler words, the gain is distributed among all the agents/players.

2. **Symmetry.** Two players with equal contributions have equal Shapley values:

$$ (\forall I \subseteq [n] \setminus \{i, j\}, \nu(I \cup \{i\}) = \nu(I \cup \{j\})) \implies \phi_i(\nu) = \phi_j(\nu). \quad (8) $$

3. **Null player.** A player that has no contribution has a Shapley value of 0:

$$ (\forall I \subseteq [n] \setminus \{i\}, \nu(I \cup \{i\}) = \nu(I)) \implies \phi_i(\nu) = 0. \quad (9) $$

4. **Linearity.** The Shapley values of a sum $\nu_1 + \nu_2$ of value functions are the sum of the Shapley values of the individual functions:

$$ \forall i \in [n], \phi_i(\nu_1 + \nu_2) = \phi_i(\nu_1) + \phi_i(\nu_2). \quad (10) $$

These properties make the Shapley value a uniquely interesting attribution method for the contribution of each variable to the prediction of a model. Although some raise doubts on its reliability (Huang & Marques-Silva, 2023) or the desirability of these four properties (Fryer et al., 2021), we do not address this debate in this paper.

Different value functions yield different Shapley values, so the choice of value function is crucial. In this paper, we wish to attribute to each variable in a dataset its contribution to the output of a model. Consider the instance space $\mathcal{X} = \mathbb{R}^n$, so that the natural numbers up to $n$, which we denoted $[n] := \{1, \dots, n\}$, identify the $n$ variables that make up an instance $x = (x_1, \dots, x_n) \in \mathcal{X}$. For some model $h: \mathcal{X} \to \mathbb{R}$, what is the contribution of variable $i \in [n]$ to the model output $h(x)$? In general, it is not possible to see what a model would output given a coalition $I \subset [n]$ of variables, since most models are undefined with missing values. To sidestep this issue, the solution most commonly used (Lundberg & Lee, 2017) is to intervene on the values of $x$ on the variables outside of the coalition. To that end, we define the **replacement function**. For any feature subset $P \subseteq [n]$, and $x, z \in \mathcal{X}$, the vector $r_P(z, x)$ is $z$, except for the variables in $P$ that are replaced by the values in $x$:

$$ r_P(z, x)_i = \begin{cases} z_i & i \notin P \\ x_i & i \in P. \end{cases} \quad (11) $$
Finally, suppose that we have an instance $x$ of interest and a **background** $S \subseteq \mathcal{X}$ consisting of a set of instances. Denote by $U(S)$ the uniform distribution over $S$. The value function $\nu: \mathcal{P}([n]) \to \mathbb{R}$ that we focus on is:

$$ \nu(I) := \mathbb{E}_{z \sim U(S)} [h(r_I(z, x))]. \qquad (12) $$

In simple words, the output of $h$ is averaged over the background, but the variables in $I$ are fixed to their value in $x$. This value function yields a so-called **interventional** Shapley value, after the act of replacing feature values with those from another instance, i.e. “intervening” on these features. This is the Shapley value calculated by SHAP (Lundberg & Lee, 2017).

Rather than Equation (6), we use an equivalent form using permutations of $[n]$ (see e.g. Equation (3) of Khorrami Chokami & Rabitti (2024)). A **permutation** is a bijective function $\pi: [n] \to [n]$; $\pi(i)$ is the position of variable $i$ in the permutation. Let $\Omega$ be the set of all permutations of $[n]$, and $U(\Omega)$ be the uniform distribution over $\Omega$. Also define:

$$ \pi_{:i} := \{j \in [n] \mid \pi(j) < \pi(i)\}, \qquad (13) $$

the set of variables preceding (but excluding $i$) in the order implied by the permutation $\pi$. Given all of these notations, the **Shapley value** of variable $i$, for example $x$, for the predictor $h$, is defined as:

$$ \phi_i^{\text{SHAP}}(h, x) := \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} [h(r_{\pi_{:i} \cup \{i\}}(z, x)) - h(r_{\pi_{:i}}(z, x))]. \qquad (14) $$

(Definition taken from Laberge et al. (2023).)

Calculating Equation (14) directly is in $O(n!)$, as it requires enumerating every permutation of $[n]$. It therefore cannot reasonably be used except on particularly low-dimensional datasets. However, some model families admit much faster formulas or algorithms for calculating the Shapley values. Such model families include linear models (Lundberg & Lee, 2017), Generalized Additive Models (GAMs), and Trees and ensembles of Trees (Lundberg et al., 2018). The common characteristic of these models is limited variable interactions, so that the contribution of a particular variable can be calculated independently of the variables with which it does not interact. For example, the Shapley values of a linear model $h(x) = \sum_j w_j x_j$ are obtained easily from its coefficients by the formula $\phi_i^{\text{SHAP}}(h, x) = w_i(x_i - \mathbb{E}_{z \sim U(S)} z_i)$ (Lundberg & Lee, 2017); the Shapley values of a GAM $h(x) = \sum_j f_j(x_j)$ are given by:

$$ \phi_i^{\text{SHAP}}(h, x) = f_i(x_i) - \mathbb{E}_{z \sim U(S)} [f_i(z_i)]. \qquad (15) $$

(The proof is in the appendix.)

Structured Additive Regression (STAR) models are a generalization of both linear models and GAMs (see e.g. Mayer et al. (2021)). A STAR model allows variable interactions by defining the predictor over subsets of the data features:

$$ h(x) = \sum_{I \subseteq [n]} f_I(x_I), \qquad (16) $$

where $x_I$ is the partial vector $(x_{I_1}, \dots, x_{I_{|I|}})$, and $f_I$ is a function defined over the feature subset $I$. Importantly, not all of these functions $f_I$ need to be nonzero. We can indeed choose which variable interactions are desired by selecting which feature subsets $I$ are associated with a nonzero $f_I$. For example, GAMs use only singleton feature subsets, i.e. individual variables, while GA2Ms (Lou et al., 2013) use up to pairwise interactions between the variables.

STAR models have various names in the literature. Bordt & von Luxburg (2023) define Generalized Additive Models of order $k$ as models of the form:

$$ h(x) = \sum_{I \subseteq [n], |I| \le k} f_I(x_I), \qquad (17) $$

and introduce $k$-Shapley Values, which attribute a contribution to each feature subset, instead of only individual variables. This generalizes Shapley values, but suffers from the same high computational cost.
Laberge et al. (2024) call Equation (16) a Functional Decomposition, following previous work on the subject (Hooker, 2004; Kuo et al., 2010), and show among other things how to obtain the Shapley values of a function from its Interventional Decomposition (Kuo et al., 2010). This, again, is a high cost procedure.

In short, there are two types of formulas or algorithms for exactly calculating the Shapley values of a model. Model-agnostic formulas like the original Equation (13) of Shapley et al. (1953), Equation 14, Equation (11) of Bordt & von Luxburg (2023); and algorithms like SHAP Lundberg & Lee (2017). These algorithms all suffer from the same delimitating complexity with regard to the total number of variables. On the other hand, model-specific formulas, such as those previously mentioned for linear models and GAMs, and algorithms, like TreeSHAP (Lundberg et al., 2018), are vastly more efficient, but are constrained to their respective model families.

In order to work with higher dimensional data and arbitrary models, model-agnostic approximation algorithms have been developed. Notably among these are KernelSHAP (Lundberg & Lee, 2017), Unbiased KernelSHAP (Covert & Lee, 2020) and SHAP-IQ (Fumagalli et al., 2023). These provide approximative Shapley values (and even k-Shapley values in the case of SHAP-IQ) without the exponential complexity in the number of variables. Although the values are approximations, rather than the exact Shapley values, Unbiased KernelSHAP and SHAP-IQ have been shown to be unbiased, i.e. they return values which, in expectation, are the correct ones, and consistent, meaning that the variance of the approximations goes to zero at the limit as the Monte Carlo sampling size increases. Approximation is therefore a valid alternative when efficient model-specific exact formulas or algorithms are not known or do not exist.

In Section 3, we present STAR-SHAP, a model-specific algorithm for calculating the exact Shapley values of any STAR model. As long as the maximal number k of interactions is limited, this is an efficient algorithm. Specifically, its computational complexity does scale exponentially in the maximal size k of the feature subsets and scales linearly in the number of subsets I used in Equation (16), but does not explicitly depend on the dimensionality n of the data. This allows using many more interaction terms, while maintaining the ability to interpret the model predictions through its Shapley values.

In Section 5.2, we compare the computation time of this algorithm to the model-agnostic methods we previously mentioned: SHAP, KernelSHAP, Unbiased KernelSHAP, and SHAP-IQ. We also study the quality of these approximations in relation to the exact values provided by our algorithm.

As a final note to this section, we wish to mention the recent work of Khorrami Chokami & Rabitti (2024), who provides an exact formula for the Shapley values of GAMs, however using a different value function. Where we seek to explain the output of a model on specific instances, Khorrami Chokami & Rabitti (2024) are interested in the variance explained by each variable across the entire data distribution. They therefore use the value function $\nu(I) := \text{Var}_x [\sum_{i \in I} f_i(x_i)]$, which yields, as the corresponding Shapley value, a global variable importance index, as opposed to a granular instance-specific explanation of the output of the model. This shows how rich the space of possible Shapley values are when considering all the possible value functions that can be used. This thought aside, we focus on the value function given by Equation (12) for the remainder of this paper.

## 2.3 RKHS Weightings

Dubé & Marchand (2024) introduced a new and flexible model family, parameterized through the tuple $(W, \phi, K, p)$ of a parameter space $W$; base predictor $\phi: W \times X \to \mathbb{R}$; kernel $K: W \times W \to \mathbb{R}$; distribution $p$ over $W$. The kernel $K$ implicitly defines a reproducing kernel Hilbert space (RKHS) $\mathcal{H}$, a space of functions of the form $\alpha: W \to \mathbb{R}$ with particular properties. (A thorough understanding of reproducing kernel Hilbert spaces is not required for this paper, but the interested reader can see Dubé & Marchand (2024) for a short primer on RKHS.) They define the operator $\Lambda: \mathcal{H} \to L^\infty(\mathcal{X})$, which takes as input a weight function $\alpha \in \mathcal{H}$, and returns the bounded predictor $\Lambda\alpha: \mathcal{X} \to \mathbb{R}$ via the equation:

$$ \Lambda\alpha(x) := \mathbb{E}_{w \sim p} [\alpha(w)\phi(w, x)]. \qquad (18) $$
In practice, the weight function $\alpha \in \mathcal{H}$ takes the form of a sum $\sum_{t=1}^{T} a_t K(w_t, \cdot)$ (a property of the RKHS), for some real coefficients $a_1, \dots, a_T$ and $w_1, \dots, w_T \in \mathcal{W}$. The model can then be written as:

$$ \Lambda\alpha(x) = \sum_{t=1}^{T} a_t \mathbb{E}_{w \sim p} [K(w_t, w)\phi(w, x)]. \qquad (19) $$

Therefore, an **RKHS Weighting** model is a weighted sum of the partial functions defined by the expectations $E_{w \sim p}[K(w_t, w)\phi(w, x)]$. Each such function of $x$ is parameterized by a $w_t \in \mathcal{W}$. While this is in theory a highly flexible family of models, the current bottleneck for its utilization is the requirement for exactly calculating these expectations. This can be done, depending on the choice of instantiation $(\mathcal{W}, \phi, \mathcal{K}, p)$. Dubé & Marchand (2024) presented two instantiations of the model, which they called I1 and I2. We will refer to these as RWSign and RWStumps respectively after their base predictor $\phi$, a simple sign function for RWSign, and a decision stump for RWStumps, with RW standing for RKHS Weighting. The reader can refer to Table 3 in Section 4 for the explicit definition of RWSign and RWStumps, and to Table 1 of Dubé & Marchand (2024) for the analytical form of the expectation $E_{w \sim p}[K(w_t, w)\phi(w, x)]$. In both cases, the entire model (Equation (19)) becomes a weighted sum of error functions<fnref n="2" />, a highly expressive form.

Since Equation (18) represents an aggregation of base predictor predictions (with regard to a distribution $p$ and a weight function $\alpha$), and since the base predictor, in both cases, can only take two values, we surmise that while these instantiations can be used for regression, they are first of all intended for classification. We introduce in Section 4 a new instantiation which we expect to be better suited to regression, using the **rectified linear unit**<fnref n="3" /> (ReLU) as the base predictor. Furthermore, we show in Section 4 how to generally build RKHS Weightings that are also STAR models. We also show specifically how some existing instantiations, which satisfy some conditions, can easily be converted into STAR models. This includes RWSign and RWRelu, the instantiation that we introduce in Section 4. (RWStumps is in fact already a GAM, and so does not require any modification.) These models are all interpretable, as their Shapley values can be calculated using Algorithm 1, which we introduce in the next section.

## 3 Shapley values of STAR models

The following theorem encapsulates the main contribution of this paper, an efficient formula for calculating the Shapley values of Structured Additive Regression (STAR) models.

**Theorem 3.1** (Shapley values of a STAR model). Consider a STAR model:

$$ h(x) = \sum_{I \subseteq [n]} f_I(x_I). \qquad (20) $$

Consider the replacement function $r_P$ defined in Equation (11), and the sets:

$$ \mathcal{A}^{+}(i, I) := \{A \subseteq I | i \in A, 1 \le |A| < |I|\} \qquad (21) $$

$$ \mathcal{A}^{-}(i, I) := \{A \subseteq I | i \in A, 1 < |A| \le |I|\}. \qquad (22) $$

Then the Shapley values of the model are:

$$ \phi_i^{\text{SHAP}}(h, x) = \mathbb{E}_{z \sim \mathcal{U}(S)} \left[ \sum_{I: f_I \neq 0, i \in I} \left( \frac{f_I(x_I) - f_I(z_I)}{|I|} + \sum_{A \in \mathcal{A}^+(i, I)} \frac{f_I(r_A(z, x)_I)}{|A|\binom{|I|}{|A|}} - \sum_{A \in \mathcal{A}^-(i, I)} \frac{f_I(r_A(z, x)_I)}{|A|\binom{|I|}{|A|}} \right) \right]. \qquad (23) $$

*Proof.* The proof is a series of combinatorial arguments, and can be found in the appendix. □

<footnote n="2">²erf($x$) := $\frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt$</footnote>

<footnote n="3">³ReLU($z$) := $\max(0, z)$</footnote>
Assuming that the feature subsets $I$ are limited in size, then Equation (23) is vastly more efficient to calculate than the definition of the Shapley values in Equation (14). Instead of going through all the permutations of $[n]$, Equation (23) requires only enumerating the subsets of each feature subset $I$ present in Equation (16). However, we should note that using Equation (23) as it is to calculate all the Shapley values for a given instance leads to many redundant computations. Indeed, observe that:

1. For any $j \in A \in \mathcal{A}^+(i, I)$, we also have $A \in \mathcal{A}^+(j, I)$. This means that the computation $f_I(r_A(z, x)_I)$ appears in the calculation of the Shapley value of each variable in $A$.

2. Consider any $A \in \mathcal{A}^+(i, I)$, and a variable $j \in I \setminus A$. Then $A \cup \{j\} \in \mathcal{A}^-(j, I)$, and so the computation $f_I(r_{(A \cup \{j\}) \setminus \{j\}}(z, x)_I)$ is equal to $f_I(r_A(z, x)_I)$.

In short, when calculating the Shapley values of all variables (rather than a single one), the same computation $f_I(r_A(z, x)_I)$ is done exactly $|I|$ times; $|A|$ times from the first observation, and $|I| - |A|$ times from the second observation. We can save a factor $|I|$ by doing this computation only once. The term $\frac{f_I(x_I) - f_I(z_I)}{|I|}$ can also be calculated once, as the value is the same for all variables. These optimizations in mind, we present Algorithm 1, which simultaneously calculates all the Shapley values of a STAR model for a given instance.

<caption>**Algorithm 1 STAR-SHAP: Shapley Values of a Structured Additive Regression Model**</caption>

**input** Input space $\mathcal{X} \subseteq \mathbb{R}^n$, instance $x \in \mathcal{X}$, STAR model $h(x) = \sum_{I \subseteq [n]} f_I(x_I)$, sample $S \in \mathcal{X}^m$
**Initialize** $\Phi = (\phi_1, \dots, \phi_n) = 0 \in \mathbb{R}^n$
**for** $z \in S$, $I \subseteq [n]$ such that $f_I \neq 0$ **do**
    **baseline** $\leftarrow \frac{f_I(x_I) - f_I(z_I)}{|I|}$
    **for** $i \in I$ **do**
        $\phi_i \leftarrow \phi_i + \text{baseline}$
    **end for**
    **for** $A \subseteq I$, $A \neq \emptyset$ **do**
        term $\leftarrow f_I(r_A(z, x)_I)$
        **for** $i \in I$ **do**
            **if** $i \in A$ and $A \in \mathcal{A}^+(i, I)$ **then**
                $\phi_i \leftarrow \phi_i + \frac{\text{term}}{|A|(|A|)}$
            **else if** $i \in I \setminus A$ and $(A \cup \{i\}) \in \mathcal{A}^-(i, I)$ **then**
                $\phi_i \leftarrow \phi_i - \frac{\text{term}}{(|A|+1)(|A|+1)}$
            **end if**
            **end for**
            **end for**
            **end for**
        **output** $\frac{1}{m}\Phi$

A highly interesting characteristic of Algorithm 1 is the lack of an explicit dependency on the dimensionality of the data in its algorithmic complexity. Indeed, the dominant operation is the computation of the partial model outputs $f_I(r_A(z, x)_I)$. Note that, at this point, we make no assumption on the partial models $f_I$ other than that they utilize $|I|$ variables. Let us therefore suppose that some nondecreasing function $F: \mathbb{N} \to [0, \infty)$ characterizes the computational complexity of the partial models, so that the cost of calculating $f_I(x_I)$ is in $O(F(|I|))$ for all $I$ (and all $x$). Denote $k$ the maximal size of the feature subsets associated with a nonzero partial model ($f_I \neq 0$), and $N$ the number of such feature subsets. Algorithm 1 is a loop over these $N$ feature subsets and over the $m$ examples in the background $S$. For each feature subset $I$, the algorithm then iterates over all its subsets $A \subseteq I$, calculating the partial output $f_I(r_A(z, x)_I)$ each time. Since there are $2^k$ subsets to a set of size $k$, this dominant operation is done a total of $mN2^k$ times, each with a cost of $O(F(|I|))$. Since $|I| \le k$, this leads to a favorable algorithmic complexity of $O(mN2^kF(k))$, independently of the dimensionality $n$ of the data. One could argue that $N$, the number of feature subsets $I$ defining the
<table><caption>Table 2: Computational complexity of various Shapley value algorithms for calculating all *n* Shapley values of a single *n*-dimensional instance. All algorithms: *m* is the size of the background sample. TreeSHAP: *T* is the number of trees, *L* is the number of leaves in a tree, and *D* is the depth of a tree. STAR-SHAP: *N* is the number of nonzero terms in the model, *k* is the maximal number of interactions per term, *F*(*k*) is the cost of computing a term with *k* interactions. SHAP: *n* is the number of variables, *F*(*n*) is the cost of computing the model. All others: *K* is the sampling size for the approximation. The reported complexities of KernelSHAP, Unbiased KernelSHAP and SHAP-IQ are those of their implementation in Muschalik et al. (2024), version 1.1.0, the software package used for our experiments in Section 5. Note that their implementations of Unbiased KernelSHAP and SHAP-IQ coincide for calculating 1-Shapley values, the use case for this paper. This is faster than the original Unbiased KernelSHAP implementation, while yielding the same values (Fumagalli et al., 2023).</caption><thead><tr><th>Algorithm</th><th>Target models</th><th>Exact?</th><th>Complexity</th></tr></thead><tbody><tr><td>TreeSHAP</td><td>Trees, ensembles of trees</td><td>Yes</td><td>O(mTLD<sup>2</sup>)</td></tr><tr><td>STAR-SHAP</td><td>STAR models</td><td>Yes</td><td>O(mN<sup>2</sup>kF(k))</td></tr><tr><td>SHAP</td><td>Any</td><td>Yes</td><td>O(m2<sup>n</sup>F(n))</td></tr><tr><td>KernelSHAP</td><td>Any</td><td>No</td><td>O(Kn<sup>2</sup> + n<sup>3</sup> + KmF(n))</td></tr><tr><td>Unbiased KernelSHAP</td><td>Any</td><td>No</td><td>O(Kn<sup>2</sup> + KmF(n))</td></tr><tr><td>SHAP-IQ</td><td>Any</td><td>No</td><td>O(Kn<sup>2</sup> + KmF(n))</td></tr></tbody></table>

model in Equation 16, can grow exponentially in the dimension *n* of *X*. However, *N* is in reality chosen by whoever builds the model. Rather than using all possible subsets, a better strategy is to use prior knowledge about the problem at hand to identify meaningful interactions. Examples include using adjacent or nearby pixels in images, or small sequences of observations in time series. In this way, the computational cost of Algorithm 1 can easily be controlled, and even be held constant as the dimensionality of the data increases. The experiments of Section 5.2 illustrate this.

## 3.1 Comparison to other Shapley value algorithms

Table 2 presents the computational complexity of the various Shapley value algorithms mentioned in this paper. As we see, SHAP has an exponential complexity in the dimensionality of the data (the number of variables), and all other algorithms must make concessions to escape this prohibitive computational cost. One way is to assume a specific form for the model (TreeSHAP, STAR-SHAP). The other is to give up the exactness of the values returned by the algorithm. The most salient difference between STAR-SHAP and the model-agnostic algorithms is that STAR-SHAP only needs to calculate the output of the partial models (*f*<sub>*t*</sub> in Equation (16)), entirely forgoing an explicit dependency in the dimensionality of the data. When this cost *F*(*k*) is much smaller than the cost *F*(*n*) of calculating the output of the entire model, for instance in large, high-dimensional ones, we should expect to see a notably faster runtime for STAR-SHAP. Combined with the exactness of the values returned by STAR-SHAP, we believe it to be the algorithm to use in the context of STAR models. We run experiments in Section 5.2 to demonstrate these advantages of STAR-SHAP over other algorithms.

## 3.2 Applications of STAR-SHAP

STAR-SHAP can calculate the Shapley values, as defined in Equation (23), for any STAR model. This includes Generalized Additive Models (GAMs) (Hastie, 2017), Generalized Additive Models with Pairwise Interactions (GA2Ms) (Lou et al., 2013), Sparse Polynomial Regression (Huang et al., 2010), Sparse Decision Trees (Lin et al., 2020), and others. GAMs and Trees already have dedicated algorithms for computing their Shapley values (see Section 2.2), but for the others, Algorithm 1 can do so as well, as long as the number of variable interactions is moderate. Finally, another family of STAR models consists of STAR RKHS Weightings, which we introduce in the next section.
<table><caption>Table 3: Three RKHS Weightings instantiations.</caption><thead><tr><th></th><th>RWSign</th><th>RWRelu</th><th>RWStumps</th></tr></thead><tbody><tr><td>W</td><td>R<sup>n</sup></td><td>R<sup>n</sup></td><td>{1, ..., n} × R</td></tr><tr><td>φ(w, x)</td><td>sign(⟨w, x⟩)</td><td>max(0, ⟨w, x⟩)</td><td>sign(x<sub>w1</sub> − w<sub>2</sub>)</td></tr><tr><td>K(u, w)</td><td>exp(−‖u − w‖<sup>2</sup> / (2γ<sup>2</sup>))</td><td>exp(−‖u − w‖<sup>2</sup> / (2γ<sup>2</sup>))</td><td>1[u<sub>1</sub> = w<sub>1</sub>] exp(−(u<sub>2</sub> − w<sub>2</sub>)<sup>2</sup> / (2γ<sup>2</sup>))</td></tr><tr><td>p</td><td>N(0, σ<sup>2</sup>I<sub>n</sub>)</td><td>N(0, σ<sup>2</sup>I<sub>n</sub>)</td><td>U({1, ..., n}) × N(0, σ<sup>2</sup>)</td></tr></tbody></table>

# 4 Structured Additive RKHS Weightings of Functions

In this section, we introduce a new RKHS Weighting instantiation better suited to regression than those found in Dubé & Marchand (2024), and show how to transform any generic instantiation of a certain form into a Structured Additive Regression (STAR) model.

## 4.1 ReLU instantiation

An RKHS Weighting instantiation is a tuple $(\mathcal{W}, \phi, \mathcal{K}, p)$ of a parameter space $\mathcal{W}$, base predictor $\phi$, kernel $\mathcal{K}$ (implying an RKHS $\mathcal{H}$) and distribution $p$. For example, instantiation I1 of Dubé & Marchand (2024) (which we call RWSign) uses $\mathcal{W} = \mathcal{X} = \mathbb{R}^n$, $\phi(w, x) = \text{sign}(\langle w, x \rangle)$, the gaussian kernel $\mathcal{K}(u, w) = \exp(-\frac{\|u-w\|^2}{2\gamma^2})$, and a normal distribution $p = \mathcal{N}(0, \sigma^2 I_n)$.

The instantiation that we propose, which we name RWRelu, consists simply in replacing the base predictor of RWSign by the rectified linear unit (ReLU), i.e. $\phi(w, x) = \text{ReLU}(\langle w, x \rangle) = \max(0, \langle w, x \rangle)$. (See Table 3 for a quick overview of each instantiation.) In order to calculate the output of the model, which takes the form:

$$ \Lambda\alpha(x) = \sum_{t=1}^{T} a_t \mathbb{E}_{w \sim p} [\mathcal{K}(w_t, w)\phi(w, x)], \quad (24) $$

for some $\alpha = \sum_{t=1}^{T} a_t \mathcal{K}(w_t, \cdot) \in \mathcal{H}$, we need an analytical form for the expectation in Equation (24). We encapsulate that formula in the following theorem.

**Theorem 4.1.** Consider the RKHS Weighting instantiation $(\mathcal{W}, \phi, \mathcal{K}, p)$ where $\mathcal{W} = \mathcal{X} = \mathbb{R}^n$, $\phi(w, x) = \max(0, \langle w, x \rangle)$, $\mathcal{K}(u, w) = \exp(-\frac{\|u-w\|^2}{2\gamma^2})$, $p = \mathcal{N}(0, \sigma^2 I_n)$. Then we have:

$$ \mathbb{E}_{w \sim p} [\mathcal{K}(u, w)\phi(w, x)] = \left(1 + \frac{\sigma^2}{\gamma^2}\right)^{-n/2} e^{\frac{-\|u\|^2}{2\sigma^2 + 2\gamma^2}} \frac{\|x\|}{2\sqrt{\pi}} \left( \sqrt{2\zeta} e^{-\frac{\langle u', x \rangle^2}{2\zeta^2 \|x\|^2}} + \sqrt{\pi} \frac{\langle u', x \rangle}{\|x\|} \left[ 1 + \text{erf} \left( \frac{\langle u', x \rangle}{\sqrt{2\zeta} \|x\|} \right) \right] \right), \quad (25) $$

where $\zeta$ is defined by the relationship $\frac{1}{2\zeta^2} = \frac{1}{2\gamma^2} + \frac{1}{2\sigma^2}$ and $u' := (1 + \frac{\gamma^2}{\sigma^2})^{-1} u$.

*Proof.* The calculus can be found in the appendix. □

Dubé & Marchand (2024) also define some instantiation-dependent theoretical constants (see Theorem 4 of Dubé & Marchand (2024)) which are then used in the theoretical guarantees. We provide those values for RWRelu and the relevant calculus in the appendix. Importantly, this means that all theoretical guarantees proved in Dubé & Marchand (2024) also apply to RWRelu.

## 4.2 RKHS Weightings as STAR models

Instantiating an RKHS Weighting model consists in choosing its four components $(\mathcal{W}, \phi, \mathcal{K}, p)$. We will show how to instantiate an RKHS Weighting which is also a STAR model, first in full generality, and then by
transforming an already existing instantiation into a STAR model. By having a STAR model which is also an RKHS Weighting, we keep all of the benefits of the RKHS Weightings, i.e. its theoretical guarantees, learning algorithms, pruning methods (Dubé & Marchand, 2024), but also gain interpretability through the computation of the Shapley values, using Algorithm 1.

The defining characteristic of STAR models is that they are built with functions that use only a subset of the variables at a time. The key to defining an RKHS Weighting which is also a STAR model is therefore the base predictor $\phi$. By choosing a base predictor which can flexibly change the features that it considers, we can obtain a STAR model. We now go through each ingredient and see how they must be defined.

1. The parameter space $\mathcal{W}$ is accompanied by the space of feature subsets $\mathcal{P}([n])$ (the family of all subsets of $[n]$).

2. The base predictor takes as input both a parameter $w$ and a feature subset $I$, as well as an instance $x$; $\phi((w, I), x) := \phi_I(w, x_I)$ could be any predictor defined on the feature subset $I$.

3. The kernel must take into account the feature subsets, meaning it is now a function of the form $\mathcal{K} : (\mathcal{W} \times \mathcal{P}([n])) \times (\mathcal{W} \times \mathcal{P}([n])) \to \mathbb{R}$.

4. The distribution $p$ on $\mathcal{W} \times \mathcal{P}([n])$ generates both parameters and feature subsets.

All in all, the instantiation takes the form $(\mathcal{W} \times \mathcal{P}([n]), \phi, \mathcal{K}, p)$. Now, to show that this is indeed a STAR model. Given a weight function $\alpha = \sum_{t=1}^{T} a_t \mathcal{K}((w_t, I_t), \cdot)$, the output of the predictor on an arbitrary instance $x$ is given by:

$$ 
\begin{aligned}
\Lambda\alpha(x) &= \mathbb{E}_{(w,I)\sim p} [\alpha(w, I)\phi((w, I), x)] \\
&= \sum_{t=1}^{T} a_t \mathbb{E}_{(w,I)\sim p} [\mathcal{K}((w_t, I_t), (w, I))\phi_I(w, x_I)] \\
&= \sum_{J \subseteq [n]} \sum_{t: I_t = J} a_t \mathbb{E}_{(w,I)\sim p} [\mathcal{K}((w_t, I_t), (w, I))\phi_I(w, x_I)].
\end{aligned} 
\quad (26) $$

We have thus expressed the predictor as a sum of functions defined on feature subsets. This is therefore a STAR model, with the output on an arbitrary $x$ given by $\Lambda\alpha(x) = \sum_{J \subseteq [n]} f_J(x_J)$, with:

$$ f_J(x_J) := \sum_{t: I_t = J} a_t \mathbb{E}_{(w,I)\sim p} [\mathcal{K}((w_t, I_t), (w, I))\phi_I(w, x_I)]. \quad (27) $$

This framework for expressing an RKHS Weighting as a STAR model is in full generality. Of course, prior knowledge specific to a given problem can be inserted into the choices for a base predictor, kernel or distribution. For instance, we might be interested only in sequences of features when presented with text, time series, or other sequential data, in which case the distribution $p$ can simply put a weight of 0 to any non-sequential feature subset. This family of STAR models is therefore highly flexible.

## 4.3 Transforming an existing instantiation into a STAR model

Perhaps the simplest way to instantiate an RKHS Weighting as a STAR model using the framework described in the previous section is to take an existing instantiation $(\mathcal{W}, \phi, \mathcal{K}_{\mathcal{W}}, p_{\mathcal{W}})$ which satisfies a few simple conditions:

1. $\mathcal{X} = \mathcal{W} = \mathbb{R}^n$.

2. The base predictor $\phi(w, x)$ is a function of $(w, x)$, i.e. $\phi(w, x) = \phi(\langle w, x \rangle)$.

3. The distribution $p_{\mathcal{W}}$ on $\mathcal{W}$ generates each component independently.
<table><caption>Table 4: Cheat sheet for transforming an existing RKHS Weighting instantiation into a Structured Additive Regression model.</caption><thead><tr><th></th><th>Existing instantiation</th><th>STAR version</th></tr></thead><tbody><tr><td>Parameter space</td><td>W</td><td>W × P([n])</td></tr><tr><td>Base predictor</td><td>φ(w, x) = φ(⟨w, x⟩)</td><td>φ((w, I), x) = φ(⟨w<sub>I</sub>, x<sub>I</sub>⟩)</td></tr><tr><td>Kernel</td><td>K<sub>W</sub></td><td>K((w<sub>1</sub>, I<sub>1</sub>), (w<sub>2</sub>, I<sub>2</sub>)) = K<sub>W</sub>(w<sub>1</sub>, w<sub>2</sub>)1[I<sub>1</sub> = I<sub>2</sub>]</td></tr><tr><td>Distribution</td><td>p<sub>W</sub></td><td>p(w, I) = p<sub>W|I</sub>(w)p<sub>[n]</sub>(I)</td></tr></tbody></table>

In that situation, modifying the base predictor to work on feature subsets is trivially easy:

$$ \phi((w, I), x) := \phi(\langle w_I, x_I \rangle). \qquad (28) $$

In other words, simply ignore all the variables not in the feature subset. As for the distribution $p$, the independence allows us to marginalize away the unwanted variables. The distribution $p(w, I)$ is the product $p_{W|I}(w)p_{[n]}(I)$, where $p_{[n]}$ is a chosen distribution over the feature subsets, and $p_{W|I}(w)$ is the original distribution $p_W$ marginalized to the features in $I$. Finally, we elect to use the indicator kernel $K_{[n]}(I_1, I_2) := 1[I_1 = I_2]$ over the feature subsets, and use the kernel product $K((w_1, I_1), (w_2, I_2)) := K_W(w_1, w_2)K_{[n]}(I_1, I_2)$.<fnref n="4" /> We end up with the new instantiation $(W \times P([n]), \phi, K_W \times K_{[n]}, p)$. (Table 4 summarizes all of this.) Assuming a weight function $\alpha = \sum_{t=1}^T a_t K((w_t, I_t), \cdot)$, the model is:

$$ \Lambda\alpha(x) = \sum_{J \subseteq [n]} \sum_{t: I_t = J} a_t p_{[n]}(I_t) \mathbb{E}_{w \sim p_{W|I_t}} [K_W(w_t, w) \phi(\langle w, x_{I_t} \rangle)]. \qquad (29) $$

The only element that still needs to be further specified is the distribution $p_{[n]}$ on the feature subsets. We list below the two most obvious candidates, but note that there is complete flexibility here for the user of the model to specify any distribution that they deem adequate for the dataset at hand.

1. The simplest option is to set $p_{[n]}$ as the uniform distribution on feature subsets of size equal to $k$, where $k$ becomes a hyperparameter of the algorithm. A slight variation of this is to use feature subsets of size up to $k$, allowing a larger variety of variable interactions.

2. Some datasets, like time series or genomic data, have structure, where adjacent features are correlated in some meaningful way. When presented with such sequential data, it is natural to consider sequences of features. We can then set $p_{[n]}$ as the uniform distribution over sequences of length $k$ (or length up to $k$).

Because of the presence of the factor $p_{[n]}(I_t)$ in Equation (29), using the first option (arbitrary subsets of size $k$) with large $k$ is not recommended. Indeed, since $p_{[n]}(I_t) = 1/(n^k) \approx \frac{1}{n^k}$, the model output tends to zero exponentially as $k$ increases. Large coefficients $a_t$ are then required to offset the small density. This leads to a large norm $\|\alpha\|_H$ of the weight function, and therefore poorer guarantees (Dubé & Marchand, 2024). More intuitively, it is important that a significant portion of the feature subsets $I_t$ in Equation (29) are meaningful interactions. Otherwise, the partial predictors built on these uninformative interactions will themselves be poor. Reducing the size of the space of feature subsets by limiting the distribution $p_{[n]}$ to mostly meaningful interactions is crucial to get a well-performing STAR RKHS Weighting model. Sequential feature subsets achieve this. Since there are only $n-k+1$ sequences of features of length $k$ in an $n$-dimensional dataset, $p_{[n]}$ will always be of the order of $\frac{1}{n}$. When considering sequences of length up to $k$, then $p_{[n]}$ will be about $\frac{1}{kn}$, which is undeniably preferable over $\frac{1}{n^k}$.

<footnote n="4">⁴We use the indicator kernel for two related but distinct reasons. The first is the assumption that, in general, the comparison between parameters, i.e. the kernel evaluation $K_W(w_1, w_2)$, only makes sense if $w_1$ and $w_2$ are parameters over the same features. Second, it ensures that $w_1$ and $w_2$ are the same dimensionality, so that the kernel evaluation $K_W(w_1, w_2)$ is well-defined. There is some abuse of notation here, as we use the same symbol $K_W$ to denote a function over potentially different spaces (e.g. $w_1$ could be from $\mathbb{R}^3$, and $w_2$ from $\mathbb{R}^4$), but deem it benign enough, as it does simplify the expressions.</footnote>
We can use this transformation to immediately obtain STAR models based on instantiations RWSign and RWRelu, which can be learned using any algorithm from Dubé & Marchand (2024). (Note that instantiation RWStumps is already a GAM, so no transformation is required). This makes up a new family of STAR models, which are easy to train, as they do not require explicitly specifying feature interactions when using one of the learning algorithms from Dubé & Marchand (2024).

(For completeness, we describe in the appendix how to obtain the theoretical constants from Dubé & Marc-
hand (2024) for a STAR RKHS model. Once again, this means that the theoretical results of Dubé &
Marchand (2024) apply to the STAR RKHS Weightings constructed in this way.)

5 Experiments

We introduced three new tools in this paper: STAR-SHAP, an algorithm for calculating the Shapley values of a STAR model; RWRelu, a new RKHS Weighting instantiation; and a new family of STAR models, obtained by transforming existing RKHS Weighting instantiations into STAR models. We present in the next few sections three experiments to test these new tools. First, we compare the computation time of STAR-SHAP to other Shapley value algorithms in Section 5.2. Then, Section 5.3 presents a simple regression performance experiment, designed to compare various regression algorithms to the new RWRelu and STAR RKHS Weighting models. Next, the experiment of Section 5.4 investigates a scenario (time series regression) in which we expect the STAR RKHS Weightings to be better suited than the simple regression problems of the Section 5.3 experiment. Finally, Section 5.5 compares the explicit Shapley values of a STAR RKHS model to those of a conventional model, both trained on the same dataset.

**5.1 Generating synthetic datasets**

To generate our synthetic datasets, we took inspiration from the sum of unanimity models used by Fumagalli et al. (2023). For given values of $k$ (the size of the feature subsets), $n$ (the total number of variables), and $T$ (the desired number terms), we generate a $k$-STAR synthetic model over boolean instances ($\mathcal{X} = \{\pm1\}^n$) as $f(x) = \sum_{i=1}^{T} a_i 1[\forall j \in I_i, x_j = 1]$, where $I_1, \dots, I_T$ are uniformly and independently generated feature subsets of size $k$, and $a_1, \dots, a_T$ are coefficients uniformly generated in $[0, 1]$. To go along with every random model, we also generate a random dataset of boolean-valued instances $S \in \mathcal{X}^m$, where each component of an instance $x$ is independently and uniformly sampled from $\{\pm1\}$.

**5.2 Computation time of the Shapley values**

Figure 1 shows the computation time of the Shapley values of 5-STAR models (that is, with five-way interactions between variables) using various algorithms. We notice, of course, the exponential growth of the computation time for SHAP, while STAR-SHAP cares little for the dimensionality when the size of the feature subsets is fixed. For completeness, we also show the computation time growth of using STAR-SHAP to calculate the Shapley values of a model which uses the full n-way interactions. We observe the expected exponential behavior, mirroring SHAP. Meanwhile, Figure 1b clearly illustrates the fact that STAR-SHAP lacks an explicit dependency on the dimensionality of the data, making it much more efficient than other algorithms at higher dimensionalities.

However, Figure 1 does not account for the quality of the approximations, and only looks at 5-STAR
models. Unbiased KernelSHAP and SHAP-IQ are not so much slower than STAR-SHAP, but how good
are their approximations of the exact Shapley values? Since Unbiased KernelSHAP and SHAP-IQ coincide
when calculating 1-Shapley values (Fumagalli et al., 2023), possess theoretical guarantees, and have better
empirical performance than KernelSHAP (according to Fumagalli et al. (2023)), we assess the quality of
SHAP-IQ in Figure 2. The figure shows how the advantage of STAR-SHAP increases with dimensionality.
Indeed, not only does the computation time of SHAP-IQ increase, but also the quality of its approximations
degrades, requiring additional computation time to return good approximations. Meanwhile, the values
returned by STAR-SHAP are always exact, and are calculated efficiently.
<figure><img src="image_3.png" /><figcaption>Figure 1: Computation time of the Shapley values of a Structured Additive Regression model with respect to the dimensionality of the dataset. For every value of *n*, a 5-STAR random synthetic model with *T* = 50 terms was generated as described in Section 5.1, along with a dataset of 10 examples. Then, all algorithms were used to calculate the Shapley values of the model for all variables and examples. Similarly, *n*-STAR models, using all *n* possible interactions, were generated, and STAR-SHAP applied to this model. SHAP and the approximation algorithms were also applied to this model for values *n* < 5, for which a 5-STAR model cannot exist. Figure 1b is the rescaled continuation of Figure 1a without the exponentially increasing algorithms. Every point is the average of 20 runs to smooth out any possible variability due to hardware or CPU conflicts. We used the SHAP-IQ implementation (Muschalik et al., 2024) of KernelSHAP, Unbiased KernelSHAP and SHAP-IQ itself, with a budget of 256. Note that the complexity of KernelSHAP appears linear at this scale, but it is in fact cubic in *n* (see Table 2).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: One minus the coefficient of determination ($R^2$) of the Shapley values of synthetic models approximated by SHAP-IQ with regard to the computation time (specified indirectly through the budget hyperparameter). The models were generated as described in Section 5.1, with a fixed number of terms in the models ($T = 100$), and varying dimensionality $n$. For each value of $k$, the size of the feature subsets in the STAR model, the Shapley values of a synthetic model on a dataset of size 10 were calculated approximately by SHAP-IQ ($k$-SHAP-IQ), and exactly by STAR-SHAP ($k$-STAR-SHAP). The $R^2$ was then calculated between the vectorized $10 \times n$ matrices of approximated and exact Shapley values. A vertical dotted line labeled $k$-STAR-SHAP indicates the time STAR-SHAP took to calculate the exact Shapley values for that value of $k$. Every point is the average of 10 independently seeded runs. Figure 2a: on low-dimensionality datasets ($n = 10$ variables, $T = 100$ terms in the model), SHAP-IQ starts generating Shapley values with an $R^2$ coefficient larger than 0 faster than STAR-SHAP at $k = 3$. This means that STAR-SHAP is the better algorithm choice for $k \le 3$, and then gradually loses its merits as $k$ increases and SHAP-IQ is allowed more computation time. In Figure 2b, more total variables were used in the models ($n = 100, T = 100$). Now, STAR-SHAP is unambiguously preferable up to about $k = 5$, after which SHAP-IQ starts becoming competitive if perfectly exact Shapley values are not required. Finally, Figure 2c used $n = 500$ total variables (with $T = 100$ terms) to build the models. Even for $k = 10$, which is an unusually high number of variable interactions for a model, SHAP-IQ generates barely acceptable Shapley values when given the same amount of time as STAR-SHAP. For any smaller value of $k$ and the same computation time as STAR-SHAP, SHAP-IQ yielded extremely poor quality approximations of the Shapley values.</figcaption></figure>
<table><caption>Table 5: Datasets of the regression experiment.</caption><thead><tr><th></th><th>Training size</th><th>Test size</th><th>Dimensionality</th><th>Source (clickable)</th></tr></thead><tbody><tr><td>abalone</td><td>3132</td><td>1045</td><td>10</td><td>UCI</td></tr><tr><td>diabetes</td><td>331</td><td>111</td><td>10</td><td>Scikit-learn</td></tr><tr><td>housing</td><td>15480</td><td>5160</td><td>8</td><td>Scikit-learn</td></tr><tr><td>concrete</td><td>772</td><td>258</td><td>8</td><td>UCI</td></tr><tr><td>conductivity</td><td>15947</td><td>5316</td><td>81</td><td>UCI</td></tr><tr><td>wine</td><td>133</td><td>45</td><td>13</td><td>UCI</td></tr></tbody></table>

Of course, these computation times and accuracies are highly dependent on the models and datasets that were generated, but they strongly hint at a general trend. If the size of the feature subsets is limited, then the computation time of STAR-SHAP will be fast. Since it always returns exact values, it will be the preferable algorithm in that scenario. This is true regardless of model size. However, the more total variables in the model, the more feature interactions can be present before needing to consider using an approximation algorithm, due to their slower convergence rate at higher dimensionalities. This makes STAR-SHAP especially valuable in the presence of high-dimensional datasets.

## 5.3 Regression performance

Three new models arise from Section 4: the RWRelu, and the STAR versions of RWSign and RWRelu. This calls for a simple regression experiment comparing the performance of these models to more bread and butter regression techniques, especially interpretable ones. Using 5-fold cross-validation to select hyperparameters, we trained a variety of models on several regression datasets taken from the UCI Machine Learning Repository (Dua & Graff, 2017) and Scikit-learn (Pedregosa et al., 2011). See Table 5 for information on the datasets, and the appendix for the detailed experimental setup.

We can look to the results in Table 6 to answer a few important questions.

**Q1.** How does RWRelu, the ReLU version of an RKHS Weighting, compare to RWSign and RWStumps, the two original instantiations, for regression tasks?

**R1.** We can see that RWRelu performs equally or better than RWSign on all but one dataset (Wine). On the other hand, RWStumps outperforms RWRelu on 4 of the 6 datasets. The STAR versions of RWRelu and RWSign outperform the other (for equal *k*) on some datasets, and not others, neither systematically standing out. These observations give some weight, but not much, to the conjecture that RWRelu is better suited to regression than the other two RKHS Weighting instantiations. Rather, they might all be complementary, each with better and worse use cases.

**Q2.** How do the STAR versions of RKHS Weighting models perform in comparison to their original counterparts?

**R2.** There is a moderate loss of performance going from RWRelu to its 1-STAR version, a cost that might be acceptable to pay for the immense gain in interpretability. However, increasing *k* does improve the performance significantly on all but one dataset, Conductivity, and 3-STAR RWRelu is the best performing STAR version of the model, which shows that adding interactions can improve the quality of the model. The extremely poor performance of 3, 4 and 5-STAR models on Conductivity is caused by the higher dimensionality of the dataset. (See this footnote<fnref n="5" /> for the in-depth explanation.) Indeed, larger feature subsets should not blindly be used except on low-dimensional datasets, or failure of learning might occur. As for the STAR versions of RWSign, they are less convincing than RWRelu, with the 1-STAR version in particular exhibiting extremely poor performance overall.

<footnote n="5">⁵We recall the presence of a factor $p_{[n]}(I_t)$ in Equation (29), the expression for a STAR RKHS Weighting. When using the uniform distribution over all possible feature subsets of size *k*, we have $p_{[n]}(I_t) = 1/(n^k) \approx \frac{1}{n^k}$, which is an extremely small value for high-dimensionality datasets, leading to vanishing outputs of the model. And indeed, the offending dataset, Conductivity, is far higher-dimensional than the other ones, proving the point.</footnote>
<table><caption>Table 6: Comparison of the regression performance of STAR RKHS Weightings and the ReLU instantiation, learned using Algorithm 3 of Dubé & Marchand (2024) (the Least Squares fit) with T = 5000, to a variety of other regression algorithms. k-STAR signifies the STAR version of the model using feature subsets of size k. The table values are the test R² (coefficient of determination). The highest (best) values for each dataset have been bolded.</caption><thead><tr><th>Dataset<br>Algorithm</th><th>abalone</th><th>diabetes</th><th>housing</th><th>concrete</th><th>conductivity</th><th>wine</th></tr></thead><tbody><tr><td>DecisionTreeRegressor</td><td>0.499</td><td>0.105</td><td>0.677</td><td>0.796</td><td>0.882</td><td>0.867</td></tr><tr><td>EBM</td><td>0.552</td><td>0.331</td><td><strong>0.822</strong></td><td><strong>0.931</strong></td><td><strong>0.903</strong></td><td>0.897</td></tr><tr><td>KernelRidge</td><td><strong>0.587</strong></td><td>0.282</td><td>0.771</td><td>0.831</td><td>0.872</td><td><strong>0.951</strong></td></tr><tr><td>LinearRegression</td><td>0.543</td><td>0.359</td><td>0.591</td><td>0.623</td><td>0.736</td><td>0.804</td></tr><tr><td>RWRelu</td><td>0.577</td><td>0.349</td><td>0.734</td><td>0.835</td><td>0.842</td><td>0.832</td></tr><tr><td>RWRelu 1-STAR</td><td>0.555</td><td>0.314</td><td>0.662</td><td>0.780</td><td>0.777</td><td>0.808</td></tr><tr><td>RWRelu 2-STAR</td><td>0.577</td><td>0.305</td><td>0.731</td><td>0.847</td><td>0.805</td><td>0.843</td></tr><tr><td>RWRelu 3-STAR</td><td>0.580</td><td>0.332</td><td>0.744</td><td>0.847</td><td>0.568</td><td>0.845</td></tr><tr><td>RWRelu 4-STAR</td><td>0.579</td><td>0.350</td><td>0.737</td><td>0.846</td><td>0.032</td><td>0.835</td></tr><tr><td>RWRelu 5-STAR</td><td>0.578</td><td>0.334</td><td>0.729</td><td>0.844</td><td>0.001</td><td>0.837</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.577</td><td>0.298</td><td>0.732</td><td>0.808</td><td>0.779</td><td>0.826</td></tr><tr><td>RWSign</td><td>0.578</td><td>0.348</td><td>0.711</td><td>0.828</td><td>0.728</td><td>0.893</td></tr><tr><td>RWSign 1-STAR</td><td>0.350</td><td>0.302</td><td>0.354</td><td>0.471</td><td>0.685</td><td>0.750</td></tr><tr><td>RWSign 2-STAR</td><td>0.579</td><td>0.369</td><td>0.696</td><td>0.835</td><td>0.816</td><td>0.868</td></tr><tr><td>RWSign 3-STAR</td><td><strong>0.586</strong></td><td>0.339</td><td>0.720</td><td>0.840</td><td>0.468</td><td>0.912</td></tr><tr><td>RWSign 4-STAR</td><td>0.585</td><td>0.349</td><td>0.708</td><td>0.794</td><td>0.003</td><td>0.896</td></tr><tr><td>RWSign 5-STAR</td><td>0.580</td><td>0.333</td><td>0.705</td><td>0.770</td><td>-0.000</td><td>0.899</td></tr><tr><td>RWSign [1, 2, 3, 4, 5]-STAR</td><td>0.580</td><td><strong>0.399</strong></td><td>0.698</td><td>0.826</td><td>0.694</td><td>0.840</td></tr><tr><td>RWStumps</td><td>0.567</td><td>0.339</td><td>0.768</td><td>0.911</td><td>0.859</td><td>0.864</td></tr><tr><td>SVR</td><td>0.564</td><td>0.293</td><td>0.770</td><td>0.799</td><td>0.862</td><td><strong>0.950</strong></td></tr></tbody></table>
<table><caption>Table 7: Datasets of the time series prediction experiment.</caption><thead><tr><th></th><th>Training size</th><th>Test size</th><th>Dimensionality</th></tr></thead><tbody><tr><td>ChlorineConcentration</td><td>6990</td><td>57585</td><td>10</td></tr><tr><td>Computers</td><td>16434</td><td>16434</td><td>10</td></tr><tr><td>ECG5000</td><td>6487</td><td>58487</td><td>10</td></tr><tr><td>FacesUCR</td><td>2388</td><td>24588</td><td>10</td></tr><tr><td>LargeKitchenAppliances</td><td>24684</td><td>24684</td><td>10</td></tr><tr><td>MelbournePedestrian</td><td>3579</td><td>7314</td><td>10</td></tr></tbody></table>

**Q3.** How does the new RWRelu instantiation perform compared to other regression algorithms?

**R3.** The two longstanding, bread and butter regression algorithms that are Kernel Ridge Regression and Support Vector Regression outperform RWRelu respectively on 4 of the 6, and 3 of the 6 datasets. This is far from an indictment of RWRelu.

**Q4.** Can interpretable RKHS Weightings perform similarly or better than other interpretable methods?

**R4.** The Explainable Boosting Machine (EBM) of (Lou et al., 2013) (software package described in Nori et al. (2019)) appears overall to be the best at both accuracy and interpretability, at least in this experiment. Kernel Ridge Regression is best on two datasets, and SVR, 3-STAR RWSign, and 1-to-5 STAR RWSign are each best on one dataset. Therefore, we can say that interpretable RKHS Weightings are at least sometimes competitive.

To conclude this experiment, we see that transforming an RKHS Weighting into a STAR model is a valid way of learning an interpretable model, at a relatively small cost of performance. However, other methods, especially Explainable Boosting Machines, often outperform these RKHS Weightings in both accuracy and interpretability.

## 5.4 Time series prediction performance

Perhaps the biggest takeaway from Table 6 is the excellence of the Explainable Boosting Machine on both the accuracy and interpretability fronts. The RKHS Weightings did not shine particularly brightly. However, the most important characteristic of RKHS Weightings is their high flexibility in instantiating the model. We can improve performance by utilizing an instantiation better adapted to the problem at hand. In the specific case of STAR RKHS Weightings, a tool we can play with is the distribution on the feature subsets. For instance, instead of using the uniform distribution over all subsets of a certain length, we can use a uniform distribution over *sequences* of adjacent features (e.g. features (2, 3, 4), or (n - 2, n - 1, n)). This makes the STAR RKHS instantiation much better suited to solve regression problems for which adjacent data features are correlated, such as time series. Therefore, a better comparison between STAR RKHS Weightings and the Explainable Boosting Machine is on such datasets, using a distribution over sequences of features. This is what we do here.

Specifically, we run both the Explainable Boosting Machine (EBM) and STAR RKHS Weightings on time series datasets taken from the UCR Archive (Dau et al., 2018). Each dataset was cut into segments of length 11, so that each instance *x* has length 10, and the value to predict *y* is the 11-th value. See Table 7. Hyperparameters were selected by 5-fold cross-validation. See the appendix for the detailed experimental setup.

The results of this experiment can be found in Table 8. The RKHS Weightings fare much better in this experiment than the previous. Here, they (slightly) outperform the EBM on 5 of the 6 datasets. On the other hand, the EBM does systematically reach the highest training $R^2$, suggesting that some overfitting has taken place, even though the parameters were chosen by cross-validation. Perhaps a different set of candidate parameters would have led to a smaller generalization gap. Finally, we notice the expected trend that larger feature subsets on the STAR RKHS Weightings lead to better results, with the full RWRelu instantiation
<table><caption>Table 8: Comparison of the time series prediction performance of STAR RKHS Weightings and the ReLU instantiation, learned using Algorithm 3 of Dubé & Marchand (2024) (the Least Squares fit) with T = 5000, to the Explainable Boosting Machine (EBM). k-STAR signifies the STAR version of the model using feature subsets of size k. The highest (best) values for each dataset have been bolded.</caption><thead><tr><th>Dataset</th><th>Algorithm</th><th>Training R<sup>2</sup></th><th>Test R<sup>2</sup></th><th>Training time (s)</th></tr></thead><tbody><tr><td rowspan="8">ChlorineConcentration</td><td>EBM</td><td><strong>0.955</strong></td><td>0.916</td><td>63.853</td></tr><tr><td>RWRelu</td><td>0.934</td><td><strong>0.928</strong></td><td>3.307</td></tr><tr><td>RWRelu 1-STAR</td><td>0.638</td><td>0.642</td><td>36.375</td></tr><tr><td>RWRelu 2-STAR</td><td>0.813</td><td>0.805</td><td>40.390</td></tr><tr><td>RWRelu 3-STAR</td><td>0.907</td><td>0.899</td><td>44.500</td></tr><tr><td>RWRelu 4-STAR</td><td>0.928</td><td>0.920</td><td>48.720</td></tr><tr><td>RWRelu 5-STAR</td><td>0.939</td><td><strong>0.927</strong></td><td>54.159</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.924</td><td>0.906</td><td>19.823</td></tr><tr><td rowspan="8">Computers</td><td>EBM</td><td><strong>0.882</strong></td><td>0.842</td><td>16.291</td></tr><tr><td>RWRelu</td><td>0.821</td><td><strong>0.851</strong></td><td>6.085</td></tr><tr><td>RWRelu 1-STAR</td><td>0.786</td><td>0.838</td><td>44.742</td></tr><tr><td>RWRelu 2-STAR</td><td>0.802</td><td>0.840</td><td>50.501</td></tr><tr><td>RWRelu 3-STAR</td><td>0.811</td><td>0.848</td><td>55.049</td></tr><tr><td>RWRelu 4-STAR</td><td>0.814</td><td>0.847</td><td>58.744</td></tr><tr><td>RWRelu 5-STAR</td><td>0.818</td><td>0.848</td><td>64.858</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.817</td><td><strong>0.849</strong></td><td>29.086</td></tr><tr><td rowspan="8">ECG5000</td><td>EBM</td><td><strong>0.943</strong></td><td><strong>0.886</strong></td><td>23.525</td></tr><tr><td>RWRelu</td><td>0.900</td><td>0.882</td><td>3.144</td></tr><tr><td>RWRelu 1-STAR</td><td>0.843</td><td>0.841</td><td>36.361</td></tr><tr><td>RWRelu 2-STAR</td><td>0.877</td><td>0.871</td><td>40.024</td></tr><tr><td>RWRelu 3-STAR</td><td>0.883</td><td>0.875</td><td>43.769</td></tr><tr><td>RWRelu 4-STAR</td><td>0.893</td><td>0.879</td><td>48.122</td></tr><tr><td>RWRelu 5-STAR</td><td>0.899</td><td>0.882</td><td>53.384</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.894</td><td>0.882</td><td>19.157</td></tr><tr><td rowspan="8">FacesUCR</td><td>EBM</td><td><strong>0.868</strong></td><td>0.661</td><td>7.674</td></tr><tr><td>RWRelu</td><td>0.715</td><td><strong>0.704</strong></td><td>1.929</td></tr><tr><td>RWRelu 1-STAR</td><td>0.679</td><td>0.698</td><td>33.011</td></tr><tr><td>RWRelu 2-STAR</td><td>0.690</td><td>0.702</td><td>35.613</td></tr><tr><td>RWRelu 3-STAR</td><td>0.699</td><td><strong>0.705</strong></td><td>38.827</td></tr><tr><td>RWRelu 4-STAR</td><td>0.700</td><td><strong>0.705</strong></td><td>42.990</td></tr><tr><td>RWRelu 5-STAR</td><td>0.709</td><td><strong>0.706</strong></td><td>49.571</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.704</td><td><strong>0.706</strong></td><td>15.570</td></tr><tr><td rowspan="8">LargeKitchenAppliances</td><td>EBM</td><td><strong>0.840</strong></td><td>0.743</td><td>19.103</td></tr><tr><td>RWRelu</td><td>0.813</td><td>0.743</td><td>8.775</td></tr><tr><td>RWRelu 1-STAR</td><td>0.764</td><td>0.724</td><td>52.104</td></tr><tr><td>RWRelu 2-STAR</td><td>0.786</td><td><strong>0.748</strong></td><td>56.938</td></tr><tr><td>RWRelu 3-STAR</td><td>0.777</td><td>0.736</td><td>60.735</td></tr><tr><td>RWRelu 4-STAR</td><td>0.803</td><td>0.746</td><td>68.894</td></tr><tr><td>RWRelu 5-STAR</td><td>0.810</td><td>0.745</td><td>75.208</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.812</td><td><strong>0.749</strong></td><td>37.741</td></tr><tr><td rowspan="8">MelbournePedestrian</td><td>EBM</td><td><strong>0.979</strong></td><td><strong>0.924</strong></td><td>9.731</td></tr><tr><td>RWRelu</td><td>0.946</td><td><strong>0.925</strong></td><td>2.335</td></tr><tr><td>RWRelu 1-STAR</td><td>0.888</td><td>0.888</td><td>35.930</td></tr><tr><td>RWRelu 2-STAR</td><td>0.921</td><td>0.902</td><td>36.519</td></tr><tr><td>RWRelu 3-STAR</td><td>0.943</td><td><strong>0.923</strong></td><td>40.523</td></tr><tr><td>RWRelu 4-STAR</td><td>0.943</td><td><strong>0.925</strong></td><td>44.806</td></tr><tr><td>RWRelu 5-STAR</td><td>0.943</td><td><strong>0.924</strong></td><td>50.667</td></tr><tr><td>RWRelu [1, 2, 3, 4, 5]-STAR</td><td>0.941</td><td>0.919</td><td>16.232</td></tr></tbody></table>
often being the best (although by a very small margin). Indeed, more complex interactions allow for greater expressivity of the model, so that it can better fit the data.

It is important to note that RWRelu is just one RKHS Weighting instantiation of a theoretically limitless number. It is certain that other, yet to be defined, instantiations would perform much better than RWRelu in this experiment. Indeed, it is still an open question how to craft instantiations that can perform particularly well on a given problem type (such as time series regression), or even for a specific dataset. More research is certainly required.

## 5.5 Shapley values comparison

Our final experiment shows how a STAR RKHS Weighting model could be used in practice, in comparison to a more traditional model, and how its predictions could be interpreted using Algorithm 1 (STAR-SHAP). Figure 3 compares the Shapley values of the Explainable Boosting Machine (EBM) model and the 3-STAR RWRelu model that were trained on the California Housing dataset in Section 5.3. We see that both models agree almost exactly on the order of importance of the variables, and that the bee swarms for each variable are similarly shaped between both models. This is despite the models being fundamentally different predictors: the EBM is a sum of Trees built on single variables, while the 3-STAR RKHS Weighting is a sum of partial predictors of the form given by Equation 25, using three variables each.

<figure><img src="image_5.png" /><figcaption>Figure 3: Beeswarm plots of the Shapley values on California Housing of the Explainable Boosting Machine (Figure 3a) and 3-STAR RWRelu (Figure 3b). Models are taken from the experiment of Section 5.3. The Shapley values of 1000 test examples were calculated using SHAP (Figure 3a) and STAR-SHAP (Figure 3b), using a background of 1000 training examples. The bee swarms, top to bottom, are in decreasing order of variable importance. This order is almost the same for both models, with only HouseAge and AveRooms being interchanged.</figcaption></figure>

Figure 4 compares the Shapley values of both models on a single instance. While the predicted values are noticeably different, both models agree on the variables to consider: Longitude has a strong negative contribution to the prediction, Latitude an almost equally strong positive contribution, and AveOccup a small negative contribution. Other variables are mostly negligible. One could therefore reach the same interpretative conclusions from seeing the RKHS Weighting plot as the EBM one. This gives further credence to RKHS Weightings as a promising model family, especially an interpretable one when cast as STAR models.
<figure><img src="image_6.png" /><figcaption>Figure 4: Waterfall plots of the Shapley values of a single instance $x$ from Figure 3. Figure 4a: Explainable Boosting Machine. Figure 4b: 3-STAR RWRelu. The true $f(x)$ value for that instance is $-0.6106$.</figcaption></figure>

## 5.6 Discussion

A number of conclusions can be reached from our experiments.

Figure 1 demonstrates that STAR-SHAP is the preferred algorithm for calculating the Shapley values of a model with moderate variable interactions, especially so the more total variables are present in the data. Figure 2 more clearly shows how the advantages of STAR-SHAP increase with the total number of variables in the model, as the quality of the approximations returned by approximation algorithms degrades with dimensionality.

Table 6 shows that RKHS Weightings and STAR RKHS Weightings can reach moderately good performance on arbitrary datasets. They have a better showing in Table 8, where prior knowledge about time series led to improved performance, though it remains questionable whether RKHS Weightings should be used over conventional models.

Figures 3 and 4 show that STAR-SHAP and STAR RKHS Weightings can yield interpretable models.

All in all, while these experiments do not show that RKHS Weightings should be used over other, more conventional, tried and tested models, they are a step forward in understanding how to use RKHS Weightings, and a proof that they be high accuracy, theoretically grounded (as per Dubé & Marchand (2024)), and even interpretable models. Further research to expand the usability of RKHS Weightings, as well as find their unique strengths and applications, is warranted.

## 6 Limitations and future work

Our work presents a number of limitations, some of which can be addressed with future research.

**Algorithm 1 complexity.** In theory, Algorithm 1 can calculate the Shapley values of any STAR model, but the computational cost is exponential in the size of the feature subsets. This limits the size of the feature subsets that can be used in the model, but Figures 1 and 2 show that, for moderate interaction sizes, it is the preferred algorithm to use to calculate the Shapley values of STAR models, especially so the more total variables are in the data.

**Instantiating RKHS Weightings.** The main weakness of RKHS Weightings remains the requirement of solving a complex integral for any new instantiation. We have introduced one new instantiation, but it
is still an open problem to figure out how to access the potential of this family of models. Indeed, it is almost certain that better results than those found in Tables 6 and 8 can be obtained by using different instantiations. As Dubé & Marchand (2024) suggested, Monte Carlo approximation of the expectations might provide a good enough solution, and warrants an in-depth look.

**Prior knowledge.** As seen in Table 6, the most generic version of STAR RKHS Weightings, i.e. using the uniform distribution on all possible feature subsets of size *k*, have somewhat limited use, since their output scales as 1/($\binom{n}{k}$) ($n$ being the dimensionality of the dataset). As we have already addressed, the best option is to insert prior knowledge into the model, as we have done with time series datasets by considering sequences of variables. Extending this principle to other types of structured data would increase the applicability of STAR RKHS Weighting models.

# 7 Conclusion

In this paper, we derived an efficient algorithm, named STAR-SHAP, for calculating the Shapley values of any Structured Additive Regression (STAR) model. We demonstrated its efficiency compared to other Shapley value algorithms when the number of variable interactions is limited. This advantage is further increased on high-dimensional datasets. We introduced a new RKHS Weightings instantiation, and showed how to obtain RKHS Weightings that are STAR models, giving rise to a new family of STAR models, all of which are now interpretable thanks to STAR-SHAP. We tested the prediction performance of the introduced models. While the STAR RKHS Weightings did not rise to the level of state of the art interpretable algorithms on generic regression datasets, they proved quite capable in the context of time series prediction when infused with adequate prior knowledge. Further work to increase the breadth of usable RKHS Weightings instantiations is warranted, as well as work to expand the types of datasets where RKHS Weightings can perform well.

# Acknowledgments

This work is supported by the DEEL Project CRDPJ 537462-18 funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Consortium for Research and Innovation in Aerospace in Québec (CRIAQ), together with its industrial partners Thales Canada inc, Bell Textron Canada Limited, CAE inc and Bombardier inc.<fnref n="6" />

Special thanks to the anonymous TMLR reviewers whose comments and suggestions led to significant improvements of the manuscript.

<footnote n="6">⁶https://deel.quebec</footnote>
## References

Alain Berlinet and Christine Thomas-Agnan. *Reproducing kernel Hilbert spaces in probability and statistics*. Springer Science & Business Media, 2011.

Sebastian Bordt and Ulrike von Luxburg. From shapley values to generalized additive models and back. In *International Conference on Artificial Intelligence and Statistics*, pp. 709–745. PMLR, 2023.

Eugenio Brusa, Luca Cibrario, Cristiana Delprete, and Luigi Gianpio Di Maggio. Explainable ai for machine fault diagnosis: understanding features' contribution in machine learning models for industrial condition monitoring. *Applied Sciences*, 13(4):2038, 2023.

Rich Caruana, Yin Lou, Johannes Gehrke, Paul Koch, Marc Sturm, and Noemie Elhadad. Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission. In *Proceedings of the 21st ACM SIGKDD international conference on knowledge discovery and data mining*, pp. 1721–1730, 2015.

Ian Covert and Su-In Lee. Improving kernelshap: Practical shapley value estimation via linear regression. *arXiv preprint arXiv:2012.01536*, 2020.

Hoang Anh Dau, Eamonn Keogh, Kaveh Kamgar, Chin-Chia Michael Yeh, Yan Zhu, Shaghayegh Gharghabi, Chotirat Ann Ratanamahatana, Yanping, Bing Hu, Nurjahan Begum, Anthony Bagnall, Abdullah Mueen, Gustavo Batista, and Hexagon-ML. The ucr time series classification archive, October 2018. https://www.cs.ucr.edu/~eamonn/time_series_data_2018/.

Dheeru Dua and Casey Graff. Uci machine learning repository, 2017. URL http://archive.ics.uci.edu/ml.

Gabriel Dubé and Mario Marchand. RKHS Weightings of Functions. working paper or preprint, July 2024. URL https://hal.science/hal-04236058.

Jerome H Friedman. Greedy function approximation: a gradient boosting machine. *Annals of statistics*, pp. 1189–1232, 2001.

Daniel Vidali Fryer, Inga Strümke, and Hien Nguyen. Shapley values for feature selection: The good, the bad, and the axioms. *IEEE Access*, PP:1–1, 2021. URL https://api.semanticscholar.org/CorpusID:231985458.

Fabian Fumagalli, Maximilian Muschalik, Patrick Kolpaczki, Eyke Hüllermeier, and Barbara Hammer. Shap- iq: Unified approximation of any-order shapley interactions, 2023. URL https://arxiv.org/abs/2303.01179.

Trevor J Hastie. Generalized additive models. In *Statistical models in S*, pp. 249–307. Routledge, 2017.

Giles Hooker. Discovering additive structure in black box functions. In *Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining*, pp. 575–580, 2004.

Ling Huang, Jinzhu Jia, Bin Yu, Byung-Gon Chun, Petros Maniatis, and Mayur Naik. Predicting execution time of computer programs using sparse polynomial regression. *Advances in neural information processing systems*, 23, 2010.

Xuanxiang Huang and Joao Marques-Silva. The inadequacy of shapley values for explainability, 2023.

Amir Khorrami Chokami and Giovanni Rabitti. An exact game-theoretic variable importance index for generalized additive models. *Journal of Computational and Graphical Statistics*, pp. 1–10, 2024.

Sotiris B Kotsiantis. Decision trees: a recent overview. *Artificial Intelligence Review*, 39:261–283, 2013.

F Kuo, I Sloan, Grzegorz Wasilkowski, and Henryk Woźniakowski. On decompositions of multivariate functions. *Mathematics of computation*, 79(270):953–966, 2010.
Gabriel Laberge, Yann Pequignot, Alexandre Mathieu, Foutse Khomh, and Mario Marchand. Partial order in chaos: Consensus on feature attributions in the rashomon set. *Journal of Machine Learning Research*, 24(364):1–50, 2023.

Gabriel Laberge, Yann Pequignot, Mario Marchand, and Foutse Khomh. Tackling the XAI Disagreement Problem with Regional Explanations. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, volume 238, Valencia, Spain, May 2024. URL https://hal.science/hal-04480870.

Jimmy Lin, Chudi Zhong, Diane Hu, Cynthia Rudin, and Margo Seltzer. Generalized and scalable optimal sparse decision trees. In *International Conference on Machine Learning*, pp. 6150–6160. PMLR, 2020.

Yin Lou, Rich Caruana, Johannes Gehrke, and Giles Hooker. Accurate intelligible models with pairwise interactions. In *Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining*, pp. 623–631, 2013.

Scott M Lundberg and Su-In Lee. A unified approach to interpreting model predictions. *Advances in neural information processing systems*, 30, 2017.

Scott M Lundberg, Gabriel G Erion, and Su-In Lee. Consistent individualized feature attribution for tree ensembles. *arXiv preprint arXiv:1802.03888*, 2018.

Michael Mayer, Steven C Bourassa, Martin Hoesli, and Donato Scognamiglio. Structured additive regression and tree boosting. *Swiss Finance Institute Research Paper*, (21-83), 2021.

Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller. Methods for interpreting and understanding deep neural networks. *Digital signal processing*, 73:1–15, 2018.

Maximilian Muschalik, Hubert Baniecki, Fabian Fumagalli, Patrick Kolpaczki, Barbara Hammer, and Eyke Hüllermeier. shapiq: Shapley interactions for machine learning. In *The Thirty-eight Conference on Neural Information Processing Systems Datasets and Benchmarks Track*, 2024. URL https://openreview.net/forum?id=knxGmi6SJi.

Harsha Nori, Samuel Jenkins, Paul Koch, and Rich Caruana. Interpretml: A unified framework for machine learning interpretability, 2019. URL https://arxiv.org/abs/1909.09223.

F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. *Journal of Machine Learning Research*, 12:2825–2830, 2011.

Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. " why should i trust you?" explaining the predictions of any classifier. In *Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining*, pp. 1135–1144, 2016.

Cynthia Rudin, Chudi Zhong, Lesia Semenova, Margo Seltzer, Ronald Parr, Jiachang Liu, Srikar Katta, Jon Donnelly, Harry Chen, and Zachery Boner. Amazing things come from having many good models. 2024.

Lloyd S Shapley et al. A value for n-person games. 1953.

Mukund Sundararajan and Amir Najmi. The many shapley values for model explanation. In *International conference on machine learning*, pp. 9269–9278. PMLR, 2020.

Eyal Winter. The shapley value. *Handbook of game theory with economic applications*, 3:2025–2054, 2002.

Simon N Wood, Yannig Goude, and Simon Shaw. Generalized additive models for large data sets. *Journal of the Royal Statistical Society Series C: Applied Statistics*, 64(1):139–155, 2015.

Liang Zhang and Zhelun Chen. Large language model-based interpretable machine learning control in building energy systems. *arXiv preprint arXiv:2402.09584*, 2024.
A Proofs

First, we quickly prove Equation (15).

Proof of Equation (15). We have:

$$
\begin{align*}
\phi_i^{\text{SHAP}}(h, x) &:= \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} [h(r_{\pi; i} \cup \{i\}(z, x)) - h(r_{\pi; i}(z, x))] \\
&= \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} \left[ \left( \sum_{j \in \pi; i \cup \{i\}} f_j(x_j) + \sum_{j \notin \pi; i \cup \{i\}} f_j(z_j) \right) - \left( \sum_{j \in \pi; i} f_j(x_j) - \sum_{j \notin \pi; i} f_j(z_j) \right) \right] \\
&= \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} \left[ \left( \sum_{j \in \pi; i \cup \{i\}} f_j(x_j) - \sum_{j \in \pi; i} f_j(x_j) \right) + \left( \sum_{j \notin \pi; i \cup \{i\}} f_j(z_j) - \sum_{j \notin \pi; i} f_j(z_j) \right) \right] \\
&= \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} \left[ (f_i(x_i)) - \left( \sum_{j \in \pi; i \cup \{i\}} f_j(z_j) - \sum_{j \in \pi; i \cup \{i\}} f_j(z_j) + f_i(z_i) \right) \right] \\
&= f_i(x_i) - \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} [f_i(z_i)] \\
&= f_i(x_i) - \mathbb{E}_{z \sim U(S)} [f_i(z_i)].
\end{align*}
$$

□

The main result to prove is Theorem 3.1, which follows from two lemmas further below.

Proof of Theorem 3.1. We are calculating the Shapley values of a STAR model $h(x) = \sum_{I \subseteq [n]} f_I(x_I)$, which can be written as:

$$
\begin{align*}
\phi_i^{\text{SHAP}}(h, x) &:= \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} [h(r_{\pi; i} \cup \{i\}(z, x)) - h(r_{\pi; i}(z, x))] \\
&= \mathbb{E}_{z \sim U(S)} \mathbb{E}_{\pi \sim U(\Omega)} \left[ \sum_{I \subseteq [n]} f_I(r_{\pi; i} \cup \{i\}(z, x)_I) - \sum_{I \subseteq [n]} f_I(r_{\pi; i}(z, x)_I) \right] \\
&= \mathbb{E}_{z \sim U(S)} \left[ \sum_{I \subseteq [n]} \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i} \cup \{i\}(z, x)_I) - f_I(r_{\pi; i}(z, x)_I)] \right].
\end{align*}
$$

We can immediately cut out most of the computations by noticing two facts:

1. Most functions $f_I$ are 0.

2. If $i \notin I$, then $r_{\pi; i} \cup \{i\}(z, x)_I = r_{\pi; i}(z, x)_I$. Hence, the entire expectation:

$$
E_{\pi \sim U(\Omega)} [f_I(r_{\pi; i} \cup \{i\}(z, x)_I) - f_I(r_{\pi; i}(z, x)_I)] \quad (30)
$$

is 0.

The Shapley values therefore boil down to the following formula:

$$
\phi_i^{\text{SHAP}}(h, x) = \mathbb{E}_{z \sim U(S)} \left[ \sum_{I \subseteq [n]: f_i \neq 0, i \in I} \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i} \cup \{i\}(z, x)_I)] - \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i}(z, x)_I)] \right]. \quad (31)
$$

Lemma A.1 gives the value of $\mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i} \cup \{i\}(z, x)_I)]$, and Lemma A.2 gives the value of $\mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i}(z, x)_I)]$. $\square$
**Lemma A.1.** Consider the definitions of Section 2, assuming that $i \in I$, and the following:

$$ \mathcal{A}^+(i, I) := \{A \subset I | i \in A, 1 \le |A| < |I|\}. \qquad (32) $$

In other words, $\mathcal{A}^+(i, I)$ is the family of strict subsets of $I$ which contain $i$. Then we have:

$$ \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i} \cup \{i\}(z, x)_I)] = \frac{f_I(x_I)}{|I|} + \sum_{A \in \mathcal{A}^+(i, I)} \frac{f_I(r_A \setminus \{i\}(z, x)_I)}{|A| \binom{|I|}{|A|}}. \qquad (33) $$

**Lemma A.2.** Consider the definitions of Section 2, assuming that $i \in I$, as well as:

$$ \mathcal{A}^-(i, I) := \{A \subseteq I | i \in A, 1 < |A| \le |I|\}. \qquad (34) $$

In other words, $\mathcal{A}^-(i, I)$ is the family of subsets of $I$ which contain $i$ and at least one other element. Then we have:

$$ \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i}(z, x)_I)] = \frac{f_I(z_I)}{|I|} + \sum_{A \in \mathcal{A}^-(i, I)} \frac{f_I(r_{A \setminus \{i\}}(z, x)_I)}{|A| \binom{|I|}{|A|}}. \qquad (35) $$

*Proof of Lemma A.1.* We can separate the expectation into two parts:

$$ \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi; i} \cup \{i\}(z, x)_I)] = \underbrace{\frac{1}{n!} \sum_{\substack{\pi \in \Omega \\ I \subseteq \pi; i \cup \{i\}}} f_I(r_{\pi; i} \cup \{i\}(z, x)_I)}_{\text{Term 1}} + \underbrace{\frac{1}{n!} \sum_{\substack{\pi \in \Omega \\ I \setminus (\pi; i \cup \{i\}) \neq \emptyset}} f_I(r_{\pi; i} \cup \{i\}(z, x)_I)}_{\text{Term 2}}. \qquad (36) $$

**Calculating Term 1.** The key to calculating the first half of the right-hand side of the previous equation, which we refer to as Term 1, is to notice that $r_{\pi; i} \cup \{i\}(z, x)_I = x_I$ for all the permutations $\pi$ that satisfy $I \subseteq \pi; i \cup \{i\}$, and so the value $f_I(r_{\pi; i} \cup \{i\}(z, x)_I)$ is constant over all those permutations. It is simply $f_I(r_{\pi; i} \cup \{i\}(z, x)_I) = f_I(x_I)$. All that we need to do is calculate the number of such permutations. There is an elegant combinatorics argument that gives us the solution immediately.

Consider any permutation $\pi$ of $[n]$. This permutation satisfies the condition $I \subseteq \pi; i \cup \{i\}$ if and only if $i$ is to the right of every other variable of $I$. This rightmost position is one possibility of out of $|I|$. This means that exactly 1 out of every $|I|$ permutations satisfies the condition. The number of permutations we seek is $\frac{n!}{|I|}$. We therefore have:

$$ \text{Term 1} = \frac{1}{n!} \frac{n!}{|I|} f_I(x_I) = \frac{f_I(x_I)}{|I|}. \qquad (37) $$

**Calculating Term 2.** Here, we consider all the permutations $\pi$ such that $I \setminus (\pi; i \cup \{i\}) \neq \emptyset$, i.e. the feature subset $I$ contains variables not in $\pi; i \cup \{i\}$. These variables will be replaced when calculating the function. Defining $A := I \cap (\pi; i \cup \{i\})$, the set of shared variables (note that we always have $i \in A$), we have:

$$ r_{\pi; i \cup \{i\}}(z, x)_I = r_A(z, x)_I. $$

The function $f_I(r_{\pi; i \cup \{i\}}(z, x)_I)$ therefore takes a different value for each possible $A$, of which there are:

$$ \begin{aligned} |\mathcal{A}^+(i, I)| &= \sum_{k=0}^{|I|-2} \binom{|I|-1}{k} \\ &= \sum_{k=0}^{|I|-1} \binom{|I|-1}{k} - \binom{|I|-1}{|I|-1} \\ &= 2^{|I|-1} - 1. \end{aligned} \qquad (38) $$
(We can choose 0 up to $|I| - 2$ variables to accompany $i$ in making up $A \in \mathcal{A}^+(i, I)$.) For a given set $A \in \mathcal{A}^+(i, I)$ of shared variables, we must calculate the number of permutations $\pi$ that are such that $A = I \cap (\pi_i \cup \{i\})$. Call that number $N(n, |I|, |A|)$. We have:

$$ \text{Term 2} = \sum_{A \in \mathcal{A}^+(i, I)} \frac{N(n, |I|, |A|)}{n!} f_I(r_A(z, x)_I). \qquad (39) $$

Let us now calculate $N(n, |I|, |A|)$.

1. First, choose the $|I|$ positions for the variables in $I$. There are $\binom{n}{|I|}$ possibilities.

2. Of these $|I|$ positions, there is only one choice for $i$ itself, since the variables in $A \setminus \{i\}$ must be to the left of it, and the variables in $I \setminus A$ to the right.

3. We can permute the $|A| - 1$ variables of $A \setminus \{i\}$, the $|I| - |A|$ variables of $I \setminus A$, and the $n - |I|$ variables in $[n] \setminus I$, for a total of $(|A| - 1)!(|I| - |A|!(n - |I|)!)$ permutations.

Assembling these facts gives us the formula:

$$ N(n, |I|, |A|) := (|A| - 1)!(|I| - |A|)!(n - |I|)! \binom{n}{|I|} \qquad (40) $$

We can simplify this further by expanding the binomial coefficient:

$$ \begin{aligned} N(n, |I|, |A|) &= (|A| - 1)!(|I| - |A|)!(n - |I|)! \binom{n}{|I|} \\ &= \frac{n!(|A| - 1)!(|I| - |A|)!(n - |I|)!}{|I|!(n - |I|)!} \\ &= \frac{n!}{|A|\binom{|I|}{|A|}}. \end{aligned} \qquad (41) $$

This gives us the result. $\square$

*Proof of Lemma A.2.* The proof is quite similar to that of Lemma A.1. We can separate the expectation into two parts:

$$ \mathbb{E}_{\pi \sim U(\Omega)} [f_I(r_{\pi_{:i}}(z, x)_I)] = \underbrace{\frac{1}{n!} \sum_{\substack{\pi \in \Omega \\ I \cap \pi_{:i} = \emptyset}} f_I(r_{\pi_{:i}}(z, x)_I)}_{\text{Term 1}} + \underbrace{\frac{1}{n!} \sum_{\substack{\pi \in \Omega \\ I \cap \pi_{:i} \neq \emptyset}} f_I(r_{\pi_{:i}}(z, x)_I)}_{\text{Term 2}} \qquad (42) $$

**Calculating Term 1.** Notice that $r_{\pi_{:i}}(z, x)_I = z_I$ for all the permutations $\pi$ that satisfy $I \cap \pi_{:i} = \emptyset$, and so the value $f_I(r_{\pi_{:i}}(z, x)_I)$ is constant over all those permutations. It is simply $f_I(r_{\pi_{:i}}(z, x)_I) = f_I(z_I)$. All that we need to do to calculate the number of such permutations. However, a simple symmetry argument gives us the answer. In the proof of Lemma A.1, we showed that the number of permutations $\pi$ such that $I \subseteq \pi_{:i} \cup \{i\}$ is $\frac{n!}{|I|}$. In fact, we can notice that this number should be the same as the number of permutations that satisfy $I \cap \pi_{:i} = \emptyset$. Indeed, in the first case, all variables of $I \setminus \{i\}$ must be to the left of $i$ in $\pi_{:i}$. In the second, they must be to the right. There are exactly as many permutations that satisfy each condition.

**Calculating Term 2.** Here, we consider all the permutations $\pi$ such that $I \cap \pi_{:i} \neq \emptyset$, i.e. the feature subset $I$ shares at least one variable with $\pi_{:i}$. In this situation, some or all variables are replaced when calculating the function. Using the same definition $A := I \cap (\pi_{:i} \cup \{i\})$ as in the proof of Lemma A.1, we have:

$$ r_{\pi_{:i}}(z, x)_I = r_{A \setminus \{i\}}(z, x)_I. $$
The function $f_I(r_{\pi_i}(z, x)_I)$ therefore takes a different value for each possible $A$. Since each $A$ must contain at least two elements $(i, \text{and one more variable so that } I \cap \pi_i \neq \emptyset)$, and also $A = I$ is now acceptable, the intersection $A$ is taken from the set $\mathcal{A}^-(i, I)$, which has the same cardinality as the set $\mathcal{A}^+(i, I)$, namely $2^{|I|-1} - 1$. We also know from the previous proof that there are $N(n, |I|, |A|) = \frac{n!}{|A|(|A|!)}$ permutations $\pi$ such that $A = I \cap (\pi_i \cup \{i\})$. We therefore have:

$$ \text{Term 2} = \sum_{A \in \mathcal{A}^-(i, I)} \frac{f_I(r_{A \setminus \{i\}}(z, x)_I)}{|A|\binom{|I|}{|A|}}, \qquad (43) $$

which concludes the proof. □

## B Calculus for RWRelu expectation

Here we calculate the expectation $E_{w \sim p}[K(u, w)\phi(w, x)]$ for RWRelu (required to calculate the output of the model) through a series of lemmas. First, we recall a pair of lemmas from Dubé & Marchand (2024).

**Lemma B.1 (Lemma 15 of Dubé & Marchand (2024)).** Consider a Hilbert space $\mathcal{W}$. Let $u, w \in \mathcal{W}$ and $a, b > 0$. Then:

$$ \frac{\|w-u\|^2}{a} + \frac{\|w\|^2}{b} = \left(\frac{1}{a} + \frac{1}{b}\right) \left\|w - \frac{1}{1+\frac{a}{b}}u\right\|^2 + \frac{1}{a+b}\|u\|^2. $$

**Lemma B.2 (Lemma 16 of Dubé & Marchand (2024)).** We have:

$$ \int_{\mathbb{R}^n} e^{-a\|w-u\|^2} dw = \left(\frac{\pi}{a}\right)^{\frac{n}{2}}. $$

**Lemma B.3.** We have:

$$ \int_{-\infty}^{\infty} e^{-(w-u)^2/2\gamma^2} \max(0, wx) dw = \frac{\gamma|x|}{\sqrt{2}} \left( \sqrt{2}\gamma e^{-u^2/2\gamma^2} + \sqrt{\pi}u \left( \operatorname{sign}(x) + \operatorname{erf}\left(\frac{u}{\sqrt{2}\gamma}\right) \right) \right). $$

*Proof.* We have:

$$ \begin{aligned} \int_{-\infty}^{\infty} e^{-(w-u)^2/2\gamma^2} \max(0, wx) dw &= \int_{-\infty}^{\infty} e^{-t^2} \max\left(0, \left(\sqrt{2}\gamma t + u\right)x\right) \sqrt{2}\gamma dw & (t := \frac{w-u}{\sqrt{2}\gamma}, \quad dt = \frac{dw}{\sqrt{2}\gamma}) \\ &= \sqrt{2}\gamma \int_{-\infty}^{\infty} e^{-t^2} \max\left(0, \left(\sqrt{2}\gamma t + u\right)x\right) dt. \end{aligned} $$

**Case 1.** If $x \ge 0$, we have:

$$ \begin{aligned} & \sqrt{2}\gamma \int_{-\infty}^{\infty} e^{-t^2} \max\left(0, \left(\sqrt{2}\gamma t + u\right)x\right) dt \\ &= \sqrt{2}\gamma \int_{-u/\sqrt{2}\gamma}^{\infty} e^{-t^2} \left(\sqrt{2}\gamma t + u\right) x dt \\ &= \sqrt{2}\gamma x \left(\sqrt{2}\gamma \int_{-u/\sqrt{2}\gamma}^{\infty} t e^{-t^2} dt + u \int_{-u/\sqrt{2}\gamma}^{\infty} e^{-t^2} dt\right) \\ &= \sqrt{2}\gamma x \left(\sqrt{2}\gamma \left(\frac{-e^{-t^2}}{2}\right) \Big|_{t=-u/\sqrt{2}\gamma}^{t=\infty} + u \left[\int_{-u/\sqrt{2}\gamma}^{0} e^{-t^2} dt + \int_{0}^{\infty} e^{-t^2} dt\right]\right) \\ &= \sqrt{2}\gamma x \left(\frac{\gamma}{\sqrt{2}} e^{-u^2/2\gamma^2} + u \left[\frac{\sqrt{\pi}}{2} \operatorname{erf}\left(\frac{u}{\sqrt{2}\gamma}\right) + \frac{\sqrt{\pi}}{2} \operatorname{erf}(\infty)\right]\right) \\ &= \frac{\gamma x}{\sqrt{2}} \left(\sqrt{2}\gamma e^{-u^2/2\gamma^2} + \sqrt{\pi}u \left[1 + \operatorname{erf}\left(\frac{u}{\sqrt{2}\gamma}\right)\right]\right). \end{aligned} $$
**Case 2.** If $x < 0$, we have instead:

$$
\begin{align*}
& \sqrt{2\gamma} \int_{-\infty}^{\infty} e^{-t^2} \max(0, (\sqrt{2\gamma}t + u)x)dt \\
&= \sqrt{2\gamma} \int_{-\infty}^{-u/\sqrt{2\gamma}} e^{-t^2} (\sqrt{2\gamma}t + u)xdt \\
&= \sqrt{2\gamma}x \left( \sqrt{2\gamma} \int_{-\infty}^{-u/\sqrt{2\gamma}} te^{-t^2}dt + u \int_{-\infty}^{-u/\sqrt{2\gamma}} e^{-t^2}dt \right) \\
&= \sqrt{2\gamma}x \left( \sqrt{2\gamma} \left( \frac{-e^{-t^2}}{2} \Big|_{t=-\infty}^{t=-u/\sqrt{2\gamma}} + u \left[ \int_{-\infty}^{0} e^{-t^2}dt - \int_{-u/\sqrt{2\gamma}}^{0} e^{-t^2}dt \right] \right) \right) \\
&= \sqrt{2\gamma}x \left( \frac{-\gamma}{\sqrt{2}} e^{-u^2/2\gamma^2} + u \left[ \frac{\sqrt{\pi}}{2} \operatorname{erf}(\infty) - \frac{\sqrt{\pi}}{2} \operatorname{erf}\left(\frac{u}{\sqrt{2\gamma}}\right) \right] \right) \\
&= \frac{\gamma x}{\sqrt{2}} \left( -\sqrt{2\gamma}e^{-u^2/2\gamma^2} + \sqrt{\pi}u \left[ 1 - \operatorname{erf}\left(\frac{u}{\sqrt{2\gamma}}\right) \right] \right).
\end{align*}
$$

Both cases can be combined using the sign of $x$, leading to the desired result. □

**Lemma B.4.** We have

$$
\int_{\mathbb{R}^n} e^{-\|w-u\|^2/2\gamma^2} \max(0, \langle w, x \rangle) dw \\
= \left(\sqrt{2\pi\gamma^2}\right)^{n-1} \left[ \frac{\gamma\|x\|}{\sqrt{2}} \left( \sqrt{2\gamma}e^{-\frac{(u,x)^2}{2\gamma^2\|x\|^2}} + \sqrt{\pi}\frac{\langle u, x \rangle}{\|x\|} \left[ 1 + \operatorname{erf}\left(\frac{\langle u, x \rangle}{\sqrt{2\gamma}\|x\|}\right) \right] \right) \right].
$$

*Proof.* Calculate the integral using an orthonormal basis {$v_1, \dots, v_n$} of $\mathbb{R}^n$ such that $v_n := \frac{x}{\|x\|}$. Write $w = (w_1, \dots, w_n)$ in this new basis (i.e. $w_i := \langle w, v_i \rangle$ for all $i$), and similarly $(u_1, \dots, u_n)$ for $u$. Under this change of coordinates, the integral becomes:

$$
\int_{\mathbb{R}^n} e^{-\|w-u\|^2/2\gamma^2} \max(0, \langle w, x \rangle) dw \\
= \int_{\mathbb{R}} \int_{\mathbb{R}^{n-1}} e^{-[\sum_{i=1}^{n-1} (w_i - u_i)^2 + (w_n - u_n)^2]/2\gamma^2} \max(0, u_n \|x\|) dw_1 \dots dw_n.
$$

We are left with a product of $n$ independent integrals:

$$
\begin{align*}
& \int_{\mathbb{R}^n} e^{-\|w-u\|^2/2\gamma^2} \max(0, \langle w, x \rangle) dw \\
&= \int_{\mathbb{R}^{n-1}} e^{-\sum_{i=1}^{n-1} (w_i - u_i)^2/2\gamma^2} dw_1 \dots dw_{n-1} \int_{\mathbb{R}} e^{-\|w_n - u_n\|^2/2\gamma^2} \max(0, w_n \|x\|) dw_n \\
&= \int_{\mathbb{R}^{n-1}} \prod_{i=1}^{n-1} e^{-\|w_i - u_i\|^2/2\gamma^2} dw_1 \dots dw_{n-1} \int_{\mathbb{R}} e^{-\|w_n - u_n\|^2/2\gamma^2} \max(0, w_n \|x\|) dw_n \\
&= \prod_{i=1}^{n-1} \int_{-\infty}^{\infty} e^{-\|w_i - u_i\|^2/2\gamma^2} dw_i \int_{\mathbb{R}} e^{-\|w_n - u_n\|^2/2\gamma^2} \max(0, w_n \|x\|) dw_n.
\end{align*}
$$

For each $i$, Lemma B.2 gives us:

$$
\int_{\mathbb{R}} e^{-(w_i - u_i)^2 / 2\gamma^2} dw_i = \sqrt{2\pi\gamma^2}.
$$
Also, Lemma B.3 gives us:

$$ \int_{\mathbb{R}} e^{-(w_n-u_n)^2/2\gamma^2} \max(0, w_n \|x\|) dw_n = \frac{\gamma \|x\|}{\sqrt{2}} \left( \sqrt{2\gamma} e^{-u_n^2/2\gamma^2} + \sqrt{\pi} u_n \left( \operatorname{sign}(\|x\|) + \operatorname{erf}\left(\frac{u_n}{\sqrt{2\gamma}}\right) \right) \right). $$

Finally, since $u_n = \frac{\langle u, x \rangle}{\|x\|}$ and $\operatorname{sign}(\|x\|)$, we have the result. □

*Proof of Theorem 4.1.* We have:

$$ \begin{aligned} \mathbb{E}_{w \sim p} [\mathcal{K}(u, w)\phi(w, x)] &= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^n \int_{\mathbb{R}^n} e^{-\|w-u\|^2/2\gamma^2} e^{-\|w\|^2/2\sigma^2} \max(0, \langle w, x \rangle) dw \\ &= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^n e^{-\frac{\|u\|^2}{2\sigma^2+2\gamma^2}} \int_{\mathbb{R}^n} e^{-\left(\frac{1}{2\gamma^2} + \frac{1}{2\sigma^2}\right) \left\|w - \frac{u}{1+\frac{\gamma^2}{\sigma^2}}\right\|^2} \max(0, \langle w, x \rangle) dw \quad (\text{Lemma B.1}) \\ &= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^n e^{-\frac{\|u\|^2}{2\sigma^2+2\gamma^2}} \int_{\mathbb{R}^n} e^{-\|w-u'\|^2/2\zeta^2} \max(0, \langle w, x \rangle) dw. \end{aligned} $$

Applying Lemma B.4, we obtain:

$$ \begin{aligned} \mathbb{E}_{w \sim p} [\mathcal{K}(u, w)\phi(w, x)] &= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^n e^{\frac{-\|u\|^2}{2\sigma^2+2\gamma^2}} (\sqrt{2\pi\zeta})^{n-1} \left[ \frac{\zeta\|x\|}{\sqrt{2}} \left( \sqrt{2\zeta}e^{-\frac{\langle u', x \rangle^2}{2\zeta^2\|x\|^2}} + \sqrt{\pi}\frac{\langle u', x \rangle}{\|x\|} \left[ 1 + \operatorname{erf}\left(\frac{\langle u', x \rangle}{\sqrt{2\zeta}\|x\|}\right) \right] \right) \right] \\ &= \left(1 + \frac{\sigma^2}{\gamma^2}\right)^{-n/2} e^{\frac{-\|u\|^2}{2\sigma^2+2\gamma^2}} \frac{\|x\|}{2\sqrt{\pi}} \left( \sqrt{2\zeta}e^{-\frac{\langle u', x \rangle^2}{2\zeta^2\|x\|^2}} + \sqrt{\pi}\frac{\langle u', x \rangle}{\|x\|} \left[ 1 + \operatorname{erf}\left(\frac{\langle u', x \rangle}{\sqrt{2\zeta}\|x\|}\right) \right] \right). \end{aligned} $$

<figure><img src="image_7.png" /></figure>

## C Theoretical constants of RWRelu

For each of the two instantiations Dubé & Marchand (2024) introduce, they provide the value of two theoretical constants, $\kappa$ and $\theta$, relevant for their theoretical guarantees. For a given instantiation $(W, \phi, K, p)$, we have:

$$ \kappa := \sup_{x \in \mathcal{X}} \sqrt{\mathbb{E}_{w \sim p} \left[ \|\phi(w, x)\mathcal{K}(w, \cdot)\|_{\mathcal{H}}^2 \right]} = \sup_{x \in \mathcal{X}} \sqrt{\mathbb{E}_{w \sim p} \left[ \mathcal{K}(w, w)\phi(w, x)^2 \right]} \quad (44) $$

$$ \theta := \sup_{x \in \mathcal{X}} \left\| \mathbb{E}_{w \sim p} [\phi(w, x) \mathcal{K}(w, \cdot)] \right\|_{\mathcal{H}} = \sup_{x \in \mathcal{X}} \sqrt{\mathbb{E}_{w \sim p} \mathbb{E}_{u \sim p} [\mathcal{K}(u, w) \phi(u, x) \phi(w, x)]}. \quad (45) $$

For completeness, we give here the exact value of $\kappa$ and an upper bound for $\theta$, and the relevant calculus, in the case of instantiation RWRelu. We begin by a couple of lemmas to help with the calculus.

**Lemma C.1.** Consider $\sigma > 0$ and $\gamma > 0$. Then:

$$ \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\|u-w\|^2/2\gamma^2} e^{-\langle u,w \rangle/\sigma^2} dw dw = (2\pi)^n \left( \frac{\gamma^2\sigma^4}{2\sigma^2 - \gamma^2} \right)^{n/2}. \quad (46) $$
Proof.

$$
\begin{align*}
\int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|u-w\|^2}{2\gamma^2}} e^{-\frac{\langle u,w \rangle}{\sigma^2}} \mathrm{d}udw &= \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|t\|^2}{2\gamma^2}} e^{-\frac{\langle t+w,w \rangle}{\sigma^2}} \mathrm{d}tdw \quad (t := u-w, dt = du) \\
&= \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|t\|^2}{2\gamma^2}} e^{-\frac{\|w\|^2}{\sigma^2}} e^{-\frac{\langle t,w \rangle}{\sigma^2}} \mathrm{d}tdw \\
&= \int_{\mathbb{R}^n} e^{-\frac{\|w\|^2}{\sigma^2}} \left[ \int_{\mathbb{R}^n} e^{-\frac{\|t\|^2}{2\gamma^2}} e^{-\frac{\langle t,w \rangle}{\sigma^2}} dt \right] dw \\
&= \int_{\mathbb{R}^n} e^{-\frac{\|w\|^2}{\sigma^2}} \left[ \int_{\mathbb{R}^n} e^{-\frac{\left\| t+\frac{\gamma^2 w}{\sigma^2} \right\|^2}{2\gamma^2}} e^{-\frac{\left\| \frac{\gamma^2 w}{\sigma^2} \right\|^2}{2\gamma^2}} dt \right] dw \\
&= \int_{\mathbb{R}^n} e^{-\frac{\|w\|^2}{\sigma^2}} e^{-\frac{\left\| \frac{\gamma^2 w}{\sigma^2} \right\|^2}{2\gamma^2}} \left[ \int_{\mathbb{R}^n} e^{-\frac{\left\| t+\frac{\gamma^2 w}{\sigma^2} \right\|^2}{2\gamma^2}} dt \right] dw \\
&= \left(\sqrt{2\pi\gamma^2}\right)^n \int_{\mathbb{R}^n} e^{-\frac{\|w\|^2}{\sigma^2}} e^{-\frac{\left\| \frac{\gamma^2 w}{\sigma^2} \right\|^2}{2\gamma^2}} dw.
\end{align*}
$$

Then, simplifying the exponent:

$$
\begin{align*}
-\frac{\|w\|^2}{\sigma^2} + \frac{\left\|\frac{\gamma^2 w}{\sigma^2}\right\|^2}{2\gamma^2} &= -\frac{\|w\|^2}{2\sigma^2} \left(2 - \frac{\gamma^2}{\sigma^2}\right) \\
&= -\frac{\|w\|^2}{2\sigma^2} \left(\frac{2\sigma^2 - \gamma^2}{\sigma^2}\right) \\
&= -\frac{\|w\|^2}{2\sigma^4} (2\sigma^2 - \gamma^2),
\end{align*}
$$

we get:

$$
\begin{align*}
\int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|u-w\|^2}{2\gamma^2}} e^{-\frac{\langle u,w \rangle}{\sigma^2}} \mathrm{d}udw &= \left(\sqrt{2\pi\gamma^2}\right)^n \int_{\mathbb{R}^n} e^{-\frac{\|w\|^2}{2\sigma^4}(2\sigma^2-\gamma^2)} \mathrm{d}w \\
&= \left(\sqrt{2\pi\gamma^2}\right)^n \left(\sqrt{2\pi \frac{\sigma^4}{2\sigma^2 - \gamma^2}}\right)^n \\
&= (2\pi)^n \left(\frac{\gamma^2 \sigma^4}{2\sigma^2 - \gamma^2}\right)^{n/2}.
\end{align*}
$$

□

**Lemma C.2.** Consider $\sigma > 0$ and $\gamma > 0$. Denote $I_n$ the identity matrix in $\mathbb{R}^n$. Then:

$$
E_{w \sim \mathcal{N}(0, \sigma^2 I_n)} E_{u \sim \mathcal{N}(0, \sigma^2 I_n)} \left[ e^{-\|u-w\|^2 / 2\gamma^2} \right] = \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-n/2}. \quad (47)
$$
Proof. The expectation is a straightforward integral:

$$
\begin{align*}
\mathbb{E}_{w \sim p, u \sim p} \mathbb{E}[\mathcal{K}(u, w)] &= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^{2n} \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|u-w\|^2}{2\gamma^2}} e^{-\frac{\|u\|^2}{2\sigma^2}} e^{-\frac{\|w\|^2}{2\sigma^2}} \mathrm{d}uw \\
&= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^{2n} \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|u-w\|^2}{2\gamma^2}} e^{-\frac{\|u\|^2}{2\sigma^2}} e^{\frac{\langle u, w \rangle}{\sigma^2}} e^{-\frac{\|w\|^2}{2\sigma^2}} e^{-\frac{\langle u, w \rangle}{\sigma^2}} \mathrm{d}uw \\
&= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^{2n} \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\frac{\|u-w\|^2}{2\zeta^2}} e^{-\frac{\langle u, w \rangle}{\sigma^2}} \mathrm{d}uw \quad (\frac{1}{2\zeta^2} = \frac{1}{2\gamma^2} + \frac{1}{2\sigma^2} = \frac{\sigma^2+\gamma^2}{2\sigma^2\gamma^2}) \\
&= \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^{2n} (2\pi)^n \left(\frac{\zeta^2\sigma^4}{2\sigma^2 - \zeta^2}\right)^{n/2} \quad (\text{Lemma C.1}) \\
&= \left(\frac{\zeta^2}{2\sigma^2 - \zeta^2}\right)^{n/2} \\
&= \left(\frac{2\sigma^2}{\zeta^2} - 1\right)^{-n/2} \\
&= \left(\frac{2\sigma^2}{\sigma^2} + \frac{2\sigma^2}{\gamma^2} - 1\right)^{-n/2} \\
&= \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-n/2}.
\end{align*}
$$

<figure><img src="image_8.png" /></figure>

Lemma C.3. Considering RWRelu, we have:

$$ \kappa = \frac{\sigma}{\sqrt{2}} \sup_{x \in \mathcal{X}} \|x\|. $$

Proof. We have:

$$
\begin{align*}
\kappa^2 &:= \sup_{x \in \mathcal{X}} \mathbb{E}_{w \sim p} [\mathcal{K}(w, w) \phi(w, x)^2] \\
&= \sup_{x \in \mathcal{X}} \mathbb{E}_{w \sim p} \left[ \max(0, \langle w, x \rangle)^2 \right] \\
&= \sup_{x \in \mathcal{X}} \left( \frac{1}{\sqrt{2\pi\sigma^2}} \right)^n \int_{\mathbb{R}^n} e^{-\|w\|^2/2\sigma^2} \max(0, \langle w, x \rangle)^2 dw \\
&= \sup_{x \in \mathcal{X}} \frac{1}{\sqrt{2\pi\sigma^2}} \int_{\mathbb{R}} e^{-w_n^2/2\sigma^2} \max(0, w_n \|x\|)^2 dw_n \quad (w_n := \frac{\langle w, x \rangle}{\|x\|}, \text{ see proof of Lemma B.4}) \\
&= \sup_{x \in \mathcal{X}} \frac{1}{\sqrt{2\pi\sigma^2}} \int_0^\infty e^{-w_n^2/2\sigma^2} w_n^2 \|x\|^2 dw_n \\
&= \sup_{x \in \mathcal{X}} \|x\|^2 \frac{1}{\sqrt{2\pi\sigma^2}} \int_0^\infty e^{-w_n^2/2\sigma^2} w_n^2 dw_n \\
&= \frac{1}{2} \sup_{x \in \mathcal{X}} \|x\|^2 \mathbb{E}_{w_n \sim \mathcal{N}(0, \sigma^2)} [w_n^2] \\
&= \frac{\sigma^2}{2} \sup_{x \in \mathcal{X}} \|x\|^2.
\end{align*}
$$

We have the result by taking the square root.

<figure><img src="image_9.png" /></figure>
**Lemma C.4.** *Considering RWRelu, we have:*

$$ \theta^2 \leq \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{\sigma^2}{2\pi}\right). \quad (48) $$

*Proof.* We have:

$$ \begin{aligned} \theta^2 &= \sup_{x \in \mathcal{X}} \mathbb{E}_{u \sim p, w \sim p} \mathbb{E} [K(u, w)\phi(u, x)\phi(w, x)] \\ &= \sup_{x \in \mathcal{X}} \left(\frac{1}{2\pi\sigma^2}\right)^n \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\|u-w\|^2/2\gamma^2} e^{-\|u\|^2/2\sigma^2} e^{-\|w\|^2/2\sigma^2} \max(0, \langle u, x \rangle) \max(0, \langle w, x \rangle) \mathrm{d}u\mathrm{d}w. \end{aligned} $$

We are free to solve these integrals using any orthonormal basis of $\mathbb{R}^n$. We choose a basis $\{v_1, \dots, v_n\}$ such that:

$$ v_n := \frac{x}{\|x\|}. \quad (49) $$

In particular, we have:

$$ u_n := \left\langle u, \frac{x}{\|x\|} \right\rangle, $$

$$ w_n := \left\langle w, \frac{x}{\|x\|} \right\rangle, $$

$$ \langle u, x \rangle = u_n \|x\|, $$

$$ \langle w, x \rangle = w_n \|x\|. $$

Denoting $u_{1:n-1} := (u_1, \dots, u_{n-1})$ and $w_{1:n-1} := (w_1, \dots, w_{n-1})$, we can rewrite the expression:

$$ e^{-\|u-w\|^2/2\gamma^2} e^{-\|u\|^2/2\sigma^2} e^{-\|w\|^2/2\sigma^2} \quad (50) $$

as:

$$ e^{-\|u_{1:n-1} - w_{1:n-1}\|^2/2\gamma^2} e^{-\|u_{1:n-1}\|^2/2\sigma^2} e^{-\|w_{1:n-1}\|^2/2\sigma^2} e^{-(u_n - w_n)^2/2\gamma^2} e^{-u_n^2/2\sigma^2} e^{-w_n^2/2\sigma^2}. \quad (51) $$

This allows us to separate the integral into two parts. The first one is the integral over $(u_1, \dots, u_{n-1})$ and $(w_1, \dots, w_{n-1})$:

$$ \left(\frac{1}{2\pi\sigma^2}\right)^{n-1} \int_{\mathbb{R}^{n-1}} \int_{\mathbb{R}^{n-1}} e^{-\|u_{1:n-1} - w_{1:n-1}\|^2/2\gamma^2} e^{-\|u_{1:n-1}\|^2/2\sigma^2} e^{-\|w_{1:n-1}\|^2/2\sigma^2} \mathrm{d}u_{1:n-1} \mathrm{d}w_{1:n-1}. \quad (52) $$

Lemma C.2 tells us that the previous expression is equal to:

$$ \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2}. \quad (53) $$

The second integral is the one over $u_n$ and $w_n$:

$$ \left(\frac{1}{2\pi\sigma^2}\right) \int_{\mathbb{R}} \int_{\mathbb{R}} e^{-(u_n - w_n)^2/2\gamma^2} e^{-u_n^2/2\sigma^2} e^{-w_n^2/2\sigma^2} \max(0, u_n \|x\|) \max(0, w_n \|x\|) \mathrm{d}u_n \mathrm{d}w_n. \quad (54) $$

We can simplify the integral by noticing that $\max(0, u_n) = u_n$ if $u_n \ge 0$, and 0 otherwise:

$$ \left(\frac{1}{2\pi\sigma^2}\right) \|x\|^2 \int_0^\infty \int_0^\infty e^{-(u_n - w_n)^2/2\gamma^2} e^{-u_n^2/2\sigma^2} e^{-w_n^2/2\sigma^2} u_n w_n \mathrm{d}u_n \mathrm{d}w_n. \quad (55) $$

The expression for $\theta^2$ has now become:

$$ \begin{aligned} \theta^2 &= \sup_{x \in \mathcal{X}} \left(\frac{1}{2\pi\sigma^2}\right)^n \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} e^{-\|u-w\|^2/2\gamma^2} e^{-\|u\|^2/2\sigma^2} e^{-\|w\|^2/2\sigma^2} \max(0, \langle u, x \rangle) \max(0, \langle w, x \rangle) \mathrm{d}u\mathrm{d}w \\ &= \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{1}{2\pi\sigma^2}\right) \int_0^\infty \int_0^\infty e^{-(u_n - w_n)^2/2\gamma^2} e^{-u_n^2/2\sigma^2} e^{-w_n^2/2\sigma^2} u_n w_n \mathrm{d}u_n \mathrm{d}w_n. \end{aligned} \quad (56) $$
The remaining integral is difficult. We can simplify by using the fact that $e^{-(u_n-w_n)^2/2\gamma^2} \le 1$, and obtain an upper bound on $\theta^2$:

$$ \begin{aligned} \theta^2 &= \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{1}{2\pi\sigma^2}\right) \int_0^\infty \int_0^\infty e^{-(u_n-w_n)^2/2\gamma^2} e^{-u_n^2/2\sigma^2} e^{-w_n^2/2\sigma^2} u_n w_n du_n dw_n \\ &\le \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{1}{2\pi\sigma^2}\right) \int_0^\infty \int_0^\infty e^{-u_n^2/2\sigma^2} e^{-w_n^2/2\sigma^2} u_n w_n du_n dw_n \\ &= \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{1}{2\pi\sigma^2}\right) \left(\int_0^\infty e^{-w_n^2/2\sigma^2} w_n dw_n\right)^2 \\ &= \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{1}{2\pi\sigma^2}\right) \left(\left(-\sigma^2 e^{-w_n^2/2\sigma^2}\right) \Big|_{w_n=0}^\infty\right)^2 \\ &= \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{1}{2\pi\sigma^2}\right) \sigma^4 \\ &= \sup_{x \in \mathcal{X}} \|x\|^2 \left(1 + \frac{2\sigma^2}{\gamma^2}\right)^{-(n-1)/2} \left(\frac{\sigma^2}{2\pi}\right). \end{aligned} $$

<figure><img src="image_10.png" /></figure>

## D Theoretical constants of STAR RKHS Weightings

We showed in Section 4 how to transform an RKHS Weighting instantiation into a STAR model, while maintaining the form of the model as an RKHS Weighting (Table 4 contains the details). In particular, all algorithms and theoretical guarantees from Dubé & Marchand (2024) still apply. The guarantees involve the theoretical constants:

$$ \kappa := \sup_{x \in \mathcal{X}} \sqrt{\mathbb{E}_{w \sim p} \left[ \|\phi(w, x) \mathcal{K}(w, \cdot)\|_{\mathcal{H}}^2 \right]} = \sup_{x \in \mathcal{X}} \sqrt{\mathbb{E}_{w \sim p} \left[ \mathcal{K}(w, w) \phi(w, x)^2 \right]} \quad (57) $$

$$ \theta := \sup_{x \in \mathcal{X}} \left\| \mathbb{E}_{w \sim p} [\phi(w, x) \mathcal{K}(w, \cdot)] \right\|_{\mathcal{H}} = \sup_{x \in \mathcal{X}} \sqrt{\mathbb{E}_{w \sim p} \mathbb{E}_{u \sim p} [\mathcal{K}(u, w) \phi(u, x) \phi(w, x)]}, \quad (58) $$

defined for any instantiation $(\mathcal{W}, \phi, \mathcal{K}, p)$.

Here, we explain how to obtain those constants, if desired, for a STAR RKHS Weighting instantiation. For $\kappa$, we have:

$$ \begin{aligned} \kappa^2 &= \sup_{x \in \mathcal{X}} \mathbb{E}_{(w, I) \sim p_{\mathcal{W}|I} \times p_{[n]}} [\mathcal{K}((w, I), (w, I)) \phi_I(w, x_I)^2] \\ &= \sup_{x \in \mathcal{X}} \mathbb{E}_{I \sim p_{[n]}} \mathbb{E}_{w \sim p_{\mathcal{W}|I}} [\mathcal{K}_W(w, w) 1[I=I] \phi_I(w, x_I)^2] \\ &= \sup_{x \in \mathcal{X}} \mathbb{E}_{I \sim p_{[n]}} \mathbb{E}_{w \sim p_{\mathcal{W}|I}} [\mathcal{K}_W(w, w) \phi_I(w, x_I)^2] \\ &\le \mathbb{E}_{I \sim p_{[n]}} \sup_{x \in \mathcal{X}} \mathbb{E}_{w \sim p_{\mathcal{W}|I}} [\mathcal{K}_W(w, w) \phi_I(w, x_I)^2] \end{aligned} $$

The expression $\sup_{x \in \mathcal{X}} \mathbb{E}_{w \sim p_{\mathcal{W}|I}} [\mathcal{K}_W(w, w) \phi_I(w, x_I)^2]$ is in fact the constant $\kappa^2$ for the partial instantiation defined only on the feature subset $I$. Let us denote it $\kappa_I^2$. Then we have:

$$ \kappa^2 \le \mathbb{E}_{I \sim p_{[n]}} \kappa_I^2. \quad (59) $$

This quantity depends on the distribution of the feature subsets $p_{[n]}$. In the simplest cases, where for example $p_{[n]}$ is the uniform distribution over all feature subsets of size $k$, then the constant $\kappa_I$ might simply be the
same for all $I$. In the worst cases, we can take the supremum of $\kappa_I$ over all possible subsets $I$ instead of the expectation. As for $\theta$:

$$
\begin{align*}
\theta^2 &= \sup_{x \in \mathcal{X}} \mathbb{E}_{(w, I) \sim p_{\mathcal{W}|I} \times p_{[n]} (u, J) \sim p_{\mathcal{W}|J} \times p_{[n]}} [\mathcal{K}((u, J), (w, I))\phi_J(u, x_J)\phi_I(w, x_I)] \\
&= \sup_{x \in \mathcal{X}} \mathbb{E}_{(w, I) \sim p_{\mathcal{W}|I} \times p_{[n]} (u, J) \sim p_{\mathcal{W}|J} \times p_{[n]}} [\mathcal{K}_{\mathcal{W}}(u, w) \mathbb{1}[I = J]\phi_J(u, x_J)\phi_I(w, x_I)] \\
&= \sup_{x \in \mathcal{X}} \mathbb{E}_{(w, I) \sim p_{\mathcal{W}|I} \times p_{[n]} (u \sim p_{\mathcal{W}|I})} [p_{[n]}(I)\mathcal{K}_{\mathcal{W}}(u, w)\phi_I(u, x_I)\phi_I(w, x_I)] \\
&= \sup_{x \in \mathcal{X}} \mathbb{E}_{I \sim p_{[n]}} p_{[n]}(I) \mathbb{E}_{w \sim p_{\mathcal{W}|I} \times p_{[n]} (u \sim p_{\mathcal{W}|I})} [\mathcal{K}_{\mathcal{W}}(u, w)\phi_I(u, x_I)\phi_I(w, x_I)] \\
&\le \sup_{I \sim p_{[n]}} \mathbb{E}_{I \sim p_{[n]}} p_{[n]}(I) \sup_{x \in \mathcal{X}} \mathbb{E}_{w \sim p_{\mathcal{W}|I} \times p_{[n]} (u \sim p_{\mathcal{W}|I})} [\mathcal{K}_{\mathcal{W}}(u, w)\phi_I(u, x_I)\phi_I(w, x_I)].
\end{align*}
$$

Exactly as before, the expression $\sup_{x \in \mathcal{X}} \mathbb{E}_{w \sim p_{\mathcal{W}|I}} \mathbb{E}_{u \sim p_{\mathcal{W}|I}} [p_{[n]}(I) \mathcal{K}_{\mathcal{W}}(u, w) \phi_I(u, x_I) \phi_I(w, x_I)]$ is the constant $\theta^2$ for the partial model defined on the subset $I$. We can write it as $\theta_I^2$ and:

$$ \theta^2 \leq \mathbb{E}_{I \sim p_{[n]}} p_{[n]}(I) \theta_I^2. \qquad (60) $$

In the event that $\theta_I \le \hat{\theta}$ for all $I$ for some $\hat{\theta}$, and that $p_{[n]}(I)$ is constant, then we have:

$$ \theta^2 \le p_{[n]}(I) \hat{\theta}^2, \qquad (61) $$

or:

$$ \theta \le \sqrt{p_{[n]}(I)\hat{\theta}}. \qquad (62) $$

# E Details of experimentation

## Preprocessing of the datasets

All datasets have been scaled to have mean 0 and standard deviation 1 on all variables, including the target labels. Means and standard variations were calculated on the training data, then the transformation applied to both training and test datasets.

## Numerical stability of the learning algorithm

To learn RKHS Weightings, we used Algorithm 3 of Dubé & Marchand (2024), the least squares fit of the coefficients. However, we have run into some rare numerical instability problems which led to failure of learning and even negative training $R^2$ scores, which is impossible in theory (at worst, the algorithm could output $0 \in \mathcal{H}$, which has an $R^2$ of 0 exactly). To solve this issue, we added a very small $\ell^2$ regularizer to the optimization objective (Equation 53 of Dubé & Marchand (2024)), which is now:

$$ \mathcal{L}_{S}^{\text{reg}}(\Lambda\alpha(x)) = \frac{1}{2m} \|\Phi a - y\|_{2}^{2} + \frac{\lambda}{2} a^{\top} G a + \frac{\epsilon}{2} a^{\top} I a, \qquad (63) $$

where $I$ is the identity matrix, and $\epsilon$ is the $\ell^2$ regularization parameter, which we have set to $10^{-10}$. The minimizer $a$ of the previous expression is the solution to the linear problem:

$$ (\Phi^{\top}\Phi + m\lambda G + m\epsilon I)a = \Phi^{\top}y. \qquad (64) $$

## Hyperparameter selection

As pointed out in Dubé & Marchand (2024), instantiation RWSign, and now also RWRelu, exhibits exponential behavior in the dimensionality $n$ of the instance space. This can immediately be seen in the coefficient in Equation (25), shared by both RWSign and RWRelu. Instead of using $\sigma$
<table><caption>Table 9: Algorithms and models used in this paper and their hyperparameters.</caption><thead><tr><th></th><th>Cross-validation parameters</th><th>Source code (clickable)</th></tr></thead><tbody><tr><td>DecisionTreeRegressor</td><td>max depth ∈ {2, 5, 10, 20}</td><td>Scikit-learn</td></tr><tr><td>EBM</td><td>max bins ∈ {512, 1024, 2048}<br>learning rate ∈ {512, 1024, 2048}<br>max rounds ∈ {15000, 25000, 35000}<br>min samples leaf ∈ {1, 2, 3}</td><td>InterpretML</td></tr><tr><td>KernelRidge</td><td>kernel = rbf<br>alpha ∈ {0.01, 0.05, 0.1, 0.5, 1, 5}</td><td>Scikit-learn</td></tr><tr><td>LinearRegression</td><td></td><td>Scikit-learn</td></tr><tr><td>SVR</td><td>C ∈ {0.5, 1, 5, 10, 50}</td><td>Scikit-learn</td></tr><tr><td>RWSign</td><td>max theta ∈ {0.1, 0.5, 0.9}<br>λ ∈ {10<sup>-9</sup>, 10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-5</sup>, 10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>}</td><td>Dubé & Marchand (2024)</td></tr><tr><td>RWRelu</td><td>max theta ∈ {0.1, 0.5, 0.9}<br>λ ∈ {10<sup>-9</sup>, 10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-5</sup>, 10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>}</td><td>This paper</td></tr><tr><td>RWStumps</td><td>σ ∈ {0.01, 0.1, 1}<br>γ ∈ {0.01, 0.1, 1}<br>λ ∈ {10<sup>-9</sup>, 10<sup>-8</sup>, 10<sup>-7</sup>, 10<sup>-6</sup>, 10<sup>-5</sup>, 10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>}</td><td>Dubé & Marchand (2024)</td></tr></tbody></table>

and $\gamma$ as hyperparameters of the model, they use $\sigma$ and $c$, and calculate $\gamma$ from $\sigma$ and $c$. This way, the effect of dimensionality is held constant regardless of $n$.

We modify this methodology slightly. We instead use $c^2 := (1 + \frac{2\sigma^2}{\gamma^2})^{-n/2}$. This is the upper bound for $\theta^2$ calculated by Dubé & Marchand (2024) for RWSign. It can also be found in Lemma C.2, and almost as is in the upper bound for $\theta$ for RWRelu (Lemma C.4). Given $\sigma$ and $c \in (0, 1)$, we get:

$$ \gamma^2 = \frac{2}{\sigma^2(c^{-4/n} - 1)}. \qquad (65) $$

The advantage of this method is that the model output $\Lambda\alpha(x)$ being upper bounded by $\theta\|\alpha\|_{\mathcal{H}} \le c\|\alpha\|_{\mathcal{H}}$, this gives meaning to the parameter $c$. For RWSign, this allows setting the maximum value of the model. For RWRelu, the relationship is functionally the same, though the upper bound for $\theta$ is slightly more complex.

Finally, Table 9 contains the hyperparameters of all the algorithms and models used in this paper.

## Reproducing the results

The code for this paper can be found at https://github.com/gadub44/star-rkhs-weightings. The repository contains a `requirements.txt` file containing the particular Python packages that were used in our experiments. They can be installed from the command line with the command `pip install -r requirements.txt`. Table 10 lists the commands to run the experiments of Section 5, and the files they generate.
<table><caption>Table 10: Commands to run the various experiments in this paper.</caption><thead><tr><th></th><th>Command</th><th>Files generated</th></tr></thead><tbody><tr><td>Fig. 1</td><td>python -m major_experiments.shapley_time</td><td>results/shapley_time.csv<br>figures/shapley_time_left.pdf<br>figures/shapley_time_right.pdf</td></tr><tr><td>Fig. 2</td><td>python -m major_experiments.shapleyendor</td><td>results/shapleyendor accuracy-n10-T100.csv<br>results/shapleyendor accuracy-n100-T100.csv<br>results/shapleyendor accuracy-n500-T100.csv<br>figures/shapleyendor accuracy-n10-T100.pdf<br>figures/shapleyendor accuracy-n100-T100.pdf<br>figures/shapleyendor accuracy-n500-T100.pdf</td></tr><tr><td>Table 6</td><td>python -m major_experiments.regression --final</td><td>results/regression-final.csv<br>tables/regression-final.tex<br>tables/regression-final-dense.tex</td></tr><tr><td>Table 8</td><td>python -m major_experiments.time_series --final</td><td>results/regression-final.csv<br>tables/time-series-final.tex<br>tables/time-series-final-dense.tex</td></tr><tr><td>Fig. 3, 4</td><td>python -m major_experiments.shapley_comparison</td><td>results/ebmExplanation.pkl<br>results/rkhsExplanation.pkl<br>figures/beeswarm plot_ebm.pdf<br>figures/beeswarm plot_rkhs.pdf<br>figures/waterfall plot_ebm.pdf<br>figures/waterfall plot_rkhs.pdf</td></tr></tbody></table>