🤖 machine learning

Amortized Linear-time Exact Shapley Value for Product-Kernel Methods

This paper introduces PKeX-Shapley, a novel algorithm that leverages the multiplicative structure of product kernels to compute exact, parameter-free Shapley values for all features in amortized linear time, thereby overcoming the computational intractability and estimation errors inherent in existing approximation methods for kernel-based explainability and statistical analysis.

Original authors: Majid Mohammadi, Siu Lun Chau, Krikamol Muandet

Published 2026-05-08

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Majid Mohammadi, Siu Lun Chau, Krikamol Muandet

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Big Problem: The "Black Box" and the "Impossible Math"

Imagine you have a very smart, but mysterious, AI model. It's like a black box that takes in a bunch of ingredients (features) and bakes a cake (makes a prediction). You want to know: Which ingredient made the cake taste the way it did? Was it the sugar? The flour? The vanilla?

In the world of AI, we use a mathematical tool called Shapley Values to answer this fairly. It's like a game where you try every possible combination of ingredients to see how much each one adds to the final flavor.

The Catch: If you have 10 ingredients, there are 1,024 combinations to check. If you have 50 ingredients, there are more combinations than there are atoms in the universe.

The Old Way: To get a "good enough" answer, people usually guess by sampling a few combinations. It's fast, but it's an estimate, and it can be wrong, especially when you have many ingredients.
The Goal: We want the exact answer, not a guess, and we want it to happen fast, even with hundreds of ingredients.

The Solution: PKeX-Shapley

The authors introduce a new method called PKeX-Shapley. Think of this as a "magic shortcut" that works specifically for a certain type of AI model called a Product-Kernel Method.

1. The "Multiplicative Team" Analogy

Most of these models work like a team of specialists where the final result is the product (multiplication) of their individual contributions.

Imagine a recipe where the final taste is: (Salt Factor) × (Sugar Factor) × (Spice Factor).
In math, this is called a Product Kernel.

The authors realized that because these models multiply things together, they have a special property: If you remove an ingredient, you don't need to re-bake the whole cake. You just replace that ingredient's factor with a "neutral" number (the number 1).

Example: If you remove the "Spice Factor," you just multiply by 1. The math stays simple and clean.
Why this matters: Old methods tried to simulate "removing" an ingredient by looking at other data or guessing what the missing data would be. This new method just says, "Let's pretend this ingredient is a neutral 1." It requires no guessing, no sampling, and no extra data.

2. The "Assembly Line" Trick (Speeding Up)

Even with the "neutral 1" trick, calculating the exact contribution for every single ingredient usually takes a long time (exponential time).

The authors found a way to organize the math like a factory assembly line.

Instead of calculating the contribution of Ingredient A, then Ingredient B, then Ingredient C separately (which is slow), they realized that the calculations for A, B, and C share a lot of the same "building blocks."
They built a system (using something called Elementary Symmetric Polynomials) that calculates all these shared blocks once and then reuses them for every ingredient.
The Result: Instead of taking hours or days for 1,000 ingredients, their method takes seconds. It scales linearly, meaning if you double the ingredients, you only double the time, not square it.

What Can It Do? (According to the Paper)

The paper claims this method works for three main things:

Predictive Models: It can explain why a model made a specific prediction (like a Support Vector Machine or Kernel Ridge Regression) by telling you exactly how much each feature contributed.
Comparing Distributions (MMD): Imagine you have two groups of people (Group A and Group B). You want to know why they are different. This method can tell you exactly which features (like age, income, or height) are driving the difference between the two groups.
Measuring Dependence (HSIC): Imagine you want to know if two things are related (e.g., "Does the weather affect ice cream sales?"). This method can break down that relationship to show you exactly which weather factors (temperature, humidity, wind) are responsible for the connection.

The "Catch" (Limitations)

The paper is very honest about its limits:

It only works for "Product" models. If your AI model mixes ingredients in a complex, non-multiplicative way (like a deep neural network with tangled layers), this specific shortcut doesn't work.
It's exact, but specific. It trades the ability to work on any model for the ability to be perfectly accurate and fast on this specific type of model.

Summary in a Nutshell

The Problem: Explaining complex AI models is usually slow and full of guessing.
The Innovation: The authors found a mathematical "cheat code" for models that multiply their inputs together.
The Magic: By treating "missing" ingredients as a neutral "1," they avoid all the guessing and sampling.
The Speed: They built an assembly line to calculate the answers for all ingredients at once, making it fast enough to handle thousands of features without losing accuracy.
The Outcome: You get a perfectly fair, exact breakdown of what matters, whether you are predicting a number, comparing two groups of data, or checking if two things are related.

Technical Summary: Amortized Linear-time Exact Shapley Value for Product-Kernel Methods

Problem Statement
Kernel methods are widely adopted in machine learning and statistics for their flexibility and expressive power, yet their "black-box" nature hinders adoption in high-stakes applications. While Shapley value-based attribution methods (e.g., SHAP, RKHS-SHAP) offer a principled framework for explainability, exact computation of Shapley values is generally intractable due to the exponential number of feature subsets ( $2^d$ ). Existing approaches rely on approximations (e.g., Monte Carlo sampling, regression-based estimators) that incur unavoidable estimation errors, which grow with dimensionality. Furthermore, standard value functions often require sampling from marginal or conditional distributions or density estimation, introducing additional computational burdens and dependencies on background data.

