Shapley Values of Structured Additive Regression Models and Application to RKHS Weightings of Functions
Authors: Gabriel Dubé
Mario Marchand
Venue: Accepted by TMLR ()

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.

Published in Transactions on Machine Learning Research (01/2025)
Shapley Values of Structured Additive Regression Models and Application to RKHS Weightings of Functions
Gabriel Dubé
gadub44 @rulaval.c
Département d'informatique et de génie logiciel
Université Laval
Mario Marchand
mario.marchaand@@....oooo
Département d' d'informatique et de génie logiciel
Université Loval
Reviewed on OpenReview: // openreiee.. net/ forum? id= aWRMuXTUP
Abstract
Shapley values are widely used in machine learning to interpret model predictions. However, they have an important drawback their computational time, which is exponential 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 familics are fairly present STAR-SHAP, algorithm for efficiently calculating
restrictedd Consequently, 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, 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 new RKHS Weighting instantiation, and show how to transform and other RKHS Weightings into STAR models. We therefore introduce new family of STAR models, as well as the means to interpret their outputs in a timely manner.
1 Introduction
The ability to interpret machine learning model is 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 (GAM-)) Hastie (017)) or Decision Trees (Kotsiantis) 2013), are inherently easier to understand, but do not always reach the same levels of performance (though recent work appears show otherwise (Rudin et al. 2024)). Work improve explanation methods, especially the interpretation of more complex models, can help bridge this gap.
Informally, we mean by interpreting model the task of explaining its output. In sensitive contexts such healthcare (Caruana et al. 2015), energy systems (Zhang & Chen 2024h 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 AT model 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
TTbble
lists all the acronyms used this paper.

---

Published in Transactions on Machine Learning Research (01/2025)
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 simple predictor, and Partial Dependance Plots (Friedman 2001), which visualize the effect of 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 cach variable.
With all of this in mind, the purpose of this paper is twofold. First, we present an algorithm for efficiently caleulating the Shapley values of large family of prediction models, Structured Additive Regression (STAR) models. STAR models are 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
A presents a new RKHS Weighting instantiation, and shows how to transform an RKHS Weighting model into a STAR model. Section 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 go over some limitations of this paper as well future work to be done, and conclude in Section
3 contains various experiments, showing the computation
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 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 X, usually X' = R", the label space y = R, and a dataset of examples S S= ((),,)))) C X'x x y, our task is to find a predictor (also called model) h X + y' that optimizes some objective function of the sample data. To measure the quality of a model, use loss function 2: y x y + R, which in this paper is the squared loss:
e(h(z).)) == (h(zz) - -
(1)
2
2

---

Published in Transactions on Machine Learning Research (01/2025)
Acronym
Meaning
Relevant link
GAM
Generalized Additive Model
Hastic
2017)
GA2M
GAM with pairwise interactions Structured Additive Regression (model)
Lou et al. [2013)
STAR (model)
Equation 16
k-STAR model
STAR model using only k-way interactions
RKHS
Reproducing Kernel Hilbert Space
Berlinct & Thomas Agomas-Agnan
2011
MSE
Mean Square Error
Equation
R²
Coefficient of determination
Equation
SHAP
SHapley Additive exPlanations
Lundberg & Lee
2017
SHAP-IQ
SHAPley Interaction Quantification Structured Additive Regression SHAP
Fumagalli et al. 2023
STAR-SHAP
Algorithm
erf
Error function, er((x) = ア So e d/ Rectified Linear Unit, ReLU(=) == max(0,
ReLU
RWSign
RKHS Weighting using sign function as base predictor
I1 in Dubé & Marchand
2024
RWStumps
RKHS Weighting using decisions stumps as base predictor
T2 in Dubé & Marchand Section Lou et al. 2013
2024
RWReLU
RKHS Weighting using ReLU as base predictor
EBM
Explainable Boosting Machine
SVR
Support Vector Regression
Pedregosa el. al. (2011
Table 1: Acronymss used in this paper.
The empirical risk is the sample average loss. Square Error (MSE):
When using the squared loss, this is called the Mean
MSEs
Σ (h()))
-1)
2
2m
The MSE has units of the labels squared. To successfully compare experimental results on different datascts, we instead report the coefficient of determination:
E
- (y h(xx))² m 1= (yi 7)
Rg(h) ;=
(3)
where 5 1 Em Mi 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 D over X y. The true risk, or true MSE, is the expected squared loss over this distribution:
(SE = E (u.N)D
[(((z) y)
(4)
We estimate this value by calculating the MSE set of examples distinet from the sample.
A learning algorithm A is function that takes in sample S, and returns predictor, i.c. A(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 algorithm. Suppose we have combinations to choose from. This implies have the candidate algorithms Ai to Ac- To select the best one, we use cross-validation. specifically k-fold cross-valle validaaii. The sample S is first divided into A equal disjoint subsets S to Sz. Each algorithm A. is applied & times, giving us the k models A.(S), A.(S),))),..... St),.... Ac(S Sz). We can estimate the true risk of each of these models through its MSE on the holdout set, e.g. MSED (S) R MSFs, (A)(S),S)) In turn, we can estimate the truc risk of Ac(S) (the model trained on the entire dataset) as the average MSE of these models:
1 k MSEs, (Ac(S) (.)).
~ (A(
(5)
k
3

---

Published in Transactions on Machine Learning Research (01/2025)
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 el al. (1553) as uniquely interesting attribution method in the context of a game with n players. Denoting [n] == [I....,m] the set of natural numbers up to n, and P(nn) its power set, the Shapley value revolves around value function P([nl) + R. Given
coalition of players I [n], is the payoff that this coalition can achieve in the game by playing together. The Shapley value of player is defined as:
M1
17| - 1)! ( n!
(:(1)= = Σ
(6)
IC
It is the average value added by player i when it joins a coalition where 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:
ov(P) = ((nl))
(7)
In simpler words, the gain is distributed among all the agents/players.
2. Symmetry. Two players with equal contributions have equal Shapley values:
(VI ( [n]
(()) v(ru = 0.(1) = +,
(8)
3. Null player. A player that has no contribution has a Shapley value of 0:
- (v) =
(9)
4. Lincarity. The Shapley values of a sum Vi 2/2 of value functions are the sum of the Shapley values of the individual functions:
Vi E (: 12) = (11))) (122).
(10)
These properties make the Shapley value a uniquely interesting attribution method for the contribution of each variable to the prediction of model. Although some raise doubts on its reliability (Huang & Marques-
val 2023 or the desirability of these four properties (Fryer el al. 2021 we do not address this debate
this paper.
Different value functions yield different Shapley values, the choice of value function is crucial. In this paper, wish to attribute to each variable in dataset its contribution to the output of a model. Consider the instance space X' = R", so that the natural numbers up to n, which we denoted [n] == [I),...m), identify the variables that make up an instance x = (II...., In.) E X. For some model h X' + R, what is the contribution variable E [r] to the model output h(>)? In general, is not possible to see what a model would output given a coalition C [m] of variables, since most models are undefined with missing values. To sidestep this issuc, the solution most commonly used (Lundberg & Leo 2017) is to intervene on the values of on the variables outside of the coalition. To that end, we define the replacement function. For any feature subset P C [no], and x, z E X, the vector rp(2,x) is 2, except for the variables in P that are replaced by the values in r::
2j i P Pi i P.
=
(11)
4

---

Published in Transactions on Machine Learning Research (01/2025)
Finally, suppose that we have instance of interest and background S X consisting of a set of instances. Denote by US) the uniform distribution over S. The value function P([n)) + R that we focus on is:
vT) = E x~U(S)
[2)r)
(12)
In simple words, the output of h is averaged over the background, but the variables in 7 are fixed to their value in æ. This value function yields a so-called interventional Shapley value, after the act of replacing feature values with those from another instance, i.c. "intervening" on these features. This is the Shapley value calculated by SHAP (Lundberg & Lee 2017
Rather than Equation ( we use an equivalent form using permutations of [n]] (see e.g. Equation (3)
of Khorrami Chokami & Rabitti (2024)). A permutation is bijective function : [n] + [m]; a(t) is the
position of variable in the permutation. Let Ω be the set of all permutations of [rz], and (S2) be the
uniform distribution over Ω. Also define:
T. {i E [n] x()) < <()).
(13)
the set of variables preceding (but excluding i) in the order implied by the permutation x. Given all of these notations, the Shapley value of variable i, for example for the predictor h, is defined as:
o SHAP (h, == E
E
[h(m, ((,,.
(14)
(N) (R)
(Definition taken from Laberge et al. (0023) .)
Calculating Equation 14] directly is in Onn)), as it requires enumcrating every permutation of [n]. It therefore cannot reasonably be used except on particularly low-dimensional datasets. However, some model famílies admit much faster formulas or algorithms for calculating the Shapley values. Such model families in- clude linear models (Lundberg & Leel 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 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)
Σx Wyzj are obtained easily
from its coefficients by the formula osHAP (h,z), x) = w/ (2; EZ-U(S) 21) Lundberg & Lee 2017;: the Shapley values of GAM h(x) = E, (xx,) are given by:
SHAP (h, x)) = [((xs) - z~U(S) [((z))]
(15)
(The proof is in the appendix.)
Structured Additive Regression (STAR) models are a generalization of both linear models and GAMs (see c.g. Mayer et al. (201))) ). STAR model allows variable interactions by defining the predictor over subsets
of the data features:
h(x)
Σ f((21), Iglm]
(16)
where xy the partial vector tantly, not all of these functions fi need to be nonzero. We can indeed choose which variable interactions are desired by selecting which feature subsets I are associated with a nonzero fr. For example, GAMs use only singleton feature subsets, i.c. individual variables, while GA2Ms (Lou et al. 2013) up to pairwise interactions between the variables.
and fi is a function defined over the feature subset I. Impor-
STAR models have various names in the literature. Bordt & von Luxburg 2023 define Generalized Additive Models of order & as models of the form:
h(x)=
Σ [((x)), IEMR),ISS
(17)
and introduce k-Shapley Values, which attribute a contribution to each feature subset, instead of only in- dividual variables. This generalizes Shapley values, but suffers from the same high computational cost.
5

---

Published in Transactions on Machine Learning Research (01/2025)
Functional Decomposition, following previous work the sub- ject [Hooker] 2004; Kuo et al. 2010 and show among other things how to obtain the Shapley values of a
Laberge et al. (0244) call Equation [66
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 (hhaley et al. (1953). Equation 14 Equation (11) of [ordt & von Luxburg (2023); and algorithms like SHAP LLundberg & Lee (2017). These algorithms all suffer from the same delibitating complexity with regard to the total number of variables. On the other hand, model-specific formulas, such 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 algo- rithms have been developed. Notably among these are KerneISHAP (Lundborg & Loe. 2017]. Unbiasod KernelSHAP (Covert & Leej (2020) and SHAP-IQ (Fumagalli al.) 2023). These provide approxinative 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.c. 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 the maximal number l 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 KerneISHAP, and SHAP-IQ. We also study the quality of these approximations in relation to the exact values provided by our algorithm.
As final note to this section, we wish to mention the recent work of Khorrami Chokami & Rabitti (024 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 v(I)= := Var, /Eerr fi(z;)], which yields, as the corresponding Shapley value, a global variable importance index, opposed to 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 122 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, Φ. K) of a parameter space W; base predictor ő : W x X + R; kernel K. W W R; distribution p over W. The kernel K. implicitly defines a reproducing kernel Hilbert space (RKHS) H. space of W + R with particular properties. (A thorough understanding of reproducing kernel Hilbert spaces is not required for this paper, but the interested reader can see Dube & Marchand (2024) for a short primer on RKHS.) They define the operator A H + (X), which takes as input a weight function a E H. and returns the bounded predictor Aa : X' + R via the equation:
functions of the form
Aa(z)= = E [o(w))(w,x)). wmp
(18)
6

---

Published in Transactions on Machine Learning Research (01/2025)
In practice, the weight function H. takes the form of sum Σ a (uw,) (a property of the RKHS). for some real coefficients (1,..., dr and w1,...;; WT E W. The model can then be written as:
Aa(x) =
Σ aE whpp
[Kwe, ())
(19)
Therefore, an RKHS Weighting model is weighted sum of the partial functions defined by the expec- tations Eu-p [[,, w))(wxr)). Each such function of x is parameterized by a W, E W. While this is in theory highly flexible family of models, the current bottleneck for its utilization is the requirement for ex- actly calculating these expectations. This can be done, depending on the choice of instantiation (W, o. p). Dubé Marchand (2024) presented two instantiations of the model, which they called I1 and 12. We will refer to these RWSign and RWStumps respectively after their base predictor Ø, a simple sign function for RWSign, and a decision stump for RWStumps, with RW standing for RKHS Weighting. The reader can refer to Table 8 in Scction 4 for the explicit definition of RWSign and RWStumps, and to Table 1 of Dubé & Marchand 2024 for the analytical form of the expectation Emrep [Ki w)}(w, In both cases, the entire model [Equation 190 becomes a weighted sum of error functions a highly expressive form.
Since Equation 18 represents an aggregation of base predictor predictions (with regard to a distribution p and a weight function a), 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 new instantiation which we expect be better suited regression, using the rectified linear unit (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 cxisting 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 I
which we introduce in the next soction.
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 STAR model). Consider STAR model:
h(x)= = Σ f((1).
(20)
I5(m)
Consider the replucement function rp defined in Equation and the sels:
A+ 1/ A,1 A < μА}
(21)
A (1 //1 A, A < [/].
(22)
Then the Shapley values of the model are:
(s(0T) f((21) +
@SHAP (h,
E INU(S)
Σ
I:f130,AEr
M((A(Z,))1) (A((4)
Σ [((A)1((Z
Σ
(23)
(A((A))
AEA 6.1)
AEA (..)
Proof. The proof is a series of combinatorial arguments, and can be found in the appendix.
Perl(x) = /ズ 1 c dt ReLLU()) =)

---

Published in Transactions on Machine Learning Research (01/2025)
Assuming that the feature subsets I are limited in size, then Equation 23 is vastly more efficient calculate
than the definition of the Shapley values in Equation (14p. Instead of going through all the permutations of
[n], Equation (23) requires only enumerating the subsets of each feature subset 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 compulations. Indeed, observe that:
1. For any E A E (i, 7); we also have A E (J, 7). This means that the computation ff(ra(z,a)1)
appears in the calculation of the Shapley value of each variable in A.
Consider any A A (i, and variable j NA. Then AU( e Ai,),,, and so the computation (1(T(aut))))))) (2, is equal to f1(r
In short, when calculating the Shapley values of all variables (rather than a single one), the same computation [r(ra(z,z)1) is done exactly ||] times; |A] times from the first observation, and ||] - |A] times from the second observation. We can save factor /| by doing this computation only once. The term [r(z))=f((rr) can also be calculated once, as the value is the same for all variables. These optimizations in mind, we present Algorithm 1 which simultancously calculates all the Shapley values of a STAR model for given instance.
Algorithm STAR-SHAP: Shapley Values of a Structured Additive Regression Model
input Input space X' R", instance a E X, STAR model h(x) Eigin) f(()), sample S E xm Initialize &- = (#1,..., on) 0 E R for E S, C [n] such that fr + 0 do baseline +- f1(21)))((ar)
for i I do
di <- Φb + baseline
end for
for A C I, A # 0 do
term +- f=)()
for i E 7 do
if i E A and
AE A (i, then
6 + 0: +
term
) [
else if i 1 A and (AU E A-(I, then
Φ 0 (1
end if
end for
end for end for output m
A highly interesting characteristic of Algorithm[] 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 fr( x))). Note that, at this point, we make no assumption on the partial models fi other than that they utilize |/] variables. Let therefore suppose that some nondecreasing function F -N + [0, oo) characterizes the computational complexity of the partial models, so that the cost of calculating (z, is in O(F(())) for all I (and all z). Denote k the maximal size of the feature subsets associated with nonzero partial model (fr * 0), and N the number of such feature subsets. Algorithm [ is loop over these N feature subsets and over the m examples in the background S. For each feature subset 1, the algorithm then iterates over all its subsets C I, calculating the partial output fr x)) each time. Since there are 2 subsets to set of size k. this dominant operation is done a total of mN2 times, cach with cost of ((F(())) Since |/| < k. this leads to favorable algorithmic complexity of O(mN2* F(k)), independently of the dimensionality of the data. One could argue that N, the number of feature subsets defining the
8

---

Published in Transactions on Machine Learning Research (01/2025)
Algorithm
Target models
Fxact?
Complexity
TrPeSHAP
Trees, ensembles of trees
Yes
OmTTDD2)
STAR-SHAP
STAR models
Yes
(MN22F((k))
SHAP
Any
Yes
O F(n))
KernelSHAP
Any
No
O(Kn² + + +mF((n)
Unbiased KerneISHAP
Any
No
O(Kn² + KmF(nn)
SHAP-IQ
Any
No
O(Kn² + Km)(n))
Table 2: Computational complexity of various Shapley value algorithms for calculating all Shapley values of 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 tree. STAR-SHAP: N is the nunber 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 & 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 KerneISHAP 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 implemen- tations 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 (Funagalli el al. 2023
model in Equation 16 can grow exponentially in the dimension 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 | 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 see, SHAP has 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 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 (fr 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,
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, defined in Equation 23 includes Generalized Additive Models (GAMs) (Hastic) 2017), Generalized Additive Models with Pairwise Interactions (GA2Ms) (Lou et al. 2013), Sparse Polynomial Regression (Huang el 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), but for the others, Algorithm [ can do so as well, as long as the number of variable interactions moderate. Finally, another family of STAR models consists of STAR RKHS Weightings, which we introduce in the next section.
for any STAR model. This
9

---

Published in Transactions on Machine Learning Research (01/2025)
RWSign
RWRelu
RWStumps
W
R
I
(I,... nn x IR
2(w,))
sign( (w, 2)) wll
max(0, (w, at) (w,x)) 1 exp
sign (Twr w22
1[u1 = wil exp (w2 (u)-
x(u, w)
exp o(
24
24 ol)
212
u((1,...,n)) n2) x 0()
Table 3: Three RKHS Weightings instantiations.
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 (022)) 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 (W, o, P) of parameter space W, base predictor o, kernel K. (implying an RKHS H) and distribution p. For example, instantiation I1 Dubé & Marchand (2024) (which we call RWSign) uses W X' R", o(w,r) x) sign( (w, z)), the gaussian kernel K(uu, w) exp | and a normal distribution p = _N ((0,
The instantiation that we propose, which we name RWRelu, consists simply in replacing the base predictor of RWSign by the rectified linear unit (ReLU), ie. æ(w,z) ReLU (W, x)) max(0, (w, (See Table 3 for quick overview of each instantiation.) In order to calculate the output of the model, which takes the form:
Aa(z)=
E [K(w, w))(w,))]
(24)
for some a = Et- a.K(w,,') 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 K.p) where W X R", o(w,z)
max(0, x)), K(u, w) exp (- (- ||===|6 o2,), Then we have:
E urpp
-n/2 11
|
(w/.a)2
(u(,z)
(/2(111)) )1) (25)
(ul,
(1+ 9²)
2 12.²
V2çe
2-1 + V
+ erf
2/T
|
and v' = (1 21)-1a
1 where Ç is defined by the relationship 2(2
1 +
2x
202
Proof. The calculus can be found in the appendix.
Dubé & Marchand (2024) also define some instantiation-dependee theoretical constants (see Theorem 4 of Dubé & Marchand (2024)) which are then used in the theoretical guarantees. We provide those values for RWRchu and the relevant calculus in the appendix. Importantly, this means that all thcoretical guarantoes 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 o, K. We will show how to instantiate RKHS Weighting which is also STAR model, first. in full generality, and then by
10

---

Published in Transactions on Machine Learning Research (01/2025)
transforming already existing instantiation into STAR model. By having STAR model which is also RKHS Weighting, we keep all of the benefits of the RKHS Weightings, i.e. its theoretical guarantees, learning algorithms,, pruning methods (Dubé & Marchand 2024p, but also gain interpretability through the computation of the Shapley values, using Algorithm !1
The defining characteristic of STAR models that they are built. with functions that use only subset of the variables time. The key defining an RKHS Weighting which is also a STAR model is therefore the base predictor Φ. By choosing base predictor which can flexibly change the features that it considers, can obtain STAR model. We now go through each ingredient and see how they must be defined.
The parameter space W is accompanied by the space of feature subsets P([x)) (the family of all subsets of [n]).
2. The base predictor takes as input both parameter w and a feature subset L, as well as an instance z; (((w,, x). == d1 (w,,1) could be any predictor defined on the feature subset I.
3. The kernel must take into account the feature subsets, mcaning it is now function of the form K. (W x P((n)) x (W x P((n))) I.
The distribution p on x P([n]) generates both parameters and feature subsets.
All in all, the instantiation takes the form (W x P([n)), o, K.K. p). Now, to show that this is indeed STAR model. Given weight function Ei= axK(((u, [t),), the output of the predictor arbitrary instance is given by:
Aa(x)=
E (u.J)<p
I)((w, [),
F. [K((we, [c), (w, ()),,
at
(w.1)ppp
Σ
Ot
E (w,J)<p
[K((,, Te), (ve, (w,,
(26)
JCIn
We have thus expressed the predictor as a sum of functions defined on feature subsets. This is therefore STAR model, with the put on an arbitrary x given by Aa((x) = EEn( (a(x),, with:
fz(ms)
Σ at E U.Lj (w,J)pp
[K(Wu, It), (w, (or
(27)
This framework for expressing an RKHS Weighting as a STAR model is in full generality. Of course, prior knowledge specific given problem can be inserted into the choices for base predictor, kernel distribution. For instance, might be interested only in sequences of features when presented with text, time series, other sequential data, in which case the distribution can simply put a weight of 0 to any non-soquential feature subset. This family of STAR models is therefore highly flexible.
4.3 Transforming an existing instantiation into STAR model
Perhaps the simplest way to instantiate an RKHS Weighting STAR model using the framework described in the previous iontion is to take an existing instantiation (W, d, Kw,Pw) which satisfies few simple conditions:
1. X W R
The base predictor pfw, .)) is a function of (w, x), ie. 6(w, = ((ww+)))
3. The distribution pw W generates each component independently.
11

---

Published in Transactions on Machine Learning Research (01/2025)
Existing instantiation
STAR version
Parameter space
W
W x P([n]) s(ww, I), x) (((w1,)))
Base predictor
Ó(w, x) - di ((),
Kernel Distribution
Kw
K(((w1, T1), (w2, 12)) = (w((,, 12]
Pw
p(w, 1) = Pwn (wpr(w))ppp (M)
Table 4: Cheat sheet for transforming an existing RKHS Weighting instantiation into a Structured Additive Regression model.
In that situation, modifying the base predictor to work on feature subsets is trivially easy:
o(ww, 1),) =((1x)))))))
(28)
In other words, simply ignore all the variables not in the feature subset. As for the distribution p. the independence allows us marginalize away the unwanted variables. The distribution p(w, I) is the product (w)) (T), where P(n] is chosen distribution over the feature subsets, and Pwiz (w) is the original distri-
bution Pw marginalized to the features in I. Finally, we elect to use the indicator kernel Knn) (4,, 12) = 1
I2] over the feature subsets, and use the kernel product K((W,, It,, (w2, I2)) == (w W2 [ma] (T, I2) I We
end up with the new instantiation (W x P([n)), o, Kw x K(n), p). (Table summarizes all of this.) Assuming weight function -ET a,K((u,, I),,), the model is:
Aa(x) Σ Σ JCIN) #I=U
apPlo)(T)) (Te)
E wo-pille
(29)
The only element that still needs to be further specified the distribution Pln] 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.
The simplest option is to set Pln as the uniform distribution on feature subsets of size equal to k where & becomes a hyperparamoter of the algorithm. A slight variation of this is to use feature subsets of size up to k allowing a larger variely of variable interactions.
2. Some datasets, like time series or genomic data, have structure, where adjacent features are corre- lated in some meaningful way. When presented with such sequential data, it is natural to consider sequences of features. We can then set Pin] as the uniform distribution over sequences of length k (or length up k).
Because of the presence of the factor (1e) in Equation (29). using the first option (arbitrary subsets of size k) with large is not recommended. Indeed, since PIx (7. 1/()) R _ the model output tends to zero exponentially as increases. Large coefficients ag are then required to offset the small density. This leads to a large norm ||α|| of the weight function, and therefore poorer guarantees (Dubé & Marchand 2024). More intuitively it is important that significant portion of the feature subsets I, 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 A in an n-dimensional dataset, Pin] will always be of the order 1 .. When considering sequences of length up to k, then P(n] will be about Fm which is undeniably preferable over 1 *
WW the indicator kernel for two related but distinct reasons. The first is the assumption that, in general, the comparissn between parameters, i.c. the keruel evaluation Kw(w1,32) only makes sense if wr and wy are parameters over the same features. that the kernel evaluation Kw(wi, wx) well-defined. There
Sccond, ensures that wi and we are the same dimensionality, some abuse of notation here, as we the same symbol Kw to denote function over potentially different spaces (e.g. wi could be from R3 and w2 from R4), ), but deem it benign enough, as does simplify the expressions.
12

---

Published in Transactions on Machine Learning Research (01/2025)
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é & Mare- hand (2024) for STAR RKHS model. Once again, this means that the thcoretical 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 STAR model; RWRelu, a new RKHS Weighting instantiation; and a new family of STAR models, obtained by transforming cxisting RKHS Weighting instantiations into STAR models. We present in the next few sections three experiments to test these new lools. First, we compare the computation time of STAR- SHAP other Shapley value algorithms in Section 5.2 Then, Section 5.3 presents 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 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 datasct.
5.1 Generating synthetic datasets
To generate our synthetic datasets, we took inspiration from the sum of unanimity models used by by/umagalli
et al. (022)) For given values of k (the size of the feature subsets), n (the total number of variables), and T (the desired number termss), we generate a k-STAR synthetic model over boolean instances (X' (=1
as f(x) Din [Vj E L1 = 1], where 4I..., /T are uniformly and independently generated feature subsets of size k, and a1....,00 dr are coefficients uniformly generated in [0, 1]. To go along with every random model, we also generate random dataset of boolean-valued instances S x where each component of an instance x is independently and uniformly sampled from {+1}.
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, 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 [b clearly illustrates the fact that STAR-SHAP lacks an explicit dependency on the dimensionality of the data, making much more efficient than other algorithms at higher dimensionalities.
also show the computation time growth of using STAR-SHAP
However, Figure [ does not account for the quality of the approximations, and only looks 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-1Q coincide when calculating 1-Shapley values (Fumagalli et al. 2023), possess thcoretical guarantecs, and have better empirical performance than KernelSHAP (according to Fumagalli et al. (2023)), we assess the quality of SHAP-IQ in Figure? [ The figure shows how the advantage of STAR-SHAP increases with dimensionality. Indood, not only does the computation time of SHAP-IQ incrcase, but also the quality of its approximations degrades, requiring additional computation time return good approximations. Meanwhile, the values returned by STAR-SHAP are always exact, and are calculated efficiently.
13

---

Published in Transactions on Machine Learning Research (01/2025)
Kernel SHAP
Kond SHAP
SHAP
12
SHAP4Q
SHAP44
Uniilosed Kamei SHAP
STAR-SHAP (n-STAR)
STAR-SHAP (5-STAR)
Unalased Kemel SHAP
10
STARL-SHAP (1-STAR))
Time (0)
Time (4)
10 12 16 18
250
420
500
803
1000
Dimensionelity
Dimens omennioonallyyoo
(a)
(b)
Figure 1: Computation time of the Shapley values of 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 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 5-STAR model cannot exist. Figure 1b is the rescaled continuation Figure Ta without the exponentially increasing algorithmss 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.) (024) of KerneISHAP, Unbiasod KernelSHAP and SHAP-IQ itself, with 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).
14

---

Published in Transactions on Machine Learning Research (01/2025)
T-Coefficiett determination
1.2
2-SHAP-IQ
1 -oeffiiien of determination
1.2
3-SHAP-10
3-SHAP-IQ
4-SHAP-IQ
2-STAR-SHA
3-STAR-SHAP
4-SHAP-IQ
4.STAR-SHA
1.0
3-STAR SHAA
4.STAR ESTAR-SHA--
SHAR..R.....
SHAP
G-STAR-SHAF
5-SHAP-IQ
6-SHAP-IQ
7-SHAP-IQ
0.8
7.STAR
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
0.2
0.6
Time (s)
Time (s)
(a)
(b)
- C--Cefiiieenn determination
6-SHAP-IQ
1.1
7-SHAP-IQ
8-SHAP
1.0
9-SHAP-IQ
10-SHAP-IQ
0.9
&-STAR-SHAP..
B.STAR.SHAP..
10-STARKSHAA
0.8
S-STAR- 7 SSTARRSHHA
0.7
25
50
75
100
125
Time (s)
(c)
Figure One minus the coefficient of dotermination (R²) of the Shapley values of synthetic models ap- proximated 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 the models (T = 100), and varying dimensionality For cach value of k. the size of the fcature subsets in the STAR model, the Shapley values of synthetic model on dataset of size 10 were calculated approximately by SHAP-TQ (A-SHAP-TQ), and exactly by STAR-SHAP (K-STAR-SHAP). The R² was then calculated between the vectorized 10 x n matrices of approximated and exact Shapley values. A vertical dotted line labeled A-STAR-SHAP indicates the time STAR-SHAP took to calculate the exact Shapley values for that value of Every point the average of 10 independently seeded runs. Figure 2a on low-dimensionality datascts (n = 10 variables,
100 terms in the model), SHAP-IQ starts generating Shapley values with k
This means that STAR-SHAP is the better algorithm choice for & < 3, and then gradually losses its merits 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 competitive perfectly exact. Shapley values are not required. Finally, Figure 2 used n = 500 total vari- ables (with T = 100 terms) to build the models. Even for = 10, which an unusually high number of variable interactions for model, SHAP-IQ generates barcly acceptable Shapley values when given the same amount of time STAR-SHAP. For any smaller value of to and the same computation time STAR-SHAP, SHAP-1Q yielded extremely poor quality approximations of the Shapley values.
R2 coelficient larger than 0 faster than STAR-SHAP
to about = 5, after which SHAP-IQ starts becoming
15

---

Published in Transactions on Machine Learning Research (01/2025)
Training size
Test size
Dimensionality
Source (clickable)
abalone
3132
1045
10
UCT
diabetes
331
111
10
Scikit- learn
housing
15480
5160
8
Scikit earn
concrete
772
258
8
UCI
conductivity
15947
5316
81
UCI
wine
133
45
13
UCI
Table 5: Datasets of the regression experiment.
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 limitod, then the computation time of STAR-SHAP will be fast. Since 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 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 E 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 hyperpa- rameters, we trained a varicty of models on several regression datascts taken from the UCI Machine Learning Repository (Dua & Grall 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 to answer a few important questions.
Q1. How does RWRelu, the ReLU version of an RKHS Weighting, compare 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 cases.
Q2. How do the STAR versions of RKHS Weighting models perform in comparison to their original coun- terparts?
R2. Thore is a modorate loss of performance going from RWRclu to its 1-STAR version, a cost that might be acceptable to pay for the immense gain 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 Conductivity is caused by the higher dimensionality of the dataset. (See this footnote 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.
We recall the presence of factor P in Equation 29 the expression for STAR RKHS Weighting. When using the uniform distribution over all possible feature subscts of size k we have (1e) = 1/( x ne, which is extremely small value high-ddmmessonallity datasets, leading vanishing outputs of the model. And indeed, the offending dataset, Conductivity, far higher-dimensional than the other ones, proving the point.
16

---

Published in Transactions on Machine Learning Research (01/2025)
Table Comparison the regression performance STAR RKHS Weightings and the ReLU instantiation, learned using Algorithm 3 of Dubé Marchand (2024) (the Least Squares fit) with T = 5000, to 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.
Dataset
abalone
diabetes
housing
concrete
conductivity
wine
Algorithu
Decision RegeRsgressor
0.499
0.105
0.677
0.796
0.882
0.867
EBM
0.552
0.331
0.822
0.931
0.903
0.897
KernelRidge
0.587
0.282
0.771
0.831
0.872
0.951
LigeesReg
0.543
0.359
0.591
0.623
0.736
0.804
RWRelu
0.577
0.349
0.734
0.835
0.842
0.832
RWRelu 1-STAR
0.555
0.314
0.662
0.780
0.777
0.808
RWRelu 2-STAR
0.577
0.305
0.731
0.847
0.805
0.843
RWRelu 3-STAR
0.580
0.332
0.744
0.847
0.568
0.845
RWRelu 4-STAR
0.579
0.350
0.737
0.846
0.032
0.835
RWRelu 5-STAR
0.578
0.334
0.729
0.844
0.001
0.837
RWRelu [1, 2, 3, 4, 5]-STAR
0.577
0.298
0.732
0.808
0.779
0.826
RWSign
0.578
0.348
0.711
0.828
0.728
0.893
RWSign 1-STAR
0.350
0.302
0.354
0.471
0.685
0.750
RWSign 2-STAR
0.579
0.369
0.696
0.835
0.816
0.868
RWSign 3-STAR
0.586
0.339
0.720
0.840
0.468
0.912
RWSign 4-STAR
0.585
0.349
0.708
0.794
0.003
0.896
RWSign 5-STAR
0.580
0.333
0.705
0.770
-0.000
0.899
RWSign [1, 2, 3, 4, 5]-STAR
0.580
0.399
0.698
0.826
0.694
0.840
RWStumps
0.567
0.339
0.768
0.911
0.859
0.864
SVR
0.564
0.293
0.770
0.799
0.862
0.950
17

---

Published in Transactions on Machine Learning Research (01/2025)
Training size
Test size
Dimensionality
ChlorineConccentration
6990
57585
10
Computers
16434
16434
10
RCG5000
6487
58487
10
FCCesUCR
2388
24588
10
LargeKitchennApplianess
24684
24684
10
MellourrnePedestrriann
3579
7314
Table 7: Datascts of the time series prediction experiment.
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 of the 6, and of the 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 least sometimes competitive.
To conclude this experiment,
that transforming RKHS Weighting into STAR model is valid way of learning an interpretable model, at relatively small cost of performance. However, other methods, especially Explainable Boosting Machines, often outperform these RKHS Weightings in both accuracy and interpretabilily.
5.4 Time series prediction performance
Perhaps the biggest takeaway from Table D is the excellence of the Explainable Boosting Machine 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 fcature 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)). This makes the STAR RKHS instantiation much better suited to solve regression problems for which adjacent data features are correlated, such time series. Therefore, 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 has length 10, and the value to predict y is the 11-th value. See Table 7 Hyperparameters were selected by 5-fold cross-validee voos-vaaiiation.. See the appendix for the detailed experimental setup.
The results of this experiment can be found in Table The RKHS Weightings fare much better in this experiment than the previous. Here, they (slightly) outperform the EBM on 5 of the datasets. On the other hand, the EBM does systematically reach the highest training R², 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 smaller generalization gap. Finally, we notice the expected trend that larger feature subsets the STAR RKHS Weightings lead to better results, with the full RWRelu instantiation
18

---

Published in Transactions on Machine Learning Research (01/2025)
Table 8: Comparison of the time scries prediction performance of STAR RKHS Weightings and the ReLU instantiation, learned using Algorithm 3 of Dubé & Marchand (2024) (the Least Squares fit) with 7 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 cach datasct have been bolded.
Training H²
Test R²²
Training time (s)
Dataset
Algorithm
ChlorineConeentration
EBM
0.955
0.916
63.853
RWRelu
0.934
0.928
3.307
RWRelu 1-STAR
0.638
0.641
36.375
RWRelu 2-STAR
0.813
0.805
40.390
RWRelu 3-STAR
0.907
0.899
44.500
RWRelu 4-STAR
0.928
0.920
48.720
RWRelu 5-STAR
0.939
0.927
54.159
RWRelu [1, 2, 3, 4, 5J-STAR
0.924
0.906
19.823
Computers
EBM
0.882
0.842
16.291
RWRelu
0.821
0.851
6.085
RWRelu 1-STAR
0.786
0.838
44.742
RWRelu 2-STAR
0.802
0.840
50.501
RWRelu 3-STAR
0.811
0.848
55.049
RWRelu 4-STAR
0.814
0.847
58.744
RWRelu 5-STAR
0.818
0.848
64.858
RWRelu [1, 2, 3, 4, 5]-STAR
0.817
0.849
29.086
ECG5000
EBM
0.943
0.886
23.525
RWReiu
0.900
0.882
3.144
RWRelu 1-STAR
0.843
0.841
36.361
RWRelu 2-STAR
0.877
0.871
40.024
RWRelu 3-STAR
0.883
0.875
43.769
RWRelu 4-STAR
0.893
0.879
48.122
RWRelu 5-STAR
0.899
0.882
53.384
RWRelu 21 3, 5]-STAR
0.894
0.882
19.157
FCcesUCR
EBM
0.868
0.661
7.674
RWRelu
0.715
0.704
1.929
RWReiu 1-STAR
0.679
0.698
33.011
RWRelu 2-STAR
0.690
0.702
35.613
RWReiu 3-STAR
0.699
0.705
38.827
RWRelu 4-STAR
0.700
0.705
42.990
RWRelu 5-STAR
0.709
0.706
49.571
RWRelu [1, 2, 3, 5T-R
0.704
0.706
15.570
LargeKitccen
EBM
0.840
0.743
19.103
RWRelu
0.813
0.743
8.775
RWRelu 1-STAR
0.764
0.724
52.104
RWRelu 2-STAR
0.786
0.748
56.938
RWRelu 3-STAR
0.777
0.736
60.735
RWReiu 4-STAA
0.803
0.746
68.894
RWRelu 5-STAR
0.810
0.745
75.208
RWRelu [1, 2, 3, 4, 5J-STAR
0.812
0.749
37.741
MelbournePedestrian
EBM
0.979
0.924
9.731
RWRelu
0.946
0.925
2.335
RWRelu 1-STAR
0.888
0.888
35.930
RWReiu 2-STAR
0.921
0.902
36.519
RWRelu 3-STAR
0.943
0.923
40.523
RWRelu 4-STAR
0.943
0.925
44.806
RWRelu 5-STAR
0.943
0.924
50.667
RWRelu [1, 2, 3, 4, 5J-STAR
0.941
0.919
16.232
19

---

Published in Transactions on Machine Learning Research (01/2025)
often being the best (although by 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 defnnd, instantiations would perform much better than RWRelu in this experiment. Indeed, it still an open question how to craft instantiations that can perform particularly well on 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 more traditional model, and how its predictions could be interpreted using Algorithm [ (STAR-SHAP). Figure 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 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.
High
High
Latitude
Latitude
Longitude
Longitude
Medinc
Medinc
Feature value
Feature value
AveOccup
AveOccup
AveRooms
HouseAge
HouseAge
AveRooms
AveBedrms
AveBedrms
Population
Population
Low
Low
-2
2
-2
2
SHAP value (impact on model output)
SHAP value (impact on model output)
(b)
Figure 3: Beeswarm plots of the Shapley values on California Housing of the Explainable Boosting Machine (Figure 3al 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 3al 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.
Figure g compares the Shapley values of both models on single instance. While the predicted values are notibeybly 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 secing the RKHS Weighting plot as the EBM one. This gives further credence to RKHS Weightings as promising model family, especially interpretable when cast STAR models.
20

---

Published in Transactions on Machine Learning Research (01/2025)
rix) 0.778
fx1 - -.255
Longitude
+1.26
-1.562 Longitude
15.25
0.007 AveOccup
0.007 AveOccup
1.808 = Medino
0.085 AveBedrms
0.998 AveRooms
0.437 Population
+0.06
0.085 AveBedmm
+0.01
0.764 HouseAge
00.05
0.764 HouseAge
10.01
0.988 AveRooms
-0.04
0.437 Population
+0
1.808 Medinc
+0.03
-0.5
0.0
0.5
1.0
1.5
-1.50
-1.25
-1.00
-0.75
-0.50
-0.25
0.00
0.25
E((X)
EU/20]
(a)
(b)
Figure 4: Waterfall plots of the Shapley values of a single instance r from Figure 3-STAR RWRelu. The true f(x) value for that instance is - 0..6066.
Figure 4a
Explainable
Boosting Machine. Figure 4b
5.6 Discussion
A number of conclusions can be reached from experiments.
Figure M demonstrates that STAR-SHAP is the preferred algorithm for calculating the Shapley values of a modol with moderate variable interactions, especially so the more total variables are present the data. Figure more clearly shows how the advantages STAR-SHAP increase with the total number variables in the model, as the quality of the approximations returned by approximation algorithms degrades with dimensionality.
Table G shows that RKHS Weightings and STAR RKHS Weightings reach moderately good performance arbitrary datasets. They have 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 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 step forward in understanding how to use RKHS Weightings, and proof that they be high accuracy, theoretically grounded (as per Dubé & Marchand (2024) ) and even interpretable models. Further rescarch to expand the usability of RKHS Weightings, as well as find their unique strengths and applications, warranted.
6 Limitations and future work
Our work presents number of limitations, some of which can be addressod with future rescarch.
Algorithm complexity. In thcory, Algorithm can calculate the Shapley values any STAR model, but the computational cost 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 [ and 2 show that, for moderate interaction sizes, is the preferred algorithm to use 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 complex integral for any new instantiation. We have introduced new instantiation, but it
21

---

Published in Transactions on Machine Learning Research (01/2025)
is still 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 g 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 G the most generic version of STAR RKHS Weightings, i.e. using the uniform distribution on all possible feature subsets of size have somewhat limited use, since their output scales as 1/(x) (n being the dimensionality of the dataset). As we have already addressed, the best option is to insert prior knowledge into the model,
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
Tn 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 higl-dimeensioal datasets. We introduced 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 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 generic regression datasots, 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 lypes 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 Engi- neering 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 Bombardicr inc G
Special thanks to the anonymous TMLR reviewers whose comments and suggestions led significant provements of the manuscript.
hhttps://deel.gue quebec
22

---

Published in Transactions on Machine Learning Research (01/2025)
References
Alain Berlinet and Christine Thomas-Agnan Reproducing kernel Hilbert spaces in probobility 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 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 21th 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 preprinl arXiv:2012.0155 2.11536, 2020.
Hoang Anh Dau, Eamonn Keogh, Kaveh Kamgar, Chin-Chia Michael Yeh, Yau Zhu, Shaghayegh Gharghabi, Chotirat Ann Ratanamahatana, Yanping, Bing Hu, Nurjahan Begum, Anthony Bagnall, Abdullah Mueen, Gustavo Batista, and Hexagon-ML. The uer time series classification archive, October 2018. https:
7/www.co cs ucr wdweeduuulaeee
_atrriie_dtt____ee 2ateeddtaa_a1aa////aaa
Dheeru Dua and Casey Graff. Uci machine learning repository, 2017. URL http::/arrcive. uci ucw.edu/
Gabriel Dubé and Mario Marchand. RKHS Weightings of Functions. working paper preprint, July 2024.
URL
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 http:///pppi.smmmmmccccoollll org/CorpusII:
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://arxx:: www.aaiwww.....aaaaa.....eaoo ab2/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 exccution time of computer programs using sparse polynomial regression. Advances 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-thcoretic variable importance index for generalized additive models. Journal of Computational and Graphical Statistics, pp. 10, 2024.
Sotiris B Kotsiantis. Decision trees: a recent overview. Artificial Intelligence Review, 39:261 283, 2013.
F Kuo, Sloan, Grzegorz Wasilkowski, and Henryk Wožniakowski. On decompositions of multivariate functions. Mathemalics of compulation, 79(270):953 966, 2010.
23

---

Published in Transactions on Machine Learning Research (01/2025)
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
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.
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 dala 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.08880 2018
Michael Mayer, Steven C Bourassa, Martin Hoesli, and Donato Scognamiglio. Structured additive regression and tree boosting. Swiss Finance Instibute Research Paper, (21-83), 2021.
Grégoire Montavon, Wojeiech 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üillemmcer. shapiq: Shapley intcractions for machine learning. In The Tighty- Conference on Newral Information Processing Syslems Dulasels and Benchmarks Track, 2024. URL //pperev/operrrviee. net/
Harsha Nori, Samuel Jenkins, Paul Koch, and Rich Caruana. Interpretml: A unified framework for machine learning interpretability, 2019. URL htttp::/aaxit: wwwpii/.wrrg/bb//99.e... 09223
F. Podregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prottenhofer, 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. 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 Goudc, and Simon Shaw. Goneralizod additive models for large data sets Journal of the Royal Statisticul Sociely 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, orXiv preprint ar 202::2202000588,, 2024.
24

---

Published in Transactions on Machine Learning Research (01/2025)
A Proofs
First, we quickly prove Equation 15
Proof of Equation
15
We have:
6
SHAP (h, x) = zS)(S) T()
[ - -(Tr, x)]]
Σ fz(z,) E b((es) - -(I
f((tz)-
Σ fs(og)
xS)(() t))
jETU
Σ fz(Øz) - Σ
Σ 1s(2))
Σ 5s(z5)
+
2~((S) t0)
(EET;U(i)
jEau
(
jE
]
E
(fi(zs)
Σ f5(z)
E
f((,) + f((zt)
=S)(S) k~)(()
JEF.Ufi)
jEm.U(i)
=f((x)) x~((S) (1)
E
E
[f((z)]
=f((x))
E [)(()))]-
z~((8)
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 STAR model h(x) Eigin fr(m)), which
can be written as:
SHAP xH)) = E
E 2~U(S) T-((
[M (( (Fx,(2,)))
E xSU(S) (0)
E'
Σ rc(m]
fi(r (2
Σ f((rx((2,2)1)
15(m)
E (S)(S)
E --((
((( - fI( (2,2))]
[C[m]
We can immediately cut out most of the computations by noticing two facts:
1. Most functions fI are O.
2. If i I, then (2),)) Fm, (2,x),, Hence, the entire expectation:
F (Tr fr(r=,(z,x:1))
(30)
is 0.
The Shapley values therefore boil down to the following formula:
SHAP (h,
E zS)(s)
Σ
E T~)(a)
E (f((r, =0(())
(31)
IgIn):fi7o.ie/
Lemma
A.1 gives the value of Ex-(()) (2,
and Lemma A.2 gives the value of
E. Fx
Ti(rx,
25

---

Published in Transactions on Machine Learning Research (01/2025)
Lemma A.1. Consider the definitions of Section 2
assuming that i I, and the following:
(AC
(32)
In other words, A (i,I) is the family of strict subsets of I which contain i. Then we have:
fT(1T) |
Σ Jr(ra(z,z)1)
F
(1
(33)
AEA (5.1) (A((2)
Lemma A.2. Consider the definitions of Section assurnnng that i I, as well as:
A ( (AC f/ A,1 [ s [1}.
(34)
In other words, A- 1i is the family of subsets of which contain and at least one other element. Then
we have:
(2,))) = f1(21)
E x~U(0)
+
Σ
(35)
1
(A((AA)
AEA (1.)
Proof of Lemma A.1 We can separate the expectation into two parts:
E x~)(a)
Σ
fi ((
1
Σ
(2
(36)
n!
n!
JER IgRAUf()
TEN (iFU(
Term
Term 2
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 (a) (2, = æ, for all the permutations æ that satisfy is constant over all those permutations. is simply
I I Ta U {i}, and so the value (E, (1))),))((2,)))) fzz)). All that we need to do is calculate the number of such permutations. There is elegant combinatorics argument that gives the solution immediately.
Consider any permutation of [n]. This permutation satisfies the condition / m. (ii if and only if i is to the right of every other variable of I. This rightmost position is one possibility of out of [Z]. This means that exactly 1 out of every |4| permutations satisfies the condition. The number of permutations we seek is n! We therefore have: jr]
f1(x1) = f1(27) Term 1 1 n! [7
(37)
[7]
Calculating Term 2. Here, we consider all the permutations such that (T. U {i}) * 0, ie. the feature subset. contains variables not in n. U {i}. These variables will be replaced when calculating the function. Defining A = In ( U (i)), the set shared variables (note that we always have E A), we have:
Ts (2, =A(,=
The function (1 (2))) therefore takes a different value for each possible A, of which there are:
[ (1-1)
L (1,
k
14111 1) ) (1 -1
k
= 2 1.
(38)
26

---

Published in Transactions on Machine Learning Research (01/2025)
(We choose 0 up to [4] - variables to accompany in making up A (t), For given set A E (+, I) of shared variables, must calculate the number of permutations that are such that A = n (Ta U {i)). Call that number (n, [A)). We have:
N(n, 14144 JT(rA(2,a)1).
Term 2 = Σ
(39)
n!
(EA+(..)
Let us now calculate (n,, [7]]
1. First, choose the [7] positions for the variables in J. There are (7)) possibilities.
Of these !|| positions, there is only one choice for itself, since the variables in A (i) must be to the left of it, and the variables in I A to the right.
We can permute the |A| variables of A \ {+), the |r] |A| variables of \ A, and the n | variables in [n] 1, for a total of (IA] 1)((T] - |A]) 7|) permutations.
Assembling these facts gives us the formula:
(n, 14, Al) = (IA) - 1)((41 (A)) 14D!
72 (1)
(40)
We can simplify this further by expanding the binomial coefficient:
n (n, [A) = (IA) 1)(((() - (n 14))
(1)
n((IA)
1)1(11-- -AI1(1 12D)
1///(1 !/D)
n!! (AIGAD
(41)
This gives us the result.
Proof of Lemma
A.2 The proof is quite similar to that of Lemma
A.1
We can separate the expectation
into two parts:
E T~)
(2,2)))
1 n!
Σ
f1(r,,(2,x)1))+
Σ
(2,2),
(42)
m!
IOR
Torm
Term 2
Calculating Term 1. Notice that Tn. (2, x)2 = 21 for all the permutations स that satisfy In n. = 0, and so the value (2, is constant over all those permutations. is simply fi (2,,,)) fzr). All that we need to do calculate the number of such permutations. However, simple symmotry argument gives the answer. In the proof of Lemma A.1 we showed that the number of permutations such that I IS [i) is a! . In fact, we can notice that this mumber should be the same the mumber of permutations that satisfy mn n. = 0. Indeed, in the first case, all variables of / (i) must to the left of in m.. In the second, they must be to the right. There are exactly as many permutations that satisfy each condition.
Calculating Term 2. Here, consider all the permutations such that In,, # 0, i.c. the feature subset I shares at least one variable with ny. In this situation, some or all variables are replaced when calculating the function. Using the same definition A = I U (i)) as in the proof of Lemma A.1 we have:
TR (2,4)r (2,2))
27

---

Published in Transactions on Machine Learning Research (01/2025)
The function (r (2) therefore takes different value for each possible A. Since each A must contain at least two elements (i, and one more variable so that In πı * 0), and also A = I is now acceptable, the intersection A is taken from the set (ii which has the same cardinality as the set (+, namely permutations a such 2 We also know from the previous proof that there are Nn,, (1,,|11,14) = (aI((A)
that A n (TT {i)). We therefore have:
(16 (4)
Term
Σ
(43)
AEA (6.1)
which concludes the proof.
B
Calculus for RWRelu expectation
Here calculate the expectation Ewop [K(u, w))(w,x) 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 (024)..
Lemma B.1 (Lemma 15 of Dubé & Marchand 2024
Consider Hilbert space W. Let u, E W and
a, > 0. Then:
|| ull²
||ll
(' + !)
1
+ 1 + b |ull
a b
1
Lemma B.2 (Lemma 16 of Dubé & Marchand (2024) ). We have:
aallw-u al dw
π
Lemma B.3. We have:
[ ewo max(0, wr)dw >/0
(22y (oyye-22222 12/2- + Virul sign(x) + erf
/2
/27
Proof. We have:
00
(w-u)//22 / max(0, wx)dw I
1 -2 max (o, ( (v2nt u))x) 22
(t (t:=
w-
dt dw
'27
/27
v2y I e -t max (0, (V221 + )
Case 1. If x 2 0, we have:
V2y -t mmx (0, (v2yt + u)
22
000
--2
(V2y1 + w w
-4/V27
= -2
22
x
l0 1e dl +u
1
- dl
2/
1=00
0
c - dt + 1ee -
V2yx
V2y
+
2
-//V22
t--\/V2
-2yx
127²
+
VT
orf
u
+ V
arf(30) 1)
2
V2y
2
V2ye
12/2?
+ Viu
1 erf
(2y ))
u _m)
V2
28

---

Published in Transactions on Machine Learning Research (01/2025)
Case 2. If x < we have instead:
227 00
e- max (0., 0. ( (v2yt + u w)))dtt
= V2 [
w/v2/
6-1 0 (V27t + u) xdt
4/V22
- u/v2y
-2
V27
- dt + te
e
e- dt
t= 2
uu/V/22
0
f _ c dt
0
V2yz
V2
+
c-
dt
-00
1=-00
22
V2
1 1272
crf(20)
v
2
erf
(2
)
7
22ye
/27
+ Viu 1 erf
11
V2
/24
Both cases can be combined using the sign of æ, leading to the desired result.
Lemma B.4. We have
=1l1/113//2 max(0, (w, )d)
Ba
(V2m²² n- [
fu
(u,
(u,
))}
V2ye
-
V
+ erf
2
|
<271|<0
Proof. Calculate the integral using an orthonormal basis (vi,..., un of IR" such that on T
Write
w (wi,,... Wx) this new basis (i.c. w; = (w, 04) for all i), and similarly (uI,..., Un) for u. Under this change of coordinates, the integral becomes:
/ J wl]lw=ul12/00 max(0, (w, x)) )dw
-(wn Ene (w (m)- (wi (wn--a)) un) / max(0, / [///=////dw du..
R
We are left with product of n independent integrals:
{ e 1//22 max(0, (w. J:)
Ein' (w: (m:-
dwi dwn-1 / -w.. max(0, wn |
n
(w) u.) n.))//20 /252 dwi dwn-1 ( -wmn-un)) max wn ww]l=//)ddw
E
Rn-1
i=1
n
e (w, -u; 127 dws / (wn (wn max Wx 1 |
For cach i, Lemma B.2 gives us:
( (w:: / dw; 2
29

---

Published in Transactions on Machine Learning Research (01/2025)
Also, Lemma B.3
gives us:
(wo unn) max(0, wn I=
J
|||=|
(22ee- u2.22? + Vaen sign([|=[|) + erf
Un
v2
V2n
Finally, since un (u.rr
and sign(1>>])), we have the result.
Proof of Theorem 4.1 We have:
E. ump
[Ku,, (( =
e- w-
1/22
/|0111/202
max (0, (w, x) )dw
V2T0²
2 z0x)
v V2
e 2202+29²
c
max (0, (w, x)) (Lemma B.1
R"
lul2
- 1|10-=' 1/0- /|2/ max (w, x) ddw.
22
2
Rm
Applying Lemma B.4
obtain:
[K(u, w) x)w,wx)
wpp
"
11013
<||=c|1 2
V25e -<-1i=1 + / (u', at)
2202-252
(V2T(2
+ erf
(v2c1100]
2
||=||
-n/2 -11013
||=||| T
(w) 22
(u',2))
1 erf ( (x),
((2(11211 )1)
(1+ x0)
V2ce
+ u
||
C
Theoretical constants of RWRelu
For each of the two instantiations Dubé Marchand (2024) introduce, they provide the value two theo- retical constants, x and 0, relevant for their theoretical guarantees. For given instantiation (W, K,p), have:
sup xEX
E wapp
| x)K(w, -| )133] xEX wap
sup
E [K(w, [((,
(44)
sup E 2' wrpp
[o(w, (0(w,x))K K)
sup #E V
E w-pp w--p
E [K(u, w)p(u, 2) ()),,, x:).
(45)
For completeness, we give here the exact value of к and an upper bound for A, and the relevant calculus, in the case of instantiation RWRelu. We begin by couple of lemmas to help with the calculus. Lemma C.1. Consider a > 0 and y > Then:
7200 n/2
[ RD R"
e (/lu-w|///21 (u,w)//? dudw = (2mrm
(46)
20 xz
30

---

Published in Transactions on Machine Learning Research (01/2025)
Proof.
[0,w]
c
e
dudw
c
2
dtdw
(t (t:= u
w, dt = du)
R2
2"
R
1
1wil?
27
dldw
11
1ell
e
25
dt
dw
1|
02
2y
27
dc
dw
11
82
2n
27
dt
dw
R
R
2uu
x2
-V22T7²
1w
22
e
2n
dw.
Ra
Then, simplifying the exponent:
w
o²
|wl
72
2
2x2
2r
0²
||/1|
1202 72
20²
02
|w|l
(20² -22),
2
wC get:
1N-w 12
1w12 (202_g2) dw
E
2y2
dudw
Jnn
= (V2mm²²
04
2T
20² 7?
7204 20
n/2
-2T)"
72
Lemma C.2. Consider > 0 and 7> > 0. Denote In the üdentity matria in R Then:
[ (lllu--sl) ///l////?? = -(1+ 20
-n/2
E 201 (0,
(47)
22
31

---

Published in Transactions on Machine Learning Research (01/2025)
Proof. The expectation is straightforward integral:
2n
N
[uj2
E E [K(,, w)] =
24
20
2n
dudw
wmpp ump
v 2
R
R
2n
1ul
Su
E
24
2
0
20
e
dudw
2 2
B
R
2n
Iu-wi
e.
24
dudw
2m
R
R"
v 2TT
IN-wi
762
22 dudw
2
22422
20
(²0
n/2
(2TH) 10
(Lemma C.1
v 2TO
20 (2)
2
12/
20 2
20
-n/2
-1 1
20
20
-n/2
-1 1)
02
7/2
20²
-1/2
72
Lemma C.3. Considering RWRelu, we have:
0
sup | ZE
v2
Proof. We have:
n² = sup E [[w,, (w),,,)))
zEx wap
sup E
max (0, (ur, x:))
x wmp
/1/w011//20
max (0, (w, x))
sup zEX
V2Ra²
Ra
-sup zEd V2Ta²
e w /202 max(O, wn
(wn == (w.r) ||l
see proof of Lemma B.4
/202 wn | dw
=sup =EX 2
-sup zEX
||x]|
2
1
e w2z www
22 2
1 sup
||ا||
E we
[w]]
xex o²
2
(0,00)
sup ||
2 ZEx
We have the result by taking the square root.
32

---

Published in Transactions on Machine Learning Research (01/2025)
Lemma C.4. Considering RWRelu, we have:
20²
-(n-1)/2
022
0² S sup
||>||
1
(48)
=E
7/
2π
Proof. We have:
8² sup I E [K(u, ()))),
sEX umpp uoup
a
sup sx
e
c-llu-wll
/2-2 -21==2
|wll /202 e-1lw12/20 max(0, (u, x)) max (0, (wr, a)) dudw.
2TO²
JR"
We are free to solve these integrals using any orthonormal basis of R We choose a basis (vi,... Uo} such that:
U/
(49)
|
In particular, we have:
Un
u.
[2]
Wn
w,
||:0|]
(u, x²) un |
(w, x)) wi, |1:|||-
Denoting Win-1 = (ui,.... 12-1) and win- = (wi,... ... Wn 1), we can rewrite the expression:
|- /2 /2o² /20²
(50)
88:
|1 11 1/1m /2??
/202 win 111 /202 (w,)- 22/200 -2//220
(51)
This allows us to separate the integral into two parts. The first one is the integral over ( un-1) and (Wi,.... Wn-1):
Mim-i -im---Wim-1lf /2?
I|um 12/22 -1ww 112 dui - dd
(52)
2
Lemma C.2 tells that the previous expression is equal to:
I 20
1mm-1)/2
(53)
72
The second integral is the one over un and Wnx
[e e e- wwn)// 22, e w2, max(0, un un||=|]) max(0, w. w.ll=|])d,,
(54)
2
We can simplify the integral by noticing that max(0, un)) Un if u, and 0 otherwise:
1000 (wnun- //=² e2222- 22,1200 22, un du,,
2To² |
(55)
The expression for 8 has now become:
0² sup =EX
cllu-ww1
- /122 /20² -11lw| - max(0, (u, x\) max(0, (w, x²) du)d
2T
JRm
sup
||>]||
1+
20²
-(n-1)/2
1
10
(u., (u,,,w,)) w /22 - /20² - w2,/20
/20 Un w du, dwn
(56)
72
2T
#Ex
33

---

Published in Transactions on Machine Learning Research (01/2025)
The remaining integral is difficult. We can simplify by using the fact that e e 120 S 1, and obtain upper bound on 0²:
0 sup =EA'
||0||
1+
20²
-(n-1)/2
1
00
00
[un [u.-1111125
20.//20 e euu,t20 un in
7/2
2
0
sup =ex
||0||
1+
2
-(n 1)/2
1
- 222/20" e
un un
7/2
2Tσ²
sup =eX
||
1 +
2
(n- 1)/2
2 )
1
(Foo
22//00" Wnddwn
2
7
sup #EX
||x;||
1+
20
-(n 1)/2
(-ope
u2//20" 2002000
2
72
2
10
20²
111-1)/2
sup =EX'
||>||
+ 7/2
a
2T
1 20 7/2
1)/1
0²
sup ||0|| zEN
2m
Theoretical constants of STAR RKHS Weightings
We showed in Section how to transform an RKHS Weighting instantiation into STAR model, while maintaining the form of the model as an RKHS Weighting (Table contains the details). In particular, all
algorithms and theoretical guarantecs from Dubé Marchand (2024] still apply. The guarantees involve the thcoretical constants:
sup ZE
E wmp
[1o(w, x))(w,, 1[
E ([W,,)) =EX wapp
sup
(57)
sup x
E wapp
(o(w K(,,
sup ZEX
E E wep wep
((,,, x) x)];
(58)
H
defined for any instantiation (W, o, K)
Here, we explain how obtain those constants, if desired, for a STAR RKHS Weighting instantiation. For K, have:
x²2 sup zEX' (w. Pwrp PP
[K((w, 1), (w, (w, 1)) x1)
sup EE E
[Kw(w, w) 1| = (1r x7)
ZEX' Imp(a) wWPwiz
[Kw (w, w) w))r(((())) (w, =]
sup
zex (mP()) an-Pwjr
< sup ImPpa sEA w~Pwir
[Kw (w, (w 21)2
The expression supzex Enersppis [Aw (w, w) (w, my) in fact the constant x² for the partial instantiation defined only on the feature subset I. Let us denote it n?. Then we have:
I>Pl KJ.
E
(59)
This quantity depends on the distribution of the feature subsets P(n]. In the simplest cases, where for example Pin is the uniform distribution over all feature subsets of size k then the constant, Ky might simply be the
34

---

Published in Transactions on Machine Learning Research (01/2025)
same for all I. In the worst cases, take the supremum of Ky over all possible subsets I instead of the expectation. As for 0:
0²
sup F
[K[(u, J), (rc, ( (n (1)
# (w, (rrpp P(d (a. (ap)) Xppl
sup E
E
[Kw(u, w 1 = (0)
=EX (w:n) P(o] (a.J) XPpl
sup E umpwi
E
[Pin (0
xE (m:
sup E (1)
E
E
[Kw(u, ()
=EX [mPp))]
w>Pwi u-pwirr
< E (1) sup E ImPpu
E
[Kw (u, (1)
xEx
w-Pwi wwpwi
Exactly before, the expression supeex Foe- Eur-pwir (T) or )or))r((() (w, is the con- stant 0 for the partial model defined on the subset I. We can write it as 0f and:
8 E (T) T~Ppa]
(60)
In the event that 0, s 0 for all I for some 0, and that P(n)(T) is constant, then we have:
g² Plu (I)
(61)
or:
0 S VPbi (T)O.
(62)
E
Details of experimentation
Preprocessing of the datasets
All datasets have been scaled have mean and standard deviation 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 of Dubé & Marchand (2024), the least squares fit of the coefficients. However, have run into some rare numerical instability problems which led to failure of learning and even negative training R² scores, which is impossible in theory (at worst, the algorithm could output 0 E H, which has an R² of 0 exactly). To solve this issuc, we added a very small & regularizer to the optimization objective (Equation 53 of Dubé & Marchand (2024)), which is now:
TG X a Ga + 1 Da y12 + 2m
(Aa(z))
2 T.a
(63)
where I is the identity matrix, and is the 8 regularization parameter, which we have to 10-10. The minimizer of the previous expression is the solution to the linear problem:
o T + mXG + +mella = TTy
(64)
Hyperparameter selection
As pointed out in Dubé & Marchann (2024), instantiation RWSign, and now also RWRelu, exhibits ex- ponential behavior in the dimensionality of the instance space. This immediately be seen in the
n/2
c: (1+
coefficient in Equation 25 shared by both Rigi and RWRelu. Instead of using o
35

---

Published in Transactions on Machine Learning Research (01/2025)
Cross-validation parameters
Source code (clickable)
Decision Reesssrr
max depth {2, 10, 20}
Scikit-learn
max bins E {512, 1024, 2048] learning rate E (512, 1024, 2048) max rounds {15000, 25000, 35000} min samples leaf E [1,2,3}
EBM
Interpret
kernel rbf
RernelRidge
Scikit-learn
alpha {0.01, 0.05, 0.1, 0.5, 1,5]
LinearRegression
Scikit-lcarn
SVR
C E {0.5, 1,5, 10,50)
Scikit-learn
max theta e (0.1, 0.5, 0.9) A [10-9, 100 100000000000 7,10-b,
RWSign
Dubé & Marchand (2024
1,00,,10-- 10 10-2}
max theta E {0.1, 0.5, 0.9}
RWRelu
This paper
A (10-9,101) 10-9, 10-1, ,10 10
10 10 10
(0.01, 0.1,
RWStumps
(0.01,
Dubé & Marchand (0024
A E {10-9) 10-8, 100 10 8, 10 1,- 100 10-3, 10103
Table 9: Algorithms and models used in this paper and their hyperparameters.
and y as hyperparameters of the model, they use a and c, and calculate from a and c. This way, the effect of dimensionality is held constant regardless of n.
We modify this methodology slightly. We instead use c² : (1 20
n/2
This is the upper bound for 0²
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 0 for RWRelu (Lemmn C.4h. Given a and c E (0,1), 1), we get:
2
(65)
0²(c 4/n -/ 1)"
The advantage of this method is that the model output Aa(x) being upper bounded by 0||a|y I ccla this gives meaning to the parameter For RWSign, this allows setting the maximum value of the model. For RWRclu, the relationship is functionally the same, though the upper bound for 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. httpp://ggttuu..oo//aaat///////3 http://///.aaa///aaaaaaaaaeeeeee seiigtttrrrrttts-eetttttnggggg The repository contains 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 TO lists the commands to run the experiments of Section 5 and the files they generate.
36

---

Published in Transactions on Machine Learning Research (01/2025)
Command
Files generated
resultts/shapley_tiin
Fig.
python -m majjorexperiments..oe shapheys_hhppypypppaaa time
Fig.
python major_ www.eeeeeeesessaaelee... stt.aaeacccpacaayaaacyr
reeutss/rggreession///n....ee
Table
python -m major_ wwwpeeeeees.eeeeeeee.eeeeeeeeeee --final
talls//rggresstoo0//nll1...xx
aable//eggesslonn/iiaadddeeeeee
reautts/agrression///iaa..eerrr
Table
python -m cenjjo_experimentt time_serice --final
talles/tiin--rrrree-ffiial.eee
cable//iieererree-ffaadddddas tex
reaulta/eea..epllaaattii....
reuuta//thee.xxplaaatttii...pp pk1
figuree/beswvrrr_pooteeem pdff
Fig. 3. 4
python major_ maajor_expeeeeeeee wwwprreeetssteeeeeeeee chapley_compartee
figuressbbeswarr_ploo_.te pdf
tigures/aaterral__oot_too pdf
Figure//aaterfal__ooo__aa pdf
Table 10: Commands to run the various experiments in this paper.
37


Reviews:
[{"review_id": "ny697MJAKt", "reviewer": "Reviewer_nmfe", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I think the article requires significant revisions to meet academic standards, and I propose a major revision for the following motivations. \n\n1. The text contains unclear sentences and at times employs overly simplistic or even confusing English. For example, after Equation (6), the authors write [...] the set of variables in $\\pi$, but $\\pi$ is a function to create permutations and it is not a set. I also understand that $U(\\Omega)$ is the Uniform distribution, but this should be stated explicitly. Also, if $S$ is a dataset and not a random variable, what does it mean $z\\sim S$ in the expression of the expected values? Moreover, the use of contractions is not suitable for the academic tone required for this publication.\n2. The language and explanations are not accessible to readers who are not already experts in the subject. The number of acronyms present in the manuscript is excessive and distracts the reader. For example, the terms *RKHS, ReLU* (even if their meaning can be guessed), are introduced without explanation.\n3. The section on Shapley values lacks sufficient detail and fails to clarify fundamental aspects. \n\n    a. The value function used in the Shapley value calculations is not explicitly defined or explained. An explicit description of this function and its role in interpreting Shapley values is necessary.\n\n    b. The text does not adequately discuss the Shapley values in the practical applications. The datasets are presented only in terms of sample dimension. Basically, the applications are presented as: \"These are the datasets, these are the $R^2$ of the various models, in Figure 1 there is the computational cost of the Shapley values of a Structured Additive Regression model with respect to the dimensionality of the dataset (which one?).\" I really find the part on the Shapley values quite poor, but if I look at the title and the abstract I expect this to be the core of the paper.\n\n    c. Page 5 - line 12 from the bottom: *[...] then this formula can be calculated in $O(mN 2^k)$, independently of the dimensionality of the data.* Since $m=|S|$ where $S$ is the dataset and $N$ is the number of features which satisfy a certain condition (with $k$ as maximum number), it is unclear why the complexity is independent of data dimensionality. Not only, is it really necessary to introduce $m$ as a new notation?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The purpose of the paper is twofold: the authors find a closed form expression to compute the Shapley values of a family of models, called Structured Additive Regression (STAR) models, and presents an algorithm for their calculation. This algorithm allows to compute in a reasonable amount of time the Shapley values, if the maximal number of variable interactions is small. The major part of the paper is devoted to providing connections between RKHS models and STAR models. Finally, the paper is completed with applications on various datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. The paper *Khorrami Chokami, A., \\& Rabitti, G. (2024). An Exact Game-Theoretic Variable Importance Index for Generalized Additive Models. Journal of Computational and Graphical Statistics, 1–10. DOI: 10.1080/10618600.2024.2327577* provides closed-form results for Shapley values for GAMs, along with confidence intervals. Please clarify and contrast your contribution (for GAMs) with respect to that paper.\n\n2. Concerning the application, consider to add the computed Shapley values to interpret the results. This would improve the overall interpretation and relevance of the findings.\n\n3. Cite *Owen, A. B. (2014), Sobol’ Indices and Shapley Value, SIAM/ASA Journal on Uncertainty Quantification, 2, 245–251. DOI: 10.1137/130936233* and *Wood, S. N., Goude, Y., and Shaw, S. (2015), Generalized Additive Models for Large Data Sets, Journal of the Royal Statistical Society, Series C, 64, 139–155. DOI: 10.1111/rssc.12068*}.", "broader_impact_concerns": "It is not clear from the article to what extent the proposed models are actually used and considered relevant in practical applications.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "PBcMAf8CXi", "reviewer": "Reviewer_au4X", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n\n1.\tThis paper proposed an algorithm for computing Shapley values for a family of models, i.e., structured additive regression models, which were the generalization of generalized additive models (GAMs).\n\n2.\tThe authors presented an instantiation of RKHS Weightings and showed how to transform an RKHS Weighting model into a STAR model.\n\n**Weaknesses:**\n\n1.\tThis paper should explicitly compare the time complexity of approximating Shapley values for all models versus STAR models, i.e., the time complexity of Algorithm 1. To highlight the contribution, the authors could list commonly used STAR models below Eq. (11) to clarify the application scenarios for this algorithm. Additionally, a table comparing typical use cases for Shapley values with RKHS-weighted models, STAR models, and all models would be helpful, along with the computational costs for full n-way interactions and k-way interactions (where $k \\le n$).\n\n2.\tThe purpose of introducing RKHS weighting in Eq. (12) is unclear. Since the authors referenced an unpublished paper [Dub. & Marchand (2024)], it would be beneficial to provide background on RKHS, RWSign, and RWStumps. Specifically, the authors should clarify why RKHS weighting is needed in the context of Shapley values for STAR models, the model family it represents, the typical models it includes, and its potential applications. Additionally, the authors should explain the fundamental differences between RWSign and RWStumps, and specify the scenarios for which each method is applicable.\n\n3.\tAlgorithm 1 did not essentially reduce the exponential complexity of Shapley values. Instead, it reduced the computational cost by constraining the number of interacting variables in the model. The authors should note this in the introduction. Experimentally, the authors should evaluate models with interactions involving more than just 5 variables in Figure 1. For example, in Section 6, they could include models with 1 to 10-way interactions, and report the model’s performance, providing an illustration of the algorithm’s limitations.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposed an efficient method for computing Shapley values of structured additive regression (STAR) models. Additionally, the authors introduced a novel instantiation of RKHS Weighting that was better suited for regressions, transforming it to fit within the STAR model framework.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "No concerns.", "broader_impact_concerns": "See weaknesses.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Bg0nGDHIDF", "reviewer": "Reviewer_7iuM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n\nThe paper is very clearly written. The paper is also appropriately positioned with respect to the related literature, with the exception of some recent developments in estimating shapley values (see below).\n\nTheorem 3.1 and the proposed algorithm could potentially be interesting to some readers of TMLR. \n\nGenerally, reading this paper could be interesting as an introduction for researchers who are not that well-aware of the connections between Shapley Values and functional decompositions. \n\n**Weaknesses:**\n\nThe proposed algorithm is a simple consequence of the structure of the model. I don't think the shapley-estimation part of this paper contains much novelty for researchers who are aware of the connections between Shapley Values and functional decompositions.\n\nThe paper contains no real performance comparison with other methods that estimate shapley values, except for SHAP. For a paper that introduces a new algorithm that is performant, I would expect a somewhat more extensive performance comparison with other algorithms. \n\nIn the experiments, it seems that EBMs perform pretty well, which makes sense because the EBM software package is fairly optimized for the application scenario considered in this paper. So why do we need the newly proposed method?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes the \"STAR-SHAP\" algorithm for computing shapely values for \"Structured Additive Regression Models\" (STAR models).\n\nAs the authors note in the paper, STAR models are already known in the literature under different names, such as GAMs with arbitrary interactions or functional decompositions.\n\nThe STAR-SHAP algorithm is a model-specific algorithm. It assumes that we can access the model's different feature-wise component functions. The Shapley values are then given by the expectation of a weighted sum involving the component functions (equation (16) in the paper).\n\nThe advantage of the STAR-SHAP algorithm is that it does not involve a sum over exponentially many terms (as in the Shapley formula) but only over the terms present in the model.\n\nThe paper also shows that particular classes of models from Dubé & Marchand (2024) can be cast as STAR models.\n\nThe paper has several empirical applications, all regression experiments.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The paper needs to discuss/compare with other more recent algorithms to compute shapley values for black-box models. For example: https://proceedings.neurips.cc/paper_files/paper/2023/hash/264f2e10479c9370972847e96107db7f-Abstract-Conference.html  \n\n- it would be great if the paper could include a table of different model-specific and model agnostic algorithms in the literature and state the computational complexity of the various algorithms. This would help the reader to situate the paper's result in the literature.\n\n- I would appreciate if the paper could elaborate a bit about the con's and pro's of model-specific / model-agnostic shapley values computation. \n\n- The paper should include examples of Shapley values of STAR models estimated with different methods: the STAR-SHAP algorithm, Kernel SHAP, and perhaps other model-agnostic approximation algorithms?\n\n- Because you use EBMs, you should cite the interpret software package https://arxiv.org/pdf/1909.09223\n\n\n**Questions:**\n\n**Q1:** What happened to the value function? Does the STAR-SHAP algorithm compute \"interventional\" shapley values, or does that depend on the type of the expectation operator in equation (16)? I'm asking this because some forms of shapley values (with marginal distribution) are usually much harder to estimate, and model-specific algorithms like treeSHAP usually implicitly rely on a specific value function. It would be great if the authors could clarify this in the paper.\n\n**Q2:** The connection between the RKHS methods and STAR models seems interesting, but to me, it also seems somewhat unmotivated. Why consider this particular, somewhat unusual class of models? How is this connected to the STAR-SHAP algorithm? Is the main reason that you need a method to obtain non-trivial STAR-models?\n\n**Q3:** Perhaps you could briefly explain the connections/differences between equation (16) in your paper and equation (11) in Bordt & von Luxburg (2023) (which appears in many other papers, too).", "broader_impact_concerns": "Theory / Algorithms paper; broader impact statement seems uneccesary.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers recommended acceptance (at least weakly) after the rebuttal and all indicated explicitly that their concerns were addressed.

I do note that Reviewer nmfe suggested adding a table of acronyms, which the authors should consider for the camera-ready revision.
Decision: Accept (poster)