diff --git a/.gitignore b/.gitignore index 641786cb8..d6f575f38 100644 --- a/.gitignore +++ b/.gitignore @@ -181,6 +181,7 @@ illustration_sources/ src/shapiq_games/datasets/data/tabarena_*.csv shapiq/games/benchmark/precomputed/ +benchmark/results/* precomputed.zip game_storage/* @@ -189,6 +190,7 @@ src/shapiq/_version.py # Claude .claude/ +.claude # C-Extension files *.so @@ -199,5 +201,8 @@ docs/source/sg_execution_times.rst docs/source/generated/ docs/source/gen_modules/ +# Cache +cache/ + # Local debug / scratch scripts (not public) scripts/ diff --git a/CHANGELOG.md b/CHANGELOG.md index 07a05ce04..5808db30f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # Changelog +## Unreleased + +### New Features + +- adds the `PolySHAP` approximator in `shapiq.approximator.regression` for Shapley value estimation via interaction-informed polynomial regression (PolySHAP, [Fumagalli et al., ICLR 2026](https://arxiv.org/abs/2601.18608)). PolySHAP generalizes `KernelSHAP` by fitting a *k*-additive polynomial surrogate of the game before reading off the Shapley values. The interaction frontier can be the full *k*-additive frontier up to `max_order` (with `max_order=1` recovering `KernelSHAP`), a budget-controlled partial frontier (`max_terms`), or a user-supplied set of interaction terms (`prior_frontier`). Registered as a Shapley-value (`SV`) approximator. + ## v1.6.0 (2026-07-06) ### Highlights of new Features diff --git a/docs/source/references.bib b/docs/source/references.bib index bfc2afafd..fd37bf2eb 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -111,6 +111,13 @@ @article{Fumagalli.2026 eprinttype = {arXiv}, eprint = {2602.01399} } +@inproceedings{Fumagalli.2026a, + title = {{PolySHAP}: Extending {KernelSHAP} with Interaction-Informed Polynomial Regression}, + author = {Fabian Fumagalli and R. Teal Witter and Christopher Musco}, + year = {2026}, + booktitle = {The Fourteenth International Conference on Learning Representations ({ICLR} 2026)}, + url = {https://arxiv.org/abs/2601.18608} +} @inproceedings{Harris.2022, title = {Joint Shapley values: a measure of joint feature importance}, author = {Chris Harris and Richard Pymar and Colin Rowat}, diff --git a/examples/approximators/plot_polyshap.py b/examples/approximators/plot_polyshap.py new file mode 100644 index 000000000..c4e452fdc --- /dev/null +++ b/examples/approximators/plot_polyshap.py @@ -0,0 +1,76 @@ +""" +PolySHAP +======== + +Shapley value approximation with interaction-informed polynomial regression +using :class:`~shapiq.approximator.PolySHAP` :footcite:t:`Fumagalli.2026a`. + +PolySHAP extends KernelSHAP by fitting a *k-additive* surrogate of the game -- a +polynomial of degree ``max_order`` over the players -- and reading the Shapley +values off that surrogate. ``max_order=1`` reduces to plain KernelSHAP, while +higher orders capture interactions explicitly. When the game's interactions do +not exceed ``max_order``, PolySHAP recovers the exact Shapley values. +""" + +from __future__ import annotations + +import numpy as np + +import shapiq +from shapiq.approximator import PolySHAP + +N_PLAYERS = 8 +BUDGET = 200 +feature_names = [f"x{i}" for i in range(N_PLAYERS)] + +weights = np.array([0.4, 0.3, 0.2, 0.1, 0.05, -0.1, -0.2, -0.3]) + + +def game_fun(coalitions: np.ndarray) -> np.ndarray: + # A linear game plus a single pairwise interaction between x0 and x1. + # This game is exactly 2-additive: its interactions never exceed order 2. + coalitions = np.atleast_2d(coalitions) + return (coalitions @ weights) + 0.5 * coalitions[:, 0] * coalitions[:, 1] + + +# %% +# Approximate Shapley values +# -------------------------- +# The game has a pairwise interaction, so we fit a second-order (``max_order=2``) +# polynomial surrogate to recover the Shapley values. + +approximator = PolySHAP(n=N_PLAYERS, max_order=2, random_state=42) +iv = approximator.approximate(BUDGET, game_fun) +print(iv) + +# %% +# Force plot +# ---------- + +iv.plot_force(feature_names=feature_names) + +# %% +# Higher order improves accuracy +# ------------------------------ +# Because the game is exactly 2-additive, raising ``max_order`` to match its +# interaction order recovers the exact Shapley values, whereas ``max_order=1`` +# (equivalent to KernelSHAP) leaves a visible approximation error. We compare +# both against the exact Shapley values from :class:`~shapiq.ExactComputer`. + +exact = np.asarray( + shapiq.ExactComputer(game_fun, n_players=N_PLAYERS)(index="SV", order=1).get_n_order_values(1) +) + +for max_order in (1, 2): + est = np.asarray( + PolySHAP(n=N_PLAYERS, max_order=max_order, random_state=42) + .approximate(BUDGET, game_fun) + .get_n_order_values(1) + ) + mse = float(np.mean((est - exact) ** 2)) + print(f"max_order={max_order}: MSE vs. exact Shapley = {mse:.2e}") + +# %% +# References +# ---------- +# .. footbibliography:: diff --git a/src/shapiq/approximator/__init__.py b/src/shapiq/approximator/__init__.py index 00a4494b9..15276985c 100644 --- a/src/shapiq/approximator/__init__.py +++ b/src/shapiq/approximator/__init__.py @@ -17,6 +17,7 @@ KernelSHAP, KernelSHAPIQ, OddSHAP, + PolySHAP, RegressionFBII, RegressionFSII, kADDSHAP, @@ -64,6 +65,7 @@ def __init__(self, *_args: object, **_kwargs: object) -> None: ProxySPEX, OddSHAP, ShaplEIG, + PolySHAP, ] # contains all SI approximators @@ -149,4 +151,5 @@ def __init__(self, *_args: object, **_kwargs: object) -> None: "STII_APPROXIMATORS", "FSII_APPROXIMATORS", "FBII_APPROXIMATORS", + "PolySHAP", ] diff --git a/src/shapiq/approximator/regression/__init__.py b/src/shapiq/approximator/regression/__init__.py index a0b7e5c1a..a8f40cc57 100644 --- a/src/shapiq/approximator/regression/__init__.py +++ b/src/shapiq/approximator/regression/__init__.py @@ -6,6 +6,7 @@ from .kernelshap import KernelSHAP from .kernelshapiq import InconsistentKernelSHAPIQ, KernelSHAPIQ from .oddshap import OddSHAP +from .polyshap import PolySHAP __all__ = [ "kADDSHAP", @@ -15,5 +16,6 @@ "InconsistentKernelSHAPIQ", "Regression", "RegressionFBII", + "PolySHAP", "OddSHAP", ] diff --git a/src/shapiq/approximator/regression/polyshap.py b/src/shapiq/approximator/regression/polyshap.py new file mode 100644 index 000000000..0fea783c3 --- /dev/null +++ b/src/shapiq/approximator/regression/polyshap.py @@ -0,0 +1,370 @@ +"""This module contains the PolySHAP approximator to compute Shapley values.""" + +from __future__ import annotations + +import warnings +from typing import TYPE_CHECKING, Any, Literal + +import numpy as np +from scipy.special import binom + +from shapiq.approximator.regression.base import Regression +from shapiq.interaction_values import InteractionValues +from shapiq.utils.sets import powerset + +if TYPE_CHECKING: + from collections.abc import Callable, Iterable + + from shapiq.game import Game + +ValidRegressionPolySHAPIndices = Literal["SV"] + + +class PolySHAP(Regression[ValidRegressionPolySHAPIndices]): + """Estimate Shapley values using the PolySHAP regression algorithm. + + Generalises KernelSHAP :cite:t:`Lundberg.2017`; the algorithm is described in + Fumagalli et al. (2026) :cite:t:`Fumagalli.2026a`. + + PolySHAP fits an interaction-informed polynomial surrogate of the game over an + *explanation frontier* (the interactions used as basis functions) and reads the + Shapley values off it. ``max_order=1`` recovers KernelSHAP exactly; higher orders + model interactions and recover the exact Shapley values when the game is genuinely + low-order. + + By default the frontier holds every interaction up to ``max_order`` (*k-additive*). + When that is too large for the available budget, cap its size with ``max_terms`` + (*partial* mode). Alternatively, ``prior_frontier`` lets you specify the interaction + terms directly, for instance when domain knowledge already identifies them. + + Args: + n: The number of players. + max_order: Maximum interaction order included in the frontier. Defaults to ``2``. + In *partial* mode it bounds the orders drawn from (``2 .. max_order``); set + ``max_order=n`` to leave that frontier unbounded. + max_terms: If set, cap the frontier at this many terms (*partial* mode). Whole + interaction orders up to ``max_order`` are included from low to high, and the + one order that does not fit in full is sampled at random; this keeps the + noise-robust lower orders complete. Must be at least ``n + 1``. ``None`` + (default) uses the full k-additive frontier. + sizes_to_exclude: Higher-order coalition sizes to omit from the frontier. + Singletons are always kept. Defaults to ``None``. + prior_frontier: An iterable of coalition tuples defining the frontier and its + column ordering (*prior* mode). Must include every singleton ``(i,)``. It + is mutually exclusive with ``max_order``, ``max_terms`` and + ``sizes_to_exclude``; passing any of them alongside it raises ``ValueError``. + pairing_trick: If ``True``, the pairing trick is applied to the sampling + procedure. Defaults to ``False``. + sampling_weights: An optional array of weights for the sampling procedure. + Must be of shape ``(n + 1,)`` and determines the probability of sampling + a coalition of a given size. Defaults to ``None``. + random_state: The random state of the estimator (also seeds the *partial* + frontier selection). Defaults to ``None``. + + Attributes: + n: The number of players. + N: The set of players (``0`` to ``n - 1``). + max_order: Interaction order of the approximation (always ``1``, since PolySHAP + targets Shapley values; distinct from the ``max_order`` constructor argument, + which sets the frontier order). + min_order: Minimum interaction order (always ``0``). + iteration_cost: The cost of a single iteration of the estimator. + explanation_frontier: The active explanation frontier dictionary. + """ + + def __init__( + self, + n: int, + *, + max_order: int | None = None, + max_terms: int | None = None, + sizes_to_exclude: set[int] | None = None, + prior_frontier: Iterable[tuple[int, ...]] | None = None, + pairing_trick: bool = False, + sampling_weights: np.ndarray | None = None, + random_state: int | None = None, + ) -> None: + """Initialize the PolySHAP approximator.""" + explanation_frontier = self._build_frontier( + n, + max_order=max_order, + max_terms=max_terms, + sizes_to_exclude=sizes_to_exclude, + prior_frontier=prior_frontier, + random_state=random_state, + ) + + super().__init__( + n, + max_order=1, + index="SV", + random_state=random_state, + pairing_trick=pairing_trick, + sampling_weights=sampling_weights, + ) + + self.projection_matrix = None + + # Every singleton must be present so that Shapley values can be read off directly. + for i in self._grand_coalition_set: + if (i,) not in explanation_frontier: + msg = "PolySHAP requires all main effects in the explanation frontier." + raise ValueError(msg) + + # Build the binary indicator matrix: rows = frontier terms, cols = players. + self.interaction_matrix_binary = np.zeros((len(explanation_frontier), self.n), dtype=bool) + self.explanation_frontier = explanation_frontier + for S, pos in explanation_frontier.items(): + self.interaction_lookup[S] = pos + self.interaction_matrix_binary[pos, S] = True + + # Exclude the empty-coalition term from the variable count. + self.n_variables = len(explanation_frontier) - 1 + + @staticmethod + def _build_frontier( + n: int, + *, + max_order: int | None, + max_terms: int | None, + sizes_to_exclude: set[int] | None, + prior_frontier: Iterable[tuple[int, ...]] | None, + random_state: int | None, + ) -> dict[tuple[int, ...], int]: + """Resolve the explanation frontier from the constructor arguments. + + The mode is selected by which arguments are supplied (see the class docstring): + ``prior_frontier`` yields the *prior* frontier; ``max_terms`` the *partial* + frontier (whole orders low-to-high, then a random sample of the boundary order); + otherwise the deterministic *k-additive* frontier up to ``max_order`` (default 2). + + Args: + n: The number of players. + max_order: Maximum interaction order of the frontier (``None`` defaults to 2). + max_terms: Optional cap on the number of frontier terms (*partial* mode). + sizes_to_exclude: Higher-order coalition sizes to omit. + prior_frontier: Optional explicit frontier (*prior* mode). + random_state: Seeds the random selection in *partial* mode. + + Returns: + A dictionary mapping each coalition tuple to its column index, with the empty + set at index ``0`` and every singleton present. + + Raises: + ValueError: If ``prior_frontier`` is combined with ``max_order``, + ``max_terms`` or ``sizes_to_exclude``, or if ``max_terms`` is smaller + than ``n + 1``. + """ + # Prior mode: use the caller-supplied frontier verbatim; it admits no other knobs. + if prior_frontier is not None: + if max_order is not None or max_terms is not None or sizes_to_exclude is not None: + msg = ( + "PolySHAP: 'prior_frontier' is mutually exclusive with 'max_order', " + "'max_terms' and 'sizes_to_exclude'." + ) + raise ValueError(msg) + return {S: pos for pos, S in enumerate(prior_frontier)} + + if max_order is None: + max_order = 2 + + def _excluded(size: int) -> bool: + return sizes_to_exclude is not None and size in sizes_to_exclude + + # The empty set and all singletons are always present so that Shapley values can + # be read off directly. + frontier: dict[tuple[int, ...], int] = {} + for S in powerset(range(n), max_size=1): + frontier[S] = len(frontier) + + higher_order = ( + S for S in powerset(range(n), min_size=2, max_size=max_order) if not _excluded(len(S)) + ) + + # Deterministic k-additive frontier: take every term up to max_order. + if max_terms is None: + for S in higher_order: + frontier[S] = len(frontier) + return frontier + + # Budget-capped partial interaction frontier of exactly max_terms terms + # (Fumagalli et al. 2026, Sec. 4): fill whole interaction orders from low to high, + # then randomly sample the single boundary order that does not fit in full. Lower + # orders are kept complete because they occur more frequently and are less + # sensitive to sampling noise. + if max_terms < n + 1: + msg = ( + f"PolySHAP: 'max_terms' ({max_terms}) must be at least n + 1 ({n + 1}) " + "to include the empty set and all singletons." + ) + raise ValueError(msg) + + rng = np.random.default_rng(random_state) + for order in range(2, max_order + 1): + if _excluded(order): + continue + remaining = max_terms - len(frontier) + if remaining <= 0: + break + order_terms = list(powerset(range(n), min_size=order, max_size=order)) + if len(order_terms) <= remaining: + for S in order_terms: # whole order fits: include it deterministically + frontier[S] = len(frontier) + else: # boundary order: randomly select the terms that still fit, then stop + rng.shuffle(order_terms) + for S in order_terms[:remaining]: + frontier[S] = len(frontier) + break + return frontier + + def _warn_if_underdefined(self, n_sampled: int, budget: int) -> None: + """Emit a UserWarning when the least-squares system has more variables than samples. + + Args: + n_sampled: Number of sampled coalitions (excluding the border pair). + budget: The budget originally passed to :meth:`approximate`. + """ + if n_sampled < self.n_variables: + warnings.warn( + f"The least-squares system is underdefined: {n_sampled} sampled coalition(s) " + f"but {self.n_variables} frontier variable(s). " + f"Increase the budget (currently {budget}) or reduce the explanation frontier " + "size for reliable estimates.", + UserWarning, + stacklevel=3, + ) + + def _init_sv_kernel_weights(self) -> np.ndarray: + """Initialise the order-1 (Shapley-value) regression kernel weights by coalition size. + + Weights are zero for the empty and grand coalitions (handled as hard + constraints) and follow the KernelSHAP formula otherwise. + + Returns: + Weight vector of shape ``(n + 1,)``. + """ + weight_vector = np.zeros(shape=self.n + 1) + for coalition_size in range(self.n + 1): + if coalition_size < 1 or coalition_size > self.n - 1: + weight_vector[coalition_size] = 0 + else: + weight_vector[coalition_size] = 1 / ( + (self.n - 1) * binom(self.n - 2, coalition_size - 1) + ) + return weight_vector + + def _transform_to_shapley( + self, input_values: np.ndarray + ) -> tuple[np.ndarray, dict[tuple[int, ...], int]]: + """Aggregate interaction-level values into Shapley values. + + Each interaction term's value is split equally among its members, so + higher-order terms contribute a ``1/|S|`` share to each player in ``S``. + + Args: + input_values: Array of per-frontier-term values (including the + empty-coalition entry at index 0). + + Returns: + Tuple ``(sv, sv_lookup)`` where *sv* is a length-``(n + 1)`` array + of Shapley values (index 0 reserved for the empty coalition) and + *sv_lookup* maps singleton tuples to their positions in *sv*. + """ + sv = np.zeros(self.n + 1) + sv_lookup = {} + for interaction, interaction_pos in self.interaction_lookup.items(): + if len(interaction) == 0: + sv[interaction_pos] = input_values[interaction_pos] + sv_lookup[()] = interaction_pos + for i in interaction: + sv[i + 1] += input_values[interaction_pos] / len(interaction) + sv_lookup[(i,)] = i + 1 + return sv, sv_lookup + + def approximate( + self, + budget: int, + game: Game | Callable[[np.ndarray], np.ndarray], + *args: Any | None, # noqa: ARG002 + **kwargs: Any, # noqa: ARG002 + ) -> InteractionValues: + """Approximate Shapley values via weighted least-squares regression. + + Draws random coalitions, queries the game, then solves a weighted + least-squares problem whose basis functions are the terms in the explanation + frontier. When the frontier contains only singletons the method reduces exactly + to KernelSHAP :cite:t:`Lundberg.2017`. For details see Fumagalli et al. (2026) + :cite:t:`Fumagalli.2026a`. + + Args: + budget: Total number of coalition evaluations (including the empty and + grand coalition, which are always queried). + game: Callable accepting a binary coalition matrix of shape + ``(budget, n)`` and returning a value array of shape ``(budget,)``. + *args: Ignored; accepted for API compatibility with other approximators. + **kwargs: Ignored; accepted for API compatibility with other approximators. + + Returns: + :class:`~shapiq.interaction_values.InteractionValues` containing the + estimated Shapley values. + """ + kernel_weights = self._init_sv_kernel_weights() + self.projection_matrix = np.identity(self.n_variables) - 1 / self.n_variables + + # Sample coalitions and query the game. + self._sampler.sample(budget) + game_values = game(self._sampler.coalitions_matrix) + + # Centre game values on the empty-coalition baseline. + empty_set_value = game_values[0] + game_values -= empty_set_value + full_set_value = game_values[1] + + sampling_normalization = np.sqrt( + kernel_weights[self._sampler.coalitions_size[2:]] + * self._sampler.sampling_adjustment_weights[2:] + ) + + # Build the weighted design matrix. + # When interactions are included (n_variables > n) each column checks + # whether the corresponding frontier set is a subset of the sampled coalition. + # Use the actual sample count: the border-trick may cap evaluations below budget. + n_sampled = len(self._sampler.coalitions_matrix) - 2 + self._warn_if_underdefined(n_sampled, budget) + if self.n_variables > self.n: + x_tilde = np.zeros((n_sampled, self.n_variables)) + for pos, row in enumerate(self.interaction_matrix_binary[1:, :]): + x_tilde[:, pos] = ( + np.all(row <= self._sampler.coalitions_matrix[2:, :], axis=1) + * sampling_normalization + ) + else: + x_tilde = sampling_normalization[:, np.newaxis] * self._sampler.coalitions_matrix[2:, :] + + y_tilde = game_values[2:] * sampling_normalization + + # Solve the weighted least-squares problem. + least_squares_solution = np.linalg.lstsq( + x_tilde @ self.projection_matrix, + y_tilde - full_set_value / self.n_variables * np.sum(x_tilde, axis=1), + rcond=None, + )[0] + + interaction_representation = np.zeros(self.n_variables + 1, dtype=float) + interaction_representation[0] = empty_set_value + interaction_representation[1:] = full_set_value / self.n_variables + least_squares_solution + + sv, sv_lookup = self._transform_to_shapley(interaction_representation) + + return InteractionValues( + values=sv, + index="SV", + interaction_lookup=sv_lookup, + baseline_value=float(empty_set_value), + min_order=self.min_order, + max_order=self.max_order, + n_players=self.n, + estimated=not budget >= 2**self.n, + estimation_budget=budget, + target_index=self.index, + ) diff --git a/tests/shapiq/tests_unit/tests_approximators/test_approximator_polyshap.py b/tests/shapiq/tests_unit/tests_approximators/test_approximator_polyshap.py new file mode 100644 index 000000000..9243bf6d3 --- /dev/null +++ b/tests/shapiq/tests_unit/tests_approximators/test_approximator_polyshap.py @@ -0,0 +1,595 @@ +"""Tests for the :class:`~shapiq.approximator.PolySHAP` approximator. + +PolySHAP selects its interaction frontier through one of three constructor arguments: + +* **k-additive** (default): ``max_order`` (with ``max_terms=None``). +* **partial**: ``max_terms`` (budget-capped frontier, whole low orders first). +* **prior**: ``prior_frontier`` supplies the interaction terms directly. + +The suite exercises every branch of the frontier builder and both design-matrix +paths of :meth:`~shapiq.approximator.PolySHAP.approximate` (the order-1 KernelSHAP path +and the higher-order interaction path), checks convergence to the ground truth from +:class:`~shapiq.ExactComputer`, and runs reproducibility / variance sweeps over a +fixed list of 50 seeds. +""" + +from __future__ import annotations + +from functools import cache +from itertools import pairwise + +import numpy as np +import pytest +from scipy.special import binom + +import shapiq.approximator as approximator_module +from shapiq import ExactComputer +from shapiq.approximator import PolySHAP +from shapiq.interaction_values import InteractionValues +from shapiq.utils.sets import powerset +from shapiq_games.synthetic import DummyGame + +# The sampler's border-trick and the underdefined-system check both emit UserWarnings +# that are irrelevant to most assertions here; individual tests opt back in where needed. +pytestmark = pytest.mark.filterwarnings("ignore::UserWarning") + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + + +def _singleton_prior(n: int) -> list[tuple[int, ...]]: + """Minimal prior: empty set + all singletons (the KernelSHAP frontier).""" + return list(powerset(range(n), max_size=1)) + + +def _kadd_frontier_size(n: int, max_order: int) -> int: + """Number of subsets of ``range(n)`` of size ``0 .. max_order``.""" + return int(sum(binom(n, k) for k in range(max_order + 1))) + + +@cache +def _exact_sv(n: int, interaction: tuple[int, ...]) -> tuple[float, ...]: + """Ground-truth Shapley values for ``DummyGame(n, interaction)`` via ExactComputer.""" + game = DummyGame(n, interaction) + values = ExactComputer(game, n_players=n)(index="SV", order=1).get_n_order_values(1) + return tuple(float(v) for v in np.asarray(values)) + + +def _estimate_vector(approx: PolySHAP, budget: int, game: DummyGame) -> np.ndarray: + """Run ``approx`` and return its per-player Shapley-value vector of length ``n``.""" + sv = approx.approximate(budget, game) + return np.array([sv[(i,)] for i in range(approx.n)]) + + +def _generate_seeds(rng_seed: int = 20260617, count: int = 50) -> list[int]: + """Generate a fixed, reproducible list of seeds spanning a wide value range.""" + return [int(s) for s in np.random.default_rng(rng_seed).integers(0, 2**31 - 1, size=count)] + + +SEEDS: list[int] = _generate_seeds() + +# Shared configuration for the seeded sweeps. +_SEEDED_N = 12 +_SEEDED_INTERACTION = (1, 2) +_SEEDED_MAX_TERMS = 30 +_SEEDED_BUDGET = 400 + +# The three modes, as (label, constructor-kwargs) pairs, reused by shared tests. +_MODES: list[tuple[str, dict]] = [ + ("kadd_order1", {"max_order": 1}), + ("kadd_order2", {"max_order": 2}), + ("partial", {"max_terms": 12}), + ("prior", {"prior_frontier": _singleton_prior(7)}), +] +_MODE_IDS = [label for label, _ in _MODES] +_MODE_KWARGS = [kwargs for _, kwargs in _MODES] + + +# --------------------------------------------------------------------------- +# Registration / exports +# --------------------------------------------------------------------------- + + +def test_polyshap_is_exported_and_registered(): + """PolySHAP must be publicly exported and registered as an SV approximator.""" + assert "PolySHAP" in approximator_module.__all__ + assert PolySHAP in approximator_module.SV_APPROXIMATORS + + +# --------------------------------------------------------------------------- +# Frontier construction: k-additive mode +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize(("n", "max_order"), [(3, 1), (3, 2), (7, 1), (7, 2), (10, 3)]) +def test_kadd_frontier_size_and_attributes(n, max_order): + """The k-additive frontier holds every subset up to ``max_order`` and sets attributes.""" + approx = PolySHAP(n, max_order=max_order) + expected_size = _kadd_frontier_size(n, max_order) + assert approx.n == n + assert approx.max_order == 1 # the *output* order: PolySHAP always targets SVs + assert approx.min_order == 0 + assert approx.top_order is False + assert approx.iteration_cost == 1 + assert approx.index == "SV" + assert len(approx.explanation_frontier) == expected_size + assert approx.n_variables == expected_size - 1 + + +def test_default_constructor_is_two_additive(): + """``PolySHAP(n)`` defaults to the 2-additive frontier.""" + n = 6 + approx = PolySHAP(n) + assert len(approx.explanation_frontier) == _kadd_frontier_size(n, 2) + + +def test_kadd_order1_frontier_is_kernelshap_frontier(): + """``max_order=1`` yields exactly the empty set plus all singletons.""" + n = 6 + frontier = PolySHAP(n, max_order=1).explanation_frontier + assert len(frontier) == n + 1 + assert max(len(S) for S in frontier) == 1 + + +@pytest.mark.parametrize("n", [4, 7]) +def test_kadd_frontier_contains_empty_set_and_all_singletons(n): + """The empty set sits at index 0 and every singleton is present.""" + frontier = PolySHAP(n, max_order=2).explanation_frontier + assert frontier[()] == 0 + for i in range(n): + assert (i,) in frontier + + +@pytest.mark.parametrize("excluded_size", [2, 3]) +def test_kadd_sizes_to_exclude_removes_those_sizes(excluded_size): + """Excluded higher-order sizes are absent, while singletons are always kept.""" + n = 6 + approx = PolySHAP(n, max_order=4, sizes_to_exclude={excluded_size}) + for coalition in approx.explanation_frontier: + assert len(coalition) != excluded_size + for i in range(n): + assert (i,) in approx.explanation_frontier + + +def test_kadd_sizes_to_exclude_singletons_is_ignored(): + """Singletons are mandatory, so excluding size 1 has no effect (no error).""" + n = 4 + approx = PolySHAP(n, max_order=2, sizes_to_exclude={1}) + for i in range(n): + assert (i,) in approx.explanation_frontier + + +def test_kadd_frontier_positions_are_unique(): + """Each frontier term maps to a unique, contiguous column position.""" + approx = PolySHAP(5, max_order=2) + positions = sorted(approx.explanation_frontier.values()) + assert positions == list(range(len(positions))) + + +def test_interaction_matrix_binary_shape(): + """``interaction_matrix_binary`` has shape ``(|frontier|, n)``.""" + n = 4 + approx = PolySHAP(n, max_order=2) + assert approx.interaction_matrix_binary.shape == (len(approx.explanation_frontier), n) + + +# --------------------------------------------------------------------------- +# Frontier construction: partial mode +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + ("n", "max_terms"), + [ + (7, 8), # empty + singletons only, no higher-order extension + (7, 12), # empty + singletons + 4 pairs + (7, 25), # empty + singletons + 17 pairs (21 pairs exist -> fits) + ], +) +def test_partial_frontier_size_when_within_capacity(n, max_terms): + """When capacity allows, the partial frontier holds exactly ``max_terms`` entries.""" + approx = PolySHAP(n, max_terms=max_terms, random_state=42) + assert len(approx.explanation_frontier) == max_terms + assert approx.n_variables == max_terms - 1 + + +def test_partial_frontier_capped_by_capacity_when_candidates_exhausted(): + """If ``max_terms`` exceeds the available terms up to ``max_order`` the frontier caps out.""" + n = 5 # capacity at order 2 is 1 + 5 + C(5, 2) = 16 + approx = PolySHAP(n, max_order=2, max_terms=100, random_state=0) + assert len(approx.explanation_frontier) == _kadd_frontier_size(n, 2) + assert max(len(S) for S in approx.explanation_frontier) == 2 + + +def test_partial_is_bounded_by_max_order(): + """Partial extension never adds terms above ``max_order``.""" + order2 = PolySHAP(6, max_order=2, max_terms=40, random_state=0) + assert max(len(S) for S in order2.explanation_frontier) == 2 + order3 = PolySHAP(6, max_order=3, max_terms=40, random_state=0) + assert max(len(S) for S in order3.explanation_frontier) == 3 + + +def test_partial_always_contains_all_singletons(): + """Every singleton appears regardless of the cap.""" + n = 7 + approx = PolySHAP(n, max_order=3, max_terms=15, random_state=1) + for i in range(n): + assert (i,) in approx.explanation_frontier + + +def test_partial_sizes_to_exclude_omits_those_sizes(): + """Excluded sizes never appear among the randomly selected terms.""" + approx = PolySHAP(7, max_order=3, max_terms=25, sizes_to_exclude={2}, random_state=0) + for coalition in approx.explanation_frontier: + assert len(coalition) != 2 + + +def test_partial_fills_whole_low_orders_before_sampling_boundary_order(): + """Partial mode realizes the paper's ``I_ell``: complete low orders, partial top order. + + With ``n=6``, ``max_order=3`` and ``max_terms=30`` the frontier holds the empty set, + 6 singletons and all 15 pairs (22 terms), then a *random* subset of the 20 triples to + reach 30. Every pair must be present while only some triples are. + """ + n, max_terms = 6, 30 + frontier = PolySHAP(n, max_order=3, max_terms=max_terms, random_state=0).explanation_frontier + all_pairs = list(powerset(range(n), min_size=2, max_size=2)) + triples = [S for S in frontier if len(S) == 3] + + assert all(pair in frontier for pair in all_pairs) # every lower-order term kept + assert 0 < len(triples) < int(binom(n, 3)) # boundary order only partially included + assert len(frontier) == max_terms + + +def test_partial_different_seeds_produce_different_frontiers(): + """Different random states select different higher-order interactions.""" + n, max_terms = 7, 20 + a1 = PolySHAP(n, max_order=3, max_terms=max_terms, random_state=1) + a2 = PolySHAP(n, max_order=3, max_terms=max_terms, random_state=2) + assert set(a1.explanation_frontier) != set(a2.explanation_frontier) + + +def test_partial_same_seed_reproduces_frontier(): + """The same random state yields an identical frontier.""" + n, max_terms = 7, 20 + a1 = PolySHAP(n, max_order=3, max_terms=max_terms, random_state=7) + a2 = PolySHAP(n, max_order=3, max_terms=max_terms, random_state=7) + assert set(a1.explanation_frontier) == set(a2.explanation_frontier) + + +# --------------------------------------------------------------------------- +# Frontier construction: prior mode +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("n", [3, 7, 10]) +def test_prior_frontier_matches_supplied_ordering(n): + """The prior frontier equals ``prior_frontier`` with positions in enumeration order.""" + prior = [*_singleton_prior(n), (0, 1), (1, 2)] + approx = PolySHAP(n, prior_frontier=prior) + assert len(approx.explanation_frontier) == len(prior) + assert approx.n_variables == len(prior) - 1 + for pos, coalition in enumerate(prior): + assert approx.explanation_frontier[coalition] == pos + + +def test_informed_prior_recovers_exact_interaction_game(): + """A prior carrying the game's true interaction recovers the exact Shapley values. + + This is the realistic use of *prior* mode: on a 2-additive game whose only + interaction is the pair ``(1, 2)``, supplying that pair (plus the singletons) as the + frontier lets PolySHAP recover the exact SVs from a subsampled budget (``n=12``, + ``budget << 2**n``), whereas a singleton-only prior (KernelSHAP) cannot. + """ + n, budget, interaction = 12, 200, (1, 2) + exact = np.asarray(_exact_sv(n, interaction)) + + informed = [*_singleton_prior(n), interaction] + est = _estimate_vector( + PolySHAP(n, prior_frontier=informed, random_state=1), budget, DummyGame(n, interaction) + ) + np.testing.assert_allclose(est, exact, atol=1e-6) + + # A singleton-only prior lacks the interaction term and leaves a visible error. + est_singletons = _estimate_vector( + PolySHAP(n, prior_frontier=_singleton_prior(n), random_state=1), + budget, + DummyGame(n, interaction), + ) + assert np.sqrt(np.mean((est_singletons - exact) ** 2)) > 1e-3 + + +# --------------------------------------------------------------------------- +# Frontier construction: validation / error paths +# --------------------------------------------------------------------------- + + +def test_prior_frontier_with_max_order_raises(): + """``prior_frontier`` combined with ``max_order`` is rejected.""" + with pytest.raises(ValueError, match="mutually exclusive"): + PolySHAP(5, prior_frontier=_singleton_prior(5), max_order=2) + + +def test_prior_frontier_with_max_terms_raises(): + """``prior_frontier`` combined with ``max_terms`` is rejected.""" + with pytest.raises(ValueError, match="mutually exclusive"): + PolySHAP(5, prior_frontier=_singleton_prior(5), max_terms=10) + + +def test_prior_frontier_with_sizes_to_exclude_raises(): + """``prior_frontier`` combined with ``sizes_to_exclude`` is rejected.""" + with pytest.raises(ValueError, match="mutually exclusive"): + PolySHAP(5, prior_frontier=_singleton_prior(5), sizes_to_exclude={2}) + + +def test_max_terms_below_minimum_raises(): + """``max_terms`` must leave room for the empty set and all singletons.""" + n = 6 + with pytest.raises(ValueError, match="at least"): + PolySHAP(n, max_terms=n) # needs >= n + 1 + + +def test_prior_missing_singleton_raises(): + """A prior that omits a singleton is rejected by the base-class check.""" + n = 5 + incomplete = [(), (0,), (1,), (2,)] # players 3 and 4 absent + with pytest.raises(ValueError, match="main effects"): + PolySHAP(n, prior_frontier=incomplete) + + +# --------------------------------------------------------------------------- +# Approximation output contract +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("kwargs", _MODE_KWARGS, ids=_MODE_IDS) +def test_output_is_valid_interaction_values(kwargs): + """Every mode returns a well-formed SV ``InteractionValues`` object.""" + n, budget = 7, 200 + game = DummyGame(n, (1, 2)) + sv = PolySHAP(n, **kwargs, random_state=0).approximate(budget, game) + assert isinstance(sv, InteractionValues) + assert sv.index == "SV" + assert sv.max_order == 1 + assert sv.min_order == 0 + assert sv.estimation_budget <= budget + for i in range(n): + assert (i,) in sv.interaction_lookup + + +@pytest.mark.parametrize("kwargs", _MODE_KWARGS, ids=_MODE_IDS) +def test_budget_never_exceeded(kwargs): + """The game is never queried more than ``budget`` times.""" + n, budget = 7, 150 + game = DummyGame(n, (1, 2)) + PolySHAP(n, **kwargs, random_state=0).approximate(budget, game) + assert game.access_counter <= budget + + +@pytest.mark.parametrize("kwargs", _MODE_KWARGS, ids=_MODE_IDS) +def test_sv_efficiency_holds(kwargs): + """The Shapley values sum to ``v(N) - v(empty) = 2.0`` for ``DummyGame(7, (1, 2))``.""" + n = 7 + game = DummyGame(n, (1, 2)) + sv = PolySHAP(n, **kwargs, random_state=0).approximate(500, game) + assert np.sum(sv.values) == pytest.approx(2.0, abs=0.05) + + +def test_pairing_trick_produces_valid_output(): + """``pairing_trick=True`` still yields accurate estimates (order-1 path).""" + n, budget = 7, 200 + game = DummyGame(n, (1, 2)) + est = _estimate_vector( + PolySHAP(n, max_order=1, pairing_trick=True, random_state=7), budget, game + ) + exact = np.asarray(_exact_sv(n, (1, 2))) + np.testing.assert_allclose(est, exact, atol=0.15) + + +# --------------------------------------------------------------------------- +# Convergence to the ExactComputer ground truth +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("kwargs", _MODE_KWARGS, ids=_MODE_IDS) +def test_converges_to_exact_sv_with_sufficient_budget(kwargs): + """With ample budget every mode recovers the exact Shapley values (n <= 10).""" + n, budget = 7, 600 + game = DummyGame(n, (1, 2)) + est = _estimate_vector(PolySHAP(n, **kwargs, random_state=42), budget, game) + exact = np.asarray(_exact_sv(n, (1, 2))) + np.testing.assert_allclose(est, exact, atol=0.1) + + +def test_matching_frontier_order_recovers_exact_two_additive_game(): + """A 2-additive game is recovered exactly by ``max_order=2`` but not ``max_order=1``. + + Uses ``n = 12`` with ``budget << 2**n`` so the estimate is genuinely subsampled + rather than a full-space enumeration. + """ + n, budget = 12, 200 + interaction = (1, 2) + exact = np.asarray(_exact_sv(n, interaction)) + + est_order2 = _estimate_vector( + PolySHAP(n, max_order=2, random_state=1), budget, DummyGame(n, interaction) + ) + np.testing.assert_allclose(est_order2, exact, atol=1e-6) + + est_order1 = _estimate_vector( + PolySHAP(n, max_order=1, random_state=1), budget, DummyGame(n, interaction) + ) + assert np.sqrt(np.mean((est_order1 - exact) ** 2)) > 1e-3 + + +def test_higher_order_reduces_error_on_three_additive_game(): + """On a 3-way interaction game, ``max_order=3`` recovers exact SVs while order 2 lags.""" + n, budget = 12, 300 + interaction = (1, 2, 3) + exact = np.asarray(_exact_sv(n, interaction)) + + rmse = {} + for order in (2, 3): + est = _estimate_vector( + PolySHAP(n, max_order=order, random_state=1), budget, DummyGame(n, interaction) + ) + rmse[order] = float(np.sqrt(np.mean((est - exact) ** 2))) + + np.testing.assert_allclose( + _estimate_vector( + PolySHAP(n, max_order=3, random_state=1), budget, DummyGame(n, interaction) + ), + exact, + atol=1e-6, + ) + assert rmse[3] < rmse[2] + + +# --------------------------------------------------------------------------- +# Edge cases +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("n", [2, 3]) +def test_tiny_games(n): + """PolySHAP works for very small player counts and stays efficient.""" + game = DummyGame(n, (0, 1)) + sv = PolySHAP(n, max_order=1, random_state=0).approximate(2**n, game) + assert np.sum(sv.values) == pytest.approx(2.0, abs=1e-6) + + +def test_budget_smaller_than_frontier_still_runs(): + """A budget far below the frontier size runs and returns a full SV vector.""" + n = 8 + game = DummyGame(n, (1, 2)) + sv = PolySHAP(n, max_order=2, random_state=0).approximate(6, game) + assert isinstance(sv, InteractionValues) + for i in range(n): + assert (i,) in sv.interaction_lookup + + +def test_underdefined_system_warns(): + """An under-sampled least-squares system emits an explanatory UserWarning.""" + n = 8 # order-2 frontier has 36 variables; budget of 10 is far too small + game = DummyGame(n, (1, 2)) + with pytest.warns(UserWarning, match="underdefined"): + PolySHAP(n, max_order=2, random_state=0).approximate(10, game) + + +def test_estimated_flag_true_when_budget_small(): + """``estimated`` is True when the budget is below the full coalition space.""" + n = 7 + sv = PolySHAP(n, max_order=1, random_state=0).approximate(50, DummyGame(n, (1, 2))) + assert sv.estimated is True + + +def test_estimated_flag_false_when_full_space_covered(): + """``estimated`` is False when the budget exceeds ``2**n``.""" + n = 5 + sv = PolySHAP(n, max_order=1, random_state=0).approximate(2**n + 10, DummyGame(n, (1, 2))) + assert sv.estimated is False + + +# --------------------------------------------------------------------------- +# Reproducibility under a fixed random_state +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("kwargs", _MODE_KWARGS, ids=_MODE_IDS) +def test_random_state_reproducibility(kwargs): + """Identical ``random_state`` yields bit-identical estimates for every mode.""" + n, budget = 7, 200 + sv1 = PolySHAP(n, **kwargs, random_state=99).approximate(budget, DummyGame(n, (1, 2))) + sv2 = PolySHAP(n, **kwargs, random_state=99).approximate(budget, DummyGame(n, (1, 2))) + np.testing.assert_array_equal(sv1.values, sv2.values) + + +def test_prior_singleton_matches_kadd_order1(): + """A singleton prior and ``max_order=1`` share the same frontier and results.""" + n, budget, seed = 7, 200, 5 + sv_prior = PolySHAP(n, prior_frontier=_singleton_prior(n), random_state=seed).approximate( + budget, DummyGame(n, (1, 2)) + ) + sv_kadd = PolySHAP(n, max_order=1, random_state=seed).approximate(budget, DummyGame(n, (1, 2))) + np.testing.assert_array_almost_equal(sv_prior.values, sv_kadd.values, decimal=12) + + +# --------------------------------------------------------------------------- +# Seeded sweeps over 50 fixed seeds +# --------------------------------------------------------------------------- + + +def test_seed_list_is_deterministic_and_unique(): + """The fixed seed list is regenerable and contains 50 unique seeds.""" + assert _generate_seeds() == SEEDS + assert len(SEEDS) == 50 + assert len(set(SEEDS)) == 50 + + +def test_seeded_runs_are_reproducible(): + """For every seed, two independent partial-mode runs are bit-identical.""" + for seed in SEEDS: + est1 = _estimate_vector( + PolySHAP(_SEEDED_N, max_terms=_SEEDED_MAX_TERMS, random_state=seed), + _SEEDED_BUDGET, + DummyGame(_SEEDED_N, _SEEDED_INTERACTION), + ) + est2 = _estimate_vector( + PolySHAP(_SEEDED_N, max_terms=_SEEDED_MAX_TERMS, random_state=seed), + _SEEDED_BUDGET, + DummyGame(_SEEDED_N, _SEEDED_INTERACTION), + ) + np.testing.assert_array_equal( + est1, est2, err_msg=f"non-reproducible result for seed {seed}" + ) + + +def test_seeded_estimates_are_unbiased_with_low_variance(): + """Across 50 seeds the partial-mode estimates cluster tightly around the true SVs.""" + exact = np.asarray(_exact_sv(_SEEDED_N, _SEEDED_INTERACTION)) + + estimates = np.empty((len(SEEDS), _SEEDED_N)) + for row, seed in enumerate(SEEDS): + estimates[row] = _estimate_vector( + PolySHAP(_SEEDED_N, max_terms=_SEEDED_MAX_TERMS, random_state=seed), + _SEEDED_BUDGET, + DummyGame(_SEEDED_N, _SEEDED_INTERACTION), + ) + + # Near-unbiased: the mean over seeds matches the exact Shapley values. + np.testing.assert_allclose(estimates.mean(axis=0), exact, atol=0.02) + # Tight spread around the truth. + assert estimates.std(axis=0).max() < 0.05 + # Efficiency holds on average. + assert estimates.sum(axis=1).mean() == pytest.approx(2.0, abs=0.01) + + +def test_error_decreases_as_budget_grows(): + """Mean RMSE over 50 seeds drops monotonically as the budget increases (order-1).""" + n, interaction = 8, (1, 2) + exact = np.asarray(_exact_sv(n, interaction)) + budgets = [40, 80, 160, 240] # all below 2**8 = 256, so genuinely estimated + + mean_rmse = [] + for budget in budgets: + rmses = [ + np.sqrt( + np.mean( + ( + _estimate_vector( + PolySHAP(n, max_order=1, random_state=seed), + budget, + DummyGame(n, interaction), + ) + - exact + ) + ** 2 + ) + ) + for seed in SEEDS + ] + mean_rmse.append(float(np.mean(rmses))) + + assert all(later < earlier for earlier, later in pairwise(mean_rmse))