Methodology
The paper introduces PKeX-Shapley, an algorithm designed to compute exact Shapley values for all $d$ features of product-kernel methods (e.g., Support Vector Machines, Kernel Ridge Regression) in quadratic time with respect to the number of features ( $O(d^2)$ ), or amortized linear time per feature.

The methodology rests on three core components:

Distribution-Free Removal Operator: The authors focus on the subclass of models employing product kernels, where the kernel function factorizes as $k(x, x') = \prod_{j=1}^d k_j(x_j, x'_j)$ . They identify a natural removal operator intrinsic to this structure: removing a feature $j$ corresponds to replacing its kernel factor $k_j$ with the multiplicative identity ($1$). This choice preserves the product structure, requires no background distribution, and induces a parameter-free value function.
- Unlike interventional or observational value functions that require sampling, this value function is defined via functional decomposition.
- The value function for a coalition $S$ is given by $v_x(S) = \alpha^\top k_S(X_S, x_S) - f_\emptyset$ , where $f_\emptyset$ is a constant term that cancels out in marginal contributions. This formulation requires no sampling or density estimation.
Functional Decomposition: Leveraging the removal operator, the model $f(x)$ is uniquely decomposed into a sum of components indexed by feature subsets: $f(x) = \sum_{S \subseteq D} f_S(x_S)$ . The components are derived via Möbius inversion on the subset lattice, yielding a closed form: $f_S(x_S) = \sum_{i=1}^n \alpha_i \prod_{j \in S} (k_j(x_j, x^{(i)}_j) - 1)$ .
Recursive Computation via Elementary Symmetric Polynomials (ESPs):
- The Shapley value for a single feature is expressed as a weighted sum of marginal contributions, which factorizes into terms involving Elementary Symmetric Polynomials (ESPs) of the kernel evaluations.
- Direct computation of ESPs for every feature would be inefficient. The authors derive a recursive formulation based on Newton's identities.
- To compute all $d$ Shapley values jointly, they employ a prefix–suffix construction. They define generating polynomials for the prefix (features before $j$ ) and suffix (features after $j$ ) of the feature set. By multiplying these polynomials, they recover the ESPs for the leave-one-out set $Z_{-j}$ .
- This allows the computation of all Shapley values in $O(d^2 n)$ time, where $n$ is the number of samples. The algorithm is numerically stable, relying only on element-wise additions and multiplications, avoiding the numerical instability associated with polynomial interpolation methods (e.g., Chebyshev or Fourier interpolation) used in other exact methods like TreeSHAP when applied to high dimensions.

Extension to Statistical Discrepancies
The framework is extended beyond predictive modeling to kernel-based statistical discrepancies, specifically:

Maximum Mean Discrepancy (MMD): Used to measure distributional closeness. The removal operator yields a value function that allocates the overall discrepancy between two distributions to individual variables.
Hilbert–Schmidt Independence Criterion (HSIC): Used to measure dependence between random variables. The method provides a functional decomposition of the HSIC, allowing for the attribution of dependence to specific features.

Results
The paper validates PKeX-Shapley through experiments on synthetic and real-world datasets:

Accuracy vs. Approximation: In local explanations, PKeX-Shapley (exact) significantly outperforms sampling-based approximations (e.g., Kernel SHAP, RKHS-SHAP). As dimensionality increases (e.g., $d=50$ ), sampling-based methods require prohibitively large coalition sizes to achieve low error, whereas PKeX-Shapley provides exact results regardless of dimension.
Feature Selection: In global sensitivity analysis using HSIC, PKeX-Shapley effectively identifies active features in synthetic regression tasks with polynomial and exponential target functions, outperforming or matching baselines like RKHS-SHAP, GEMFIX, and BiSHAP.
Computational Efficiency: The algorithm scales to 1,000 features in seconds, whereas brute-force computation is limited to ~20 features. It is consistently faster than sampling-based methods, even when those methods use a small number of coalition samples.
Numerical Stability: Experiments comparing different ESP aggregation strategies (quadratic recursion vs. FFT/Chebyshev interpolation) show that while interpolation methods are fast, they suffer from catastrophic numerical instability at moderate dimensions ( $d \approx 50$ ). The proposed quadratic recursion remains stable up to $d=1000$ .

Significance and Claims
The paper claims that PKeX-Shapley addresses the two primary challenges of Shapley value computation for kernel methods:

Exactness without Approximation: It provides a closed-form, parameter-free value function that eliminates the need for sampling and density estimation, thereby removing estimation error.
Scalability: It achieves amortized linear time per feature ($O(dn)$) for exact computation, a significant improvement over the exponential cost of brute-force or the high variance of sampling methods in high-dimensional settings.

The authors emphasize that while the method is restricted to product kernels, this structural constraint is necessary for tractable exact computation. The framework's extension to MMD and HSIC provides new tools for interpretable statistical analysis, enabling the attribution of distributional differences and dependencies to specific variables without relying on heuristic approximations. The authors note that the method does not naturally extend to stochastic attribution (e.g., variance of Shapley values in Gaussian Processes) in the same way methods for FANOVA GPs do, as orthogonality does not generally hold for product kernels.