From 163807b7a0afbc5663824cd1fbb4c0f7fd6bc456 Mon Sep 17 00:00:00 2001 From: YUANLUO WU Date: Wed, 10 Jun 2026 10:23:54 +0200 Subject: [PATCH 1/3] fix(tree): guard the ill-conditioned Chebyshev-Vandermonde solves in TreeSHAP-IQ The polynomial TreeSHAP machinery computed inv(vander(points).T) @ rhs at four call sites (TreeSHAPIQ N/N_cii/N_id matrices and LinearTreeSHAP's N_v2). The interpolation nodes are the first i entries of a depth-sized Chebyshev grid, so the Vandermonde systems are ill-conditioned even for moderate sizes: precision degrades silently from roughly size 30 (residuals > 1e-6 around 29, > 1e-3 around 37, O(1) around 45 on a depth-grid basis) and the matrix becomes exactly singular for very deep trees, which crashed with an unexplained LinAlgError (observed for trees of depth ~55-60). This change centralises the solve in tree/_numerics.solve_vandermonde(): - np.linalg.solve instead of forming an explicit inverse, - a RuntimeWarning when the condition number exceeds 1e12 (precision loss), - a least-squares fallback instead of a hard crash when the system is singular, with a warning that values may be inaccurate. Deep trees now explain (approximately) instead of crashing, and users are told when the exactness contract can no longer be honoured. Minimal reproduction: fit DecisionTreeRegressor(max_depth=60) on noise and construct LinearTreeSHAP(tree) - previously LinAlgError, now a warning. --- CHANGELOG.md | 14 ++ src/shapiq/tree/_numerics.py | 207 ++++++++++++++++++ src/shapiq/tree/linear/explainer.py | 9 +- src/shapiq/tree/treeshapiq.py | 27 ++- .../test_tree_numerics.py | 94 ++++++++ 5 files changed, 346 insertions(+), 5 deletions(-) create mode 100644 src/shapiq/tree/_numerics.py create mode 100644 tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 76cd143d2..bd262a639 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,19 @@ # Changelog +## Unreleased + +### Bugfix + +- Hardens the Vandermonde interpolation solves in the polynomial TreeSHAP machinery (`TreeSHAPIQ`, `LinearTreeSHAP`): + - **surfaces a pre-existing silent-wrong-answer bug**: for interpolation degrees ~27-31 (for `TreeSHAPIQ` the degree is `min(depth, n_features)`), the previous explicit matrix inversion silently lost precision - a coalesced `RuntimeWarning` now reports the worst condition number; from degree ~32 the prefix systems become rank-deficient and the returned values can violate the completeness axiom by errors up to orders of magnitude beyond the prediction itself - such cases now emit a clear `RuntimeWarning` stating the values are NOT reliable, + - very deep trees no longer crash with an unexplained `LinAlgError` (previously at degree ~60); a least-squares fallback is returned together with the reliability warning, + - no explicit matrix inverse is formed anymore (`np.linalg.solve` / SVD-based least squares), + - interpolation grids are certified by measurement (cached per grid, hoisted out of the construction loops), so custom `LinearTreeSHAP(base_func=...)` grids are handled correctly and ordinary trees pay no per-solve diagnostic work. [#545](https://github.com/mmschlk/shapiq/issues/545) + +### Changed + +- `LinearTreeSHAP` (trees deeper than ~26) and `TreeSHAPIQ` (interpolation degree above ~26) now emit a `RuntimeWarning` on paths that previously were silent (precision loss / unreliable deep-tree values); test suites running with `filterwarnings = error` will surface this intentionally. Likely warrants a minor (not patch) version bump. + ## v1.5.1 (2026-05-30) ### Bugfix diff --git a/src/shapiq/tree/_numerics.py b/src/shapiq/tree/_numerics.py new file mode 100644 index 000000000..edb3c8f86 --- /dev/null +++ b/src/shapiq/tree/_numerics.py @@ -0,0 +1,207 @@ +"""Numerically guarded solves for the Vandermonde systems of TreeSHAP-IQ. + +The polynomial TreeSHAP machinery repeatedly needs ``inv(vander(points).T) @ rhs``, +where ``points`` is a prefix ``D[:i]`` of a fixed interpolation grid ``D`` whose size +equals the interpolation degree (``min(tree depth, n_features_in_tree)`` for +TreeSHAPIQ, the full tree depth for LinearTreeSHAP). These Vandermonde systems are +severely ill-conditioned, and the conditioning is **non-monotonic in the prefix +length**: the first ``i`` points of an ``n``-point grid typically cluster near one +endpoint, so the condition number peaks around ``i ~ n/2``. For the default +Chebyshev grids, measured peak prefix condition numbers are ``4.1e11`` at grid size +26, ``1.3e12`` at 27 and ``4.2e13`` at 30 — precision degrades for grids larger +than ~26 (at *interior* prefixes, not only at the full size), and the systems become +singular to machine precision for very deep trees (sizes around 55-60), which +previously surfaced as an unexplained ``LinAlgError``. + +This module centralises the solve: + +- **Grid certification instead of assumptions.** For each distinct grid, the peak + condition number over all prefixes is measured once and cached + (:func:`grid_is_certified`); call sites hoist this decision out of their + construction loops, so certified grids take a plain ``np.linalg.solve`` with no + per-solve diagnostic work. The certification is computed on the actual grid, so + custom interpolation bases (``LinearTreeSHAP(base_func=...)``) are handled + correctly rather than assumed to behave like Chebyshev nodes. +- **Checked path with one SVD.** Uncertified grids are solved by SVD-based least + squares; the singular values of that same decomposition give the exact condition + number of the matrix actually solved, and rank-deficient systems degrade into a + least-squares solution instead of crashing. +- **One warning per matrix, not per prefix.** Call sites construct whole N-matrices + in a loop over prefixes; per-prefix warnings would flood the user (a depth-60 + construction touches dozens of degenerate prefixes). The loop passes a + :class:`VandermondeDiagnostics` collector and emits a single summary warning. + +Design note: rank-deficient systems deliberately **warn and return** the +least-squares fallback rather than raising. Before this module, the same depths +either crashed (``LinAlgError``, degree ~60) or — worse — **silently returned +wrong values** (rank-deficient prefixes appear from grid size ~32; completeness-axiom errors up to +orders of magnitude beyond the prediction itself). Warning loudly while degrading is strictly +more informative than either previous behaviour, and pipelines that want hard +failures can promote the warning with the ``warnings`` module. +""" + +from __future__ import annotations + +import warnings +from functools import lru_cache + +import numpy as np + +# Above this condition number, double precision has lost ~12 of its ~16 significant +# digits; results may no longer be exact. (Empirical anchor for the default +# Chebyshev grids: the peak prefix condition number stays below this threshold up +# to grid size 26 and exceeds it from size 27 on; see module docstring. The +# certification below measures each actual grid, so the constant is a warning +# threshold, not a structural assumption.) +_COND_WARN_THRESHOLD = 1e12 + +# Certification is only attempted for grids up to this size; larger grids go +# straight to the checked path (their peak prefix conditioning exceeds the +# threshold for every node family used in practice, and bounding the cache work +# keeps certification O(small)). +_CERTIFY_MAX_SIZE = 32 + + +@lru_cache(maxsize=128) +def _max_prefix_condition(grid_key: tuple[float, ...]) -> float: + """Measured peak condition number over all prefixes of the given grid.""" + grid = np.asarray(grid_key) + worst = 1.0 + for i in range(2, len(grid) + 1): + try: + worst = max(worst, float(np.linalg.cond(np.vander(grid[:i]).T))) + except np.linalg.LinAlgError: # pragma: no cover - cond rarely raises + return float("inf") + return worst + + +class VandermondeDiagnostics: + """Collects conditioning diagnostics across the solves of one N-matrix. + + Use one instance per matrix-construction loop and call :meth:`emit` once after + the loop, so the user sees a single summary warning instead of one warning per + degenerate prefix. + """ + + def __init__(self) -> None: + self.worst_condition: float = 0.0 + self.rank_deficient_sizes: list[int] = [] + self.worst_residual: float = 0.0 + + def emit(self, stacklevel: int = 4) -> None: + """Emit at most one summary warning for the collected diagnostics. + + Args: + stacklevel: Forwarded to :func:`warnings.warn`; pass the depth of the + user's call site relative to ``emit`` (the construction loops sit + at different depths in LinearTreeSHAP vs TreeSHAPIQ). + """ + if self.rank_deficient_sizes: + warnings.warn( + f"Interpolation systems of size {self.rank_deficient_sizes} are " + "singular to machine precision (very deep tree). Least-squares " + "fallbacks were used (worst relative residual " + f"{self.worst_residual:.1e}) - the resulting explanation values are " + "NOT reliable. Limit the tree depth to obtain exact values.", + RuntimeWarning, + stacklevel=stacklevel, + ) + elif self.worst_condition > _COND_WARN_THRESHOLD: + warnings.warn( + "The interpolation systems for this tree reach a condition number " + f"of {self.worst_condition:.1e}; TreeSHAP-IQ values may lose " + "precision. Consider limiting the tree depth (the interpolation " + "degree is min(tree depth, number of features in the tree)).", + RuntimeWarning, + stacklevel=stacklevel, + ) + + +def grid_is_certified(grid: np.ndarray) -> bool: + """Whether every prefix system of ``grid`` is measured well-conditioned. + + Note that "certified" bounds, but does not eliminate, precision loss: the + threshold admits condition numbers up to 1e12 (a size-26 Chebyshev grid peaks + at ~4e11). Certification costs one SVD per prefix on first touch of a grid + (cached afterwards, and shared across the N-matrix builders), roughly doubling + the linear-algebra work of the first construction for that grid size. + + Call once per interpolation grid (the measurement is cached per grid) and pass + the result to :func:`solve_vandermonde` as ``certified=`` so the per-prefix + solves carry no diagnostic work at all for ordinary trees. + """ + grid = np.asarray(grid) + if len(grid) > _CERTIFY_MAX_SIZE: + return False + return _max_prefix_condition(tuple(grid.tolist())) <= _COND_WARN_THRESHOLD + + +def solve_vandermonde( + points: np.ndarray, + rhs: np.ndarray, + *, + certified: bool = False, + diagnostics: VandermondeDiagnostics | None = None, +) -> np.ndarray: + """Solve ``vander(points).T @ x = rhs`` without forming an explicit inverse. + + Args: + points: Interpolation nodes (a prefix of the interpolation grid). + rhs: Right-hand side vector of the same length. + certified: Result of :func:`grid_is_certified` for the *full* interpolation + grid that ``points`` is a prefix of (conditioning peaks at interior + prefixes, so safety is a property of the whole grid). When ``True``, a + plain solve is used with no diagnostic work. Defaults to ``False`` + (checked SVD-based path). + diagnostics: Optional collector. If given, conditioning findings are + recorded on it (call :meth:`VandermondeDiagnostics.emit` after the + construction loop) instead of warning per call. + + Returns: + The solution vector ``x``. + + Warns: + RuntimeWarning: Only when ``diagnostics`` is ``None``: if the solved + system's condition number exceeds ``1e12`` (possible precision loss), + or if it is rank-deficient and a least-squares fallback is returned + whose values are not reliable. + """ + V = np.vander(points).T + size = len(points) + if certified: + # Certified grid: every prefix system measured well-conditioned. + return np.linalg.solve(V, rhs) + + # Checked path: one SVD yields the solution, the exact condition number of this + # very matrix, and rank-deficiency detection. + solution, _, rank, singular_values = np.linalg.lstsq(V, rhs, rcond=None) + if rank < size: + residual = float(np.linalg.norm(V @ solution - rhs)) + scale = float(np.linalg.norm(rhs)) or 1.0 + if diagnostics is not None: + diagnostics.rank_deficient_sizes.append(size) + diagnostics.worst_residual = max(diagnostics.worst_residual, residual / scale) + else: + warnings.warn( + f"The interpolation system of size {size} is singular to machine " + "precision (very deep tree). A least-squares fallback is returned " + f"with relative residual {residual / scale:.1e} - the resulting " + "explanation values are NOT reliable. Limit the tree depth to " + "obtain exact values.", + RuntimeWarning, + stacklevel=2, + ) + return solution + cond = float(singular_values[0] / singular_values[-1]) + if diagnostics is not None: + diagnostics.worst_condition = max(diagnostics.worst_condition, cond) + elif cond > _COND_WARN_THRESHOLD: + warnings.warn( + f"The interpolation system of size {size} has condition number " + f"{cond:.1e}; TreeSHAP-IQ values for trees this deep may lose " + "precision. Consider limiting the tree depth (the interpolation " + "degree is min(tree depth, number of features in the tree)).", + RuntimeWarning, + stacklevel=2, + ) + return solution diff --git a/src/shapiq/tree/linear/explainer.py b/src/shapiq/tree/linear/explainer.py index 2c565dc77..9c469ba14 100644 --- a/src/shapiq/tree/linear/explainer.py +++ b/src/shapiq/tree/linear/explainer.py @@ -8,6 +8,7 @@ import scipy.special as sp from shapiq.interaction_values import InteractionValues +from shapiq.tree._numerics import VandermondeDiagnostics, grid_is_certified, solve_vandermonde from shapiq.tree.conversion.edges import create_edge_tree from shapiq.tree.validation import validate_tree_model from shapiq.utils.sets import generate_interaction_lookup, powerset @@ -38,8 +39,14 @@ def get_N_v2(D: np.ndarray) -> np.ndarray: """Get N_v2 matrix for Linear Tree Shap.""" depth = D.shape[0] Ns = np.zeros((depth + 1, depth)) + diagnostics = VandermondeDiagnostics() + certified = grid_is_certified(D) for i in range(1, depth + 1): - Ns[i, :i] = np.linalg.inv(np.vander(D[:i]).T).dot(1.0 / get_norm_weight(i - 1)) + Ns[i, :i] = solve_vandermonde( + D[:i], 1.0 / get_norm_weight(i - 1), certified=certified, diagnostics=diagnostics + ) + # default stacklevel traces emit <- get_N_v2 <- __init__ <- user + diagnostics.emit() return Ns diff --git a/src/shapiq/tree/treeshapiq.py b/src/shapiq/tree/treeshapiq.py index 450e9fe84..defe6b055 100644 --- a/src/shapiq/tree/treeshapiq.py +++ b/src/shapiq/tree/treeshapiq.py @@ -13,6 +13,7 @@ from shapiq.interaction_values import InteractionValues from shapiq.utils.sets import generate_interaction_lookup, powerset +from ._numerics import VandermondeDiagnostics, grid_is_certified, solve_vandermonde from .conversion.edges import create_edge_tree from .validation import validate_tree_model @@ -711,20 +712,32 @@ def _get_n_matrix(interpolated_poly: np.ndarray) -> np.ndarray: """ depth = interpolated_poly.shape[0] Ns = np.zeros((depth + 1, depth)) + diagnostics = VandermondeDiagnostics() + certified = grid_is_certified(interpolated_poly) for i in range(1, depth + 1): - Ns[i, :i] = np.linalg.inv(np.vander(interpolated_poly[:i]).T).dot( - 1.0 / np.array([sp.special.binom(i - 1, k) for k in range(i)]) + Ns[i, :i] = solve_vandermonde( + interpolated_poly[:i], + 1.0 / np.array([sp.special.binom(i - 1, k) for k in range(i)]), + certified=certified, + diagnostics=diagnostics, ) + diagnostics.emit(stacklevel=5) return Ns def _get_n_cii_matrix(self, interpolated_poly: np.ndarray, order: int) -> np.ndarray: """Computes the N matrix for the CII index.""" depth = interpolated_poly.shape[0] Ns = np.zeros((depth + 1, depth)) + diagnostics = VandermondeDiagnostics() + certified = grid_is_certified(interpolated_poly) for i in range(1, depth + 1): - Ns[i, :i] = np.linalg.inv(np.vander(interpolated_poly[:i]).T).dot( + Ns[i, :i] = solve_vandermonde( + interpolated_poly[:i], i * np.array([self._get_subset_weight_cii(j, order) for j in range(i)]), + certified=certified, + diagnostics=diagnostics, ) + diagnostics.emit(stacklevel=5) return Ns def _get_subset_weight_cii(self, t: int, order: int) -> float | None: @@ -758,8 +771,14 @@ def _get_n_id_matrix(D: np.ndarray) -> np.ndarray: """Computes N_id matrix.""" depth = D.shape[0] Ns_id = np.zeros((depth + 1, depth)) + diagnostics = VandermondeDiagnostics() + certified = grid_is_certified(D) for i in range(1, depth + 1): - Ns_id[i, :i] = np.linalg.inv(np.vander(D[:i]).T).dot(np.ones(i)) + Ns_id[i, :i] = solve_vandermonde( + D[:i], np.ones(i), certified=certified, diagnostics=diagnostics + ) + # stacklevel traces emit <- builder <- _init_summary_polynomials <- __init__ <- user + diagnostics.emit(stacklevel=5) return Ns_id @staticmethod diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py new file mode 100644 index 000000000..f2615c502 --- /dev/null +++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py @@ -0,0 +1,94 @@ +"""Tests for the guarded Chebyshev-Vandermonde solves in ``shapiq.tree._numerics``.""" + +from __future__ import annotations + +import warnings + +import numpy as np +import pytest +from numpy.polynomial.chebyshev import chebpts2 + +from shapiq.tree._numerics import grid_is_certified, solve_vandermonde +from shapiq.tree.linear.explainer import LinearTreeSHAP, get_N_v2 + + +@pytest.mark.parametrize("grid_size", [4, 8, 12, 20, 26]) +def test_matches_explicit_inverse_on_well_conditioned_grids(grid_size): + """The fast path agrees with the previous ``inv(V) @ rhs`` formulation.""" + D = chebpts2(grid_size) + rng = np.random.default_rng(0) + for i in range(2, grid_size + 1): + rhs = rng.normal(size=i) + V = np.vander(D[:i]).T + expected = np.linalg.inv(V).dot(rhs) + with warnings.catch_warnings(): + warnings.simplefilter("error") # no warning may fire on the fast path + result = solve_vandermonde(D[:i], rhs, certified=grid_is_certified(D)) + np.testing.assert_allclose(result, expected, rtol=1e-8, atol=1e-10) + + +def test_precision_warning_fires_at_interior_prefixes_of_large_grids(): + """Conditioning peaks at interior prefixes (i ~ n/2): for a size-30 grid the + warning must fire even though the degenerate prefix is shorter than the safe + grid size — this is the regression test for the mis-calibrated size gate.""" + D = chebpts2(28) + i = 14 # peak conditioning region; full-rank with cond ~4e12 (comfortably inside + # the precision-warning band: rank-deficient prefixes only appear from grid ~32) + with pytest.warns(RuntimeWarning, match="condition number"): + solve_vandermonde(D[:i], np.ones(i)) + + +def test_no_warning_for_ordinary_grids(): + D = chebpts2(12) + with warnings.catch_warnings(): + warnings.simplefilter("error") + for i in range(2, 13): + solve_vandermonde(D[:i], np.ones(i)) + + +def test_singular_system_returns_least_squares_fallback_with_warning(): + """A rank-deficient system must not crash; it returns a least-squares solution + and warns that the values are not reliable.""" + points = np.concatenate([chebpts2(30), chebpts2(30)[:30]]) # duplicated nodes + size = len(points) + with pytest.warns(RuntimeWarning, match="NOT reliable"): + result = solve_vandermonde(points, np.ones(size), certified=grid_is_certified(points)) + assert result.shape == (size,) + assert np.all(np.isfinite(result)) + + +def test_deep_tree_constructs_and_explains_with_warning(): + """End-to-end: LinearTreeSHAP on a deep tree no longer raises LinAlgError.""" + sklearn = pytest.importorskip("sklearn.tree") + + rng = np.random.default_rng(0) + X = rng.normal(size=(20000, 10)) + y = rng.normal(size=20000) + tree = sklearn.DecisionTreeRegressor(max_depth=60, min_samples_leaf=1, random_state=0).fit(X, y) + if tree.get_depth() < 50: + pytest.skip("random data did not produce a deep enough tree") + with pytest.warns(RuntimeWarning): + explainer = LinearTreeSHAP(tree) + values = explainer.explain_function(X[0]) + assert np.all(np.isfinite(values.values)) + + +def test_one_summary_warning_per_matrix_construction(): + """A deep grid touches many degenerate prefixes; the construction loop must + emit exactly one coalesced warning, not one per prefix.""" + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + get_N_v2(chebpts2(30)) + runtime_warnings = [w for w in caught if issubclass(w.category, RuntimeWarning)] + assert len(runtime_warnings) == 1 + + +def test_custom_clustered_grid_is_not_silently_fast_pathed(): + """Certification measures the actual grid: a small but badly conditioned custom + grid (legal via LinearTreeSHAP's base_func) must not take the unchecked path.""" + grid = np.linspace(0.0, 1.0, 26) ** 8 # clustered near zero, size <= 26 + assert not grid_is_certified(grid) + assert grid_is_certified(chebpts2(26)) + with pytest.warns(RuntimeWarning): + for i in range(2, len(grid) + 1): + solve_vandermonde(grid[:i], np.ones(i), certified=grid_is_certified(grid)) From 2828ba43cb8fb620538cf15b74b5d48c4d8d405f Mon Sep 17 00:00:00 2001 From: YUANLUO WU Date: Wed, 10 Jun 2026 23:44:00 +0200 Subject: [PATCH 2/3] fix(tree): solve the TreeSHAP Vandermonde systems exactly at every depth Replaces the guarded float64 solves with exact ones, removing all warning and least-squares fallback paths: - O(n^2) Bjorck-Pereyra dual recursion in scaled-integer fixed-point arithmetic (standard-library int), with a convergence certificate: each system is solved at increasing precision until two independent computations agree bitwise, so the result is the float64 rounding of the exact rational solution at any depth (a depth-100 grid's full prefix workload takes ~0.3 s). Nodes that collapse onto the same scaled integer at the first precision rung climb the ladder instead of being misreported as coincident. - The previous inversion drifted at the ~1e-7 level from interpolation degree ~20, returned silently wrong values from ~32, and crashed (LinAlgError) at ~60+; the exact solves carry no conditioning error at any degree. - The one remaining limit is representational: the monomial-basis N entries grow exponentially with the interpolation degree and the downstream float64 pipeline cancels them, so beyond a measured magnitude bound (degree ~29 on the default grids, pinned by a boundary test) construction raises an explanatory RepresentationLimitError instead of returning silently wrong values; for order-1 Shapley values TreeExplainer re-routes affected trees to TreeSHAPIQ when its feature-bounded degree still fits. - Tests: bitwise agreement with an exact rational-elimination oracle (including the previously rank-deficient sizes 28-45 and a clustered custom grid), the exact 29/30 accept/refuse boundary, deep-chain completeness at depths 20/24/28, refusal at 35/60 on both the LinearTreeSHAP and TreeSHAPIQ paths, re-routing for SV and order-1 SII, precision-ladder climbing for nodes finer than the first rung's resolution, and cache isolation. --- CHANGELOG.md | 12 +- docs/source/_templates/autosummary/module.rst | 11 + src/shapiq/__init__.py | 7 +- src/shapiq/tree/__init__.py | 3 + src/shapiq/tree/_numerics.py | 427 +++++++++-------- src/shapiq/tree/explainer.py | 103 ++--- src/shapiq/tree/linear/explainer.py | 19 +- src/shapiq/tree/treeshapiq.py | 53 +-- src/shapiq/tree/validation.py | 9 +- src/shapiq/utils/__init__.py | 3 +- src/shapiq/utils/errors.py | 14 + .../test_tree_explainer.py | 2 +- .../test_tree_explainer_validate.py | 6 + .../test_tree_numerics.py | 432 +++++++++++++++--- 14 files changed, 716 insertions(+), 385 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index bd262a639..07cdc910b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,15 +4,15 @@ ### Bugfix -- Hardens the Vandermonde interpolation solves in the polynomial TreeSHAP machinery (`TreeSHAPIQ`, `LinearTreeSHAP`): - - **surfaces a pre-existing silent-wrong-answer bug**: for interpolation degrees ~27-31 (for `TreeSHAPIQ` the degree is `min(depth, n_features)`), the previous explicit matrix inversion silently lost precision - a coalesced `RuntimeWarning` now reports the worst condition number; from degree ~32 the prefix systems become rank-deficient and the returned values can violate the completeness axiom by errors up to orders of magnitude beyond the prediction itself - such cases now emit a clear `RuntimeWarning` stating the values are NOT reliable, - - very deep trees no longer crash with an unexplained `LinAlgError` (previously at degree ~60); a least-squares fallback is returned together with the reliability warning, - - no explicit matrix inverse is formed anymore (`np.linalg.solve` / SVD-based least squares), - - interpolation grids are certified by measurement (cached per grid, hoisted out of the construction loops), so custom `LinearTreeSHAP(base_func=...)` grids are handled correctly and ordinary trees pay no per-solve diagnostic work. [#545](https://github.com/mmschlk/shapiq/issues/545) +- Solves the Vandermonde interpolation systems of the polynomial TreeSHAP machinery (`TreeSHAPIQ`, `LinearTreeSHAP`) exactly at every depth, replacing the explicit float64 inverse: + - an O(n^2) Björck-Pereyra recursion in scaled-integer arithmetic with a bitwise convergence certificate returns the float64 rounding of the exact rational solution for any grid of distinct nodes, + - this **fixes a pre-existing silent-wrong-answer bug**: the previous inversion drifted at the ~1e-7 level from interpolation degree ~20, returned silently wrong values from ~32, and crashed with an unexplained `LinAlgError` at ~60+; the exact solves carry no conditioning error at any degree, + - beyond the float64 representation limit of the monomial-basis pipeline (interpolation degree ~29 for `LinearTreeSHAP` and ~25 for `TreeSHAPIQ` on the default grids, where the exact coefficients exceed ~3e10 and downstream double-precision evaluation provably cancels more than ~1% of the result) construction raises an explanatory `RepresentationLimitError` (a `ValueError`) instead of returning silently wrong values. [#545](https://github.com/mmschlk/shapiq/issues/545) ### Changed -- `LinearTreeSHAP` (trees deeper than ~26) and `TreeSHAPIQ` (interpolation degree above ~26) now emit a `RuntimeWarning` on paths that previously were silent (precision loss / unreliable deep-tree values); test suites running with `filterwarnings = error` will surface this intentionally. Likely warrants a minor (not patch) version bump. +- Models whose interpolation degree exceeds the float64 representation limit (~29 for `LinearTreeSHAP`, ~25 for `TreeSHAPIQ`, on the default grids) now raise `RepresentationLimitError` (a `ValueError`) at construction; previously they returned silently wrong values (degrees ~32-59) or crashed with `LinAlgError` (~60+). The degree depends on the configuration: the tree depth for `LinearTreeSHAP`, `min(depth, features in the tree)` for `TreeSHAPIQ` with SV/SII/k-SII, and the number of features in the tree for STII/FSII/BII. For order-1 SV/SII configurations `TreeExplainer` re-routes the affected trees to `TreeSHAPIQ` when its feature-bounded degree still fits; other configurations propagate the error. +- Shapley values for interpolation degrees above ~20 may shift at the ~1e-7 level relative to previous releases: the coefficients are now the correctly rounded exact solutions rather than the output of a float64 inverse. Below the representation limit the remaining end-to-end completeness error is bounded by the downstream float64 evaluation (about `max|N| * 1e-13`, i.e. up to ~1e-4 at degree 24 and ~1e-2 at degree 28 on adversarial chain trees). ## v1.5.1 (2026-05-30) diff --git a/docs/source/_templates/autosummary/module.rst b/docs/source/_templates/autosummary/module.rst index 40dbdeff9..ec19e6c25 100644 --- a/docs/source/_templates/autosummary/module.rst +++ b/docs/source/_templates/autosummary/module.rst @@ -14,6 +14,17 @@ {%- endfor %} {% endif %} +{% if exceptions %} +.. rubric:: Exceptions + +.. autosummary:: + :toctree: + :nosignatures: +{% for item in exceptions %} + {{ item }} +{%- endfor %} +{% endif %} + {% if functions %} .. rubric:: Functions diff --git a/src/shapiq/__init__.py b/src/shapiq/__init__.py index 140804b74..d3ecf02ed 100644 --- a/src/shapiq/__init__.py +++ b/src/shapiq/__init__.py @@ -79,8 +79,9 @@ ) from .tree import TreeExplainer -# public utils functions -from .utils import ( # sets.py # tree.py +# public utils functions and errors +from .utils import ( + RepresentationLimitError, get_explicit_subsets, powerset, safe_isinstance, @@ -114,6 +115,8 @@ "ProxySHAP", "ProxySPEX", "RegressionMSR", + # errors + "RepresentationLimitError", # explainers "Explainer", "TabularExplainer", diff --git a/src/shapiq/tree/__init__.py b/src/shapiq/tree/__init__.py index 547da8690..968fa19be 100644 --- a/src/shapiq/tree/__init__.py +++ b/src/shapiq/tree/__init__.py @@ -4,6 +4,9 @@ along with supporting data structures and algorithm variants. """ +# re-exported for convenience; documented (and listed in __all__) only in shapiq.utils +from shapiq.utils.errors import RepresentationLimitError # noqa: F401 + from .base import TreeModel from .interventional import InterventionalGame, InterventionalTreeExplainer from .linear import LinearTreeSHAP diff --git a/src/shapiq/tree/_numerics.py b/src/shapiq/tree/_numerics.py index edb3c8f86..924a9b3ce 100644 --- a/src/shapiq/tree/_numerics.py +++ b/src/shapiq/tree/_numerics.py @@ -1,207 +1,262 @@ -"""Numerically guarded solves for the Vandermonde systems of TreeSHAP-IQ. - -The polynomial TreeSHAP machinery repeatedly needs ``inv(vander(points).T) @ rhs``, -where ``points`` is a prefix ``D[:i]`` of a fixed interpolation grid ``D`` whose size -equals the interpolation degree (``min(tree depth, n_features_in_tree)`` for -TreeSHAPIQ, the full tree depth for LinearTreeSHAP). These Vandermonde systems are -severely ill-conditioned, and the conditioning is **non-monotonic in the prefix -length**: the first ``i`` points of an ``n``-point grid typically cluster near one -endpoint, so the condition number peaks around ``i ~ n/2``. For the default -Chebyshev grids, measured peak prefix condition numbers are ``4.1e11`` at grid size -26, ``1.3e12`` at 27 and ``4.2e13`` at 30 — precision degrades for grids larger -than ~26 (at *interior* prefixes, not only at the full size), and the systems become -singular to machine precision for very deep trees (sizes around 55-60), which -previously surfaced as an unexplained ``LinAlgError``. - -This module centralises the solve: - -- **Grid certification instead of assumptions.** For each distinct grid, the peak - condition number over all prefixes is measured once and cached - (:func:`grid_is_certified`); call sites hoist this decision out of their - construction loops, so certified grids take a plain ``np.linalg.solve`` with no - per-solve diagnostic work. The certification is computed on the actual grid, so - custom interpolation bases (``LinearTreeSHAP(base_func=...)``) are handled - correctly rather than assumed to behave like Chebyshev nodes. -- **Checked path with one SVD.** Uncertified grids are solved by SVD-based least - squares; the singular values of that same decomposition give the exact condition - number of the matrix actually solved, and rank-deficient systems degrade into a - least-squares solution instead of crashing. -- **One warning per matrix, not per prefix.** Call sites construct whole N-matrices - in a loop over prefixes; per-prefix warnings would flood the user (a depth-60 - construction touches dozens of degenerate prefixes). The loop passes a - :class:`VandermondeDiagnostics` collector and emits a single summary warning. - -Design note: rank-deficient systems deliberately **warn and return** the -least-squares fallback rather than raising. Before this module, the same depths -either crashed (``LinAlgError``, degree ~60) or — worse — **silently returned -wrong values** (rank-deficient prefixes appear from grid size ~32; completeness-axiom errors up to -orders of magnitude beyond the prediction itself). Warning loudly while degrading is strictly -more informative than either previous behaviour, and pipelines that want hard -failures can promote the warning with the ``warnings`` module. +"""Exact solves for the Vandermonde systems of the polynomial TreeSHAP machinery. + +``TreeSHAPIQ`` and ``LinearTreeSHAP`` repeatedly need ``inv(vander(D[:i]).T) @ rhs`` +for prefixes ``D[:i]`` of a fixed interpolation grid. These Vandermonde systems are +severely ill-conditioned in double precision — the conditioning is non-monotonic in +the prefix length and peaks at interior prefixes (``i ~ n/2``), so explicit +``float64`` inversion drifts at the ~1e-7 level from interpolation degree ~20, +returns wrong values from ~32 on the default grids, and raises ``LinAlgError`` +for very deep trees (~60+). + +The conditioning is purely a floating-point artifact: the interpolation nodes are +distinct, so the systems are exactly solvable over the rationals. This module +therefore solves them **exactly** and returns the correctly rounded ``float64`` +result, at every depth: + +- **O(n^2) Björck-Pereyra recursion** instead of an O(n^3) inverse. The dual + Vandermonde recursion (Golub & Van Loan, Algorithm 4.6.2) factors the inverse + into bidiagonal steps, which keeps the exact-arithmetic operand sizes small. +- **Scaled-integer fixed-point arithmetic** (standard-library ``int``): every node + and intermediate value is represented as ``round(v * 2**bits)``. With ``bits`` + large enough the rounded result equals the exact rational solution rounded to + ``float64``; integer arithmetic keeps a depth-100 grid's full prefix workload in + the sub-second range, where naive ``fractions.Fraction`` elimination needs minutes. +- **A convergence certificate instead of an error analysis**: each system is solved + at increasing precision (128, 256, ... bits) until two consecutive precision + levels produce bitwise-identical ``float64`` outputs. Agreement of two scaled + truncation grids that fine means the value has converged to the rounding of the + exact solution, so the returned coefficients carry no conditioning error at all. + +Every path either returns the exact result or refuses loudly: degenerate inputs +(coincident, non-finite, or extreme custom nodes) raise ``ValueError``, and +N-matrix rows whose magnitude exceeds what the downstream float64 pipeline can +evaluate raise :class:`~shapiq.utils.errors.RepresentationLimitError`. """ from __future__ import annotations -import warnings +from fractions import Fraction from functools import lru_cache +from typing import TYPE_CHECKING import numpy as np -# Above this condition number, double precision has lost ~12 of its ~16 significant -# digits; results may no longer be exact. (Empirical anchor for the default -# Chebyshev grids: the peak prefix condition number stays below this threshold up -# to grid size 26 and exceeds it from size 27 on; see module docstring. The -# certification below measures each actual grid, so the constant is a warning -# threshold, not a structural assumption.) -_COND_WARN_THRESHOLD = 1e12 - -# Certification is only attempted for grids up to this size; larger grids go -# straight to the checked path (their peak prefix conditioning exceeds the -# threshold for every node family used in practice, and bounding the cache work -# keeps certification O(small)). -_CERTIFY_MAX_SIZE = 32 - - -@lru_cache(maxsize=128) -def _max_prefix_condition(grid_key: tuple[float, ...]) -> float: - """Measured peak condition number over all prefixes of the given grid.""" - grid = np.asarray(grid_key) - worst = 1.0 - for i in range(2, len(grid) + 1): - try: - worst = max(worst, float(np.linalg.cond(np.vander(grid[:i]).T))) - except np.linalg.LinAlgError: # pragma: no cover - cond rarely raises - return float("inf") - return worst +from shapiq.utils.errors import RepresentationLimitError +if TYPE_CHECKING: + from collections.abc import Callable -class VandermondeDiagnostics: - """Collects conditioning diagnostics across the solves of one N-matrix. +# Precision ladder for the convergence certificate. 128 bits converges for every +# Chebyshev-grid prefix up to at least depth 100; the 256-bit rung confirms it +# whenever the cheap float64 cross-check cannot (ill-conditioned prefixes), and +# the higher rungs exist for adversarial custom grids (LinearTreeSHAP(base_func=...)). +_PRECISION_BITS = (128, 256, 512, 1024, 2048, 4096) - Use one instance per matrix-construction loop and call :meth:`emit` once after - the loop, so the user sees a single summary warning instead of one warning per - degenerate prefix. - """ +# Largest N-matrix entry the downstream float64 pipeline can absorb. The +# coefficients themselves are exact, but the explainers consume them in float64 +# inner products whose summands scale with the largest entry while the result is +# O(prediction). Measured end to end on chain trees of depth 20-40, the relative +# cancellation error tracks ``max|N| * 1e-13`` within one order of magnitude, so +# at 3e10 the expected loss is in the 0.3%-3% range. Beyond that the explainers +# refuse instead of returning values whose error may exceed ~1% of the +# prediction. For the default Chebyshev grids this corresponds to an +# interpolation degree of ~29 for LinearTreeSHAP's right-hand sides and ~25 +# for TreeSHAPIQ's (its identity N matrix, built for every index, saturates +# first). The bound is empirical, not a worst-case proof. +_REPRESENTATION_LIMIT = 3.0e10 + + +def _check_row_magnitude(row: np.ndarray) -> None: + """Refuse N-matrix rows whose magnitude exceeds the float64 pipeline's limit. - def __init__(self) -> None: - self.worst_condition: float = 0.0 - self.rank_deficient_sizes: list[int] = [] - self.worst_residual: float = 0.0 - - def emit(self, stacklevel: int = 4) -> None: - """Emit at most one summary warning for the collected diagnostics. - - Args: - stacklevel: Forwarded to :func:`warnings.warn`; pass the depth of the - user's call site relative to ``emit`` (the construction loops sit - at different depths in LinearTreeSHAP vs TreeSHAPIQ). - """ - if self.rank_deficient_sizes: - warnings.warn( - f"Interpolation systems of size {self.rank_deficient_sizes} are " - "singular to machine precision (very deep tree). Least-squares " - "fallbacks were used (worst relative residual " - f"{self.worst_residual:.1e}) - the resulting explanation values are " - "NOT reliable. Limit the tree depth to obtain exact values.", - RuntimeWarning, - stacklevel=stacklevel, - ) - elif self.worst_condition > _COND_WARN_THRESHOLD: - warnings.warn( - "The interpolation systems for this tree reach a condition number " - f"of {self.worst_condition:.1e}; TreeSHAP-IQ values may lose " - "precision. Consider limiting the tree depth (the interpolation " - "degree is min(tree depth, number of features in the tree)).", - RuntimeWarning, - stacklevel=stacklevel, - ) - - -def grid_is_certified(grid: np.ndarray) -> bool: - """Whether every prefix system of ``grid`` is measured well-conditioned. - - Note that "certified" bounds, but does not eliminate, precision loss: the - threshold admits condition numbers up to 1e12 (a size-26 Chebyshev grid peaks - at ~4e11). Certification costs one SVD per prefix on first touch of a grid - (cached afterwards, and shared across the N-matrix builders), roughly doubling - the linear-algebra work of the first construction for that grid size. - - Call once per interpolation grid (the measurement is cached per grid) and pass - the result to :func:`solve_vandermonde` as ``certified=`` so the per-prefix - solves carry no diagnostic work at all for ordinary trees. + The Vandermonde solves themselves are exact at any depth, but the + monomial-basis N-matrix entries grow exponentially with the interpolation + degree, and the downstream float64 evaluation cancels them against each + other. Beyond ``_REPRESENTATION_LIMIT`` the empirically calibrated loss may + exceed ~1% of the result, so this raises instead of letting the explainer + return silently wrong values. + + Raises: + RepresentationLimitError: If the row magnitude exceeds the limit. """ - grid = np.asarray(grid) - if len(grid) > _CERTIFY_MAX_SIZE: - return False - return _max_prefix_condition(tuple(grid.tolist())) <= _COND_WARN_THRESHOLD - - -def solve_vandermonde( - points: np.ndarray, - rhs: np.ndarray, - *, - certified: bool = False, - diagnostics: VandermondeDiagnostics | None = None, + peak = float(np.max(np.abs(row))) + if peak > _REPRESENTATION_LIMIT: + msg = ( + "The interpolation degree of this tree is too large for the float64 " + f"TreeSHAP pipeline: the exact interpolation coefficients reach {peak:.1e}, " + "and evaluating them in double precision would lose more than ~1% of the " + "result to cancellation. Reduce the interpolation degree — depending on " + "the explainer and index this is the tree depth, min(tree depth, features " + "in the tree), or the number of features in the tree; the default grids " + "support up to ~29 with LinearTreeSHAP and ~25 with TreeSHAPIQ. This is a " + "structural limit of the monomial-basis representation, not of the solver." + ) + raise RepresentationLimitError(msg) + + +def build_n_matrix( + grid: np.ndarray, + rhs_for_row: Callable[[int], np.ndarray], ) -> np.ndarray: - """Solve ``vander(points).T @ x = rhs`` without forming an explicit inverse. + """Build an interpolation N matrix row by row from exact Vandermonde solves. + + Row ``i`` is the solution of ``vander(grid[:i]).T @ x = rhs_for_row(i)``. + Every row is checked against the float64 representation limit as soon as it + is computed, so an over-deep tree is refused after the first offending row + instead of after the full (expensive) construction. + + Args: + grid: The full interpolation grid (its length is the maximum degree). + rhs_for_row: Callback returning the length-``i`` right-hand side of row ``i``. + + Returns: + The ``(len(grid) + 1, len(grid))`` N matrix (row 0 stays zero). + + Raises: + RepresentationLimitError: If any row exceeds the representation limit. + ValueError: For degenerate grids (coincident, non-finite, or extreme + nodes), propagated from :func:`solve_vandermonde`. + """ + depth = len(grid) + n_matrix = np.zeros((depth + 1, depth)) + for i in range(1, depth + 1): + n_matrix[i, :i] = solve_vandermonde(grid[:i], rhs_for_row(i)) + _check_row_magnitude(n_matrix[i, :i]) + return n_matrix + + +def _to_fixed_point(value: float, bits: int) -> int: + """Return ``value * 2**bits`` as an integer (truncated toward minus infinity). + + The truncation direction is irrelevant for the convergence certificate: the + per-step error stays below ``2**-bits`` in either case, and the certificate + only accepts results that are stable under a 128-bit precision increase. + """ + fraction = Fraction(value) + return (fraction.numerator << bits) // fraction.denominator + + +def _bjorck_pereyra_fixed_point(points: np.ndarray, rhs: np.ndarray, bits: int) -> np.ndarray: + """Solve ``vander(points).T @ x = rhs`` in scaled-integer arithmetic. + + Implements the dual Björck-Pereyra recursion on values scaled by ``2**bits``. + ``np.vander`` uses decreasing powers, so the right-hand side is reversed to + obtain the increasing-power moment problem the recursion expects. + """ + n = len(points) + nodes = [_to_fixed_point(p, bits) for p in points.tolist()] + x = [_to_fixed_point(v, bits) for v in rhs.tolist()[::-1]] + for k in range(n - 1): + node_k = nodes[k] + for i in range(n - 1, k, -1): + x[i] = x[i] - (node_k * x[i - 1] >> bits) + for k in range(n - 2, -1, -1): + for i in range(k + 1, n): + x[i] = (x[i] << bits) // (nodes[i] - nodes[i - k - 1]) + for i in range(k + 1, n): + x[i - 1] = x[i - 1] - x[i] + scale = 1 << bits + return np.array([float(Fraction(value, scale)) for value in x]) + + +def solve_vandermonde(points: np.ndarray, rhs: np.ndarray) -> np.ndarray: + """Solve ``vander(points).T @ x = rhs`` exactly, rounded to ``float64``. + + The system is solved in scaled-integer arithmetic and certified by agreement + with an independent computation: first against a plain float64 LAPACK solve + (which corroborates the well-conditioned common case at negligible cost), and + otherwise against the next rung of the precision ladder, until two + computations agree bitwise. An agreeing result equals the exact rational + solution rounded to ``float64`` up to an overwhelmingly improbable + coincidence of correlated rounding (both computations would have to land in + the same wrong rounding bin), regardless of the conditioning of the system. + Works for any grid of distinct interpolation nodes. + + Solves are memoized on the byte representation of ``(points, rhs)``: the + explainers issue identical solves for every tree of equal depth in an + ensemble, which therefore costs one tree's worth of work. Args: - points: Interpolation nodes (a prefix of the interpolation grid). + points: Interpolation nodes (a prefix of the interpolation grid). The nodes + must be pairwise distinct. rhs: Right-hand side vector of the same length. - certified: Result of :func:`grid_is_certified` for the *full* interpolation - grid that ``points`` is a prefix of (conditioning peaks at interior - prefixes, so safety is a property of the whole grid). When ``True``, a - plain solve is used with no diagnostic work. Defaults to ``False`` - (checked SVD-based path). - diagnostics: Optional collector. If given, conditioning findings are - recorded on it (call :meth:`VandermondeDiagnostics.emit` after the - construction loop) instead of warning per call. Returns: - The solution vector ``x``. + The solution vector ``x`` (the ``float64`` rounding of the exact solution). - Warns: - RuntimeWarning: Only when ``diagnostics`` is ``None``: if the solved - system's condition number exceeds ``1e12`` (possible precision loss), - or if it is rank-deficient and a least-squares fallback is returned - whose values are not reliable. + Raises: + ValueError: If the inputs are degenerate — coincident or non-finite + nodes, a non-finite right-hand side, or custom nodes so close + together (or so large) that the exact solution leaves the float64 + range or does not stabilize within the precision ladder. The + default Chebyshev grids never trigger any of these. """ - V = np.vander(points).T - size = len(points) - if certified: - # Certified grid: every prefix system measured well-conditioned. - return np.linalg.solve(V, rhs) - - # Checked path: one SVD yields the solution, the exact condition number of this - # very matrix, and rank-deficiency detection. - solution, _, rank, singular_values = np.linalg.lstsq(V, rhs, rcond=None) - if rank < size: - residual = float(np.linalg.norm(V @ solution - rhs)) - scale = float(np.linalg.norm(rhs)) or 1.0 - if diagnostics is not None: - diagnostics.rank_deficient_sizes.append(size) - diagnostics.worst_residual = max(diagnostics.worst_residual, residual / scale) - else: - warnings.warn( - f"The interpolation system of size {size} is singular to machine " - "precision (very deep tree). A least-squares fallback is returned " - f"with relative residual {residual / scale:.1e} - the resulting " - "explanation values are NOT reliable. Limit the tree depth to " - "obtain exact values.", - RuntimeWarning, - stacklevel=2, - ) - return solution - cond = float(singular_values[0] / singular_values[-1]) - if diagnostics is not None: - diagnostics.worst_condition = max(diagnostics.worst_condition, cond) - elif cond > _COND_WARN_THRESHOLD: - warnings.warn( - f"The interpolation system of size {size} has condition number " - f"{cond:.1e}; TreeSHAP-IQ values for trees this deep may lose " - "precision. Consider limiting the tree depth (the interpolation " - "degree is min(tree depth, number of features in the tree)).", - RuntimeWarning, - stacklevel=2, - ) - return solution + points = np.asarray(points, dtype=float) + rhs = np.asarray(rhs, dtype=float) + if points.ndim != 1 or points.shape != rhs.shape: + msg = "The interpolation nodes and right-hand side must be 1-D arrays of equal length." + raise ValueError(msg) + if not (np.isfinite(points).all() and np.isfinite(rhs).all()): + msg = "The interpolation nodes and right-hand side must be finite." + raise ValueError(msg) + if len(points) == 1: + return rhs.copy() + solution = _solve_vandermonde_cached(points.tobytes(), rhs.tobytes()) + return solution.copy() # the cached array must stay immutable + + +_DEGENERATE_NODES_MSG = ( + "The exact solution leaves the float64 range or does not stabilize within " + f"{_PRECISION_BITS[-1]} fractional bits; the interpolation nodes are too close " + "together or too large in magnitude. Use better-separated nodes of moderate " + "magnitude." +) + + +def clear_solver_cache() -> None: + """Drop all memoized Vandermonde solutions (for cold-timing tests or memory-sensitive hosts).""" + _solve_vandermonde_cached.cache_clear() + + +@lru_cache(maxsize=4096) +def _solve_vandermonde_cached(points_bytes: bytes, rhs_bytes: bytes) -> np.ndarray: + n = len(points_bytes) // 8 # float64 buffers; lengths are equal by construction + points = np.frombuffer(points_bytes, dtype=float, count=n) + rhs = np.frombuffer(rhs_bytes, dtype=float, count=n) + + previous = None + for bits in _PRECISION_BITS: + try: + current = _bjorck_pereyra_fixed_point(points, rhs, bits) + except ZeroDivisionError: + # Two nodes collapsed onto the same scaled integer at this precision. + # Genuinely distinct nodes separate at a finer rung (and stay separated + # at every rung above it), so climb instead of misreporting them as + # coincident; truly coincident nodes collapse at every rung and fall + # through to the error below. + previous = None + continue + except OverflowError: + # float(Fraction(...)) overflows when the exact solution leaves the + # float64 range (nearly-coincident or very large custom nodes). + raise ValueError(_DEGENERATE_NODES_MSG) from None + if previous is None: + # Cheap independent certificate: a float64 LAPACK solve corroborates + # the well-conditioned common case at a fraction of the cost of a + # second fixed-point solve. It must solve the same decreasing-power + # system the fixed-point recursion solves by reversing the rhs. + try: + float_solution = np.linalg.solve(np.vander(points).T, rhs) + if np.array_equal(current, float_solution): + return current + except np.linalg.LinAlgError: # float64 considers the system singular + pass + elif np.array_equal(previous, current): + return current + previous = current + if previous is None: + msg = "The interpolation nodes must be pairwise distinct." + raise ValueError(msg) + # Reachable only for custom grids with nearly-coincident nodes, where no + # float64 rounding of the (astronomically large) exact solution stabilizes. + raise ValueError(_DEGENERATE_NODES_MSG) diff --git a/src/shapiq/tree/explainer.py b/src/shapiq/tree/explainer.py index 6435233d3..12335cafc 100644 --- a/src/shapiq/tree/explainer.py +++ b/src/shapiq/tree/explainer.py @@ -11,6 +11,7 @@ from shapiq.explainer.base import Explainer from shapiq.tree.interventional.explainer import InterventionalTreeExplainer +from shapiq.utils.errors import RepresentationLimitError from .linear import LinearTreeSHAP from .treeshapiq import TreeSHAPIQ, TreeSHAPIQIndices @@ -97,6 +98,10 @@ def __init__( **kwargs: Additional keyword arguments are ignored. + Raises: + RepresentationLimitError: If a tree's interpolation degree exceeds the + float64 representation limit and no re-route applies (see + :class:`~shapiq.utils.errors.RepresentationLimitError`). """ super().__init__(model, index=index, max_order=max_order) @@ -116,28 +121,42 @@ def __init__( self.mode = mode self._reference_dataset: np.ndarray | None = reference_dataset - # In ``"pathdependent"`` mode, build exactly one per-tree explainer list — either - # ``LinearTreeSHAP`` (cheap, order-1 only) or ``TreeSHAPIQ`` (any order). The dispatch - # decision is fixed at construction time so callers can mutate the chosen list (e.g. - # ``_tree.thresholds`` rounding in tests) before calling :meth:`explain`. In + # In ``"pathdependent"`` mode, build per-tree explainers at construction time: + # ``LinearTreeSHAP`` (cheap, order-1 only) where possible, ``TreeSHAPIQ`` (any + # order) otherwise. The two lists can BOTH be populated when individual trees + # are re-routed (a tree whose depth-based interpolation degree exceeds the + # float64 representation limit while its feature-bounded TreeSHAPIQ degree + # does not); :meth:`explain` then aggregates across both lists. In # ``"interventional"`` mode no per-tree list is created — the - # :class:`~shapiq.tree.interventional.explainer.InterventionalTreeExplainer` handles the - # full ensemble in one shot, so a per-tree list would be meaningless. + # :class:`~shapiq.tree.interventional.explainer.InterventionalTreeExplainer` + # handles the full ensemble in one shot. self._treeshapiq_explainers: list[TreeSHAPIQ] = [] self._lineartreeshap_explainers: list[LinearTreeSHAP] = [] self._interventional_explainer: InterventionalTreeExplainer | None = None + # ``index`` (the local parameter) is already narrowed to ``TreeSHAPIQIndices``; + # ``self.index`` is the broader ``ExplainerIndices`` and would not type-check. + def _build_treeshapiq(tree: TreeModel, tree_index: TreeSHAPIQIndices) -> TreeSHAPIQ: + return TreeSHAPIQ(model=tree, max_order=self._max_order, index=tree_index) + if self.mode == "pathdependent": if self._can_use_lineartreeshap(): - self._lineartreeshap_explainers = [ - LinearTreeSHAP(model=tree) for tree in self._trees - ] + for tree in self._trees: + try: + self._lineartreeshap_explainers.append(LinearTreeSHAP(model=tree)) + except RepresentationLimitError: + # LinearTreeSHAP's interpolation degree is the full tree + # depth, while TreeSHAPIQ's is min(depth, features in the + # tree), so a deep tree over few features can still be + # explained there — at the same order-1 Shapley values. + # Only the affected tree leaves the fast path; if its + # TreeSHAPIQ degree also exceeds the limit, the error + # propagates. ``index="SV"`` matches LinearTreeSHAP's + # output label (order-1 indices all reduce to SV). + self._treeshapiq_explainers.append(_build_treeshapiq(tree, "SV")) else: - # ``index`` (the local parameter) is already narrowed to ``TreeSHAPIQIndices``; - # ``self.index`` is the broader ``ExplainerIndices`` and would not type-check. self._treeshapiq_explainers = [ - TreeSHAPIQ(model=tree, max_order=self._max_order, index=index) - for tree in self._trees + _build_treeshapiq(tree, index) for tree in self._trees ] elif self.mode == "interventional": if self._reference_dataset is None: @@ -173,45 +192,6 @@ def _can_use_lineartreeshap(self) -> bool: and all(tree.n_features_in_tree >= 2 for tree in self._trees) ) - def _explain_function_lineartreeshap( - self, - x: np.ndarray, - **kwargs: Any, # noqa: ARG002 - ) -> InteractionValues: - """Compute first-order Shapley values for ``x`` by aggregating the per-tree LinearTreeSHAP results. - - Mirrors the per-tree aggregation done by ``_explain_function_treeshapiq``: each - ``LinearTreeSHAP`` in ``self._lineartreeshap_explainers`` runs against ``x``, the - resulting :class:`~shapiq.interaction_values.InteractionValues` are summed (which also - sums ``baseline_value`` and the ``()`` entry), and ``min_order`` is finally enforced via - :meth:`InteractionValues.get_n_order` when the user asked for a stricter minimum. - - Args: - x: The instance to explain as a 1-dimensional array. - **kwargs: Additional keyword arguments are ignored. - - Returns: - The aggregated Shapley values for the instance. - """ - if len(x.shape) != 1: - msg = "explain expects a single instance, not a batch." - raise TypeError(msg) - - interaction_values: list[InteractionValues] = [ - lts.explain_function(x) for lts in self._lineartreeshap_explainers - ] - final_explanation = interaction_values[0] - for iv in interaction_values[1:]: - final_explanation += iv - - if self._min_order > final_explanation.min_order: - final_explanation = final_explanation.get_n_order( - min_order=self._min_order, - max_order=self._max_order, - ) - - return final_explanation - def _explain_function_interventionaltreeshapiq( self, x: np.ndarray, @@ -231,7 +211,7 @@ def _explain_function_interventionaltreeshapiq( raise RuntimeError(msg) return self._interventional_explainer.explain_function(x) - def _explain_function_treeshapiq( + def _explain_function_pathdependent( self, x: np.ndarray, **kwargs: Any, # noqa: ARG002 @@ -249,17 +229,19 @@ def _explain_function_treeshapiq( if len(x.shape) != 1: msg = "explain expects a single instance, not a batch." raise TypeError(msg) - # run treeshapiq for all trees - interaction_values: list[InteractionValues] = [] + # run the per-tree explainers; both lists can be populated when trees + # were re-routed at construction time (the other list is then empty) + interaction_values: list[InteractionValues] = [ + lts.explain_function(x) for lts in self._lineartreeshap_explainers + ] for explainer in self._treeshapiq_explainers: tree_explanation = explainer.explain(x) interaction_values.append(tree_explanation) # combine the explanations for all trees final_explanation = interaction_values[0] - if len(interaction_values) > 1: - for i in range(1, len(interaction_values)): - final_explanation += interaction_values[i] + for tree_explanation in interaction_values[1:]: + final_explanation += tree_explanation if self._min_order == 0 and final_explanation.min_order == 1: final_explanation.min_order = 0 @@ -295,8 +277,5 @@ def explain_function( # type: ignore[override] The computed interaction index for the instance. """ if self.mode == "pathdependent": - # Dispatch on whichever per-tree list __init__ chose to populate. - if self._lineartreeshap_explainers: - return self._explain_function_lineartreeshap(x, **kwargs) - return self._explain_function_treeshapiq(x, **kwargs) + return self._explain_function_pathdependent(x, **kwargs) return self._explain_function_interventionaltreeshapiq(x, **kwargs) diff --git a/src/shapiq/tree/linear/explainer.py b/src/shapiq/tree/linear/explainer.py index 9c469ba14..b639c6b1a 100644 --- a/src/shapiq/tree/linear/explainer.py +++ b/src/shapiq/tree/linear/explainer.py @@ -8,7 +8,7 @@ import scipy.special as sp from shapiq.interaction_values import InteractionValues -from shapiq.tree._numerics import VandermondeDiagnostics, grid_is_certified, solve_vandermonde +from shapiq.tree._numerics import build_n_matrix from shapiq.tree.conversion.edges import create_edge_tree from shapiq.tree.validation import validate_tree_model from shapiq.utils.sets import generate_interaction_lookup, powerset @@ -37,17 +37,7 @@ def get_N_prime(max_size: int = 10) -> np.ndarray: def get_N_v2(D: np.ndarray) -> np.ndarray: """Get N_v2 matrix for Linear Tree Shap.""" - depth = D.shape[0] - Ns = np.zeros((depth + 1, depth)) - diagnostics = VandermondeDiagnostics() - certified = grid_is_certified(D) - for i in range(1, depth + 1): - Ns[i, :i] = solve_vandermonde( - D[:i], 1.0 / get_norm_weight(i - 1), certified=certified, diagnostics=diagnostics - ) - # default stacklevel traces emit <- get_N_v2 <- __init__ <- user - diagnostics.emit() - return Ns + return build_n_matrix(D, lambda i: 1.0 / get_norm_weight(i - 1)) class LinearTreeSHAP: @@ -83,6 +73,11 @@ def __init__( :func:`~shapiq.tree.validation.validate_tree_model`. base_func: Callable ``int -> np.ndarray`` returning the interpolation base for the given depth. Defaults to :func:`numpy.polynomial.chebyshev.chebpts2`. + + Raises: + RepresentationLimitError: If the tree depth exceeds the float64 + representation limit of the interpolation pipeline (see + :class:`~shapiq.utils.errors.RepresentationLimitError`). """ self.clf = model self._tree = validate_tree_model(model, class_label=None)[0] diff --git a/src/shapiq/tree/treeshapiq.py b/src/shapiq/tree/treeshapiq.py index defe6b055..555fcd5a0 100644 --- a/src/shapiq/tree/treeshapiq.py +++ b/src/shapiq/tree/treeshapiq.py @@ -11,9 +11,10 @@ from shapiq.game_theory.indices import get_computation_index from shapiq.interaction_values import InteractionValues +from shapiq.utils.errors import RepresentationLimitError from shapiq.utils.sets import generate_interaction_lookup, powerset -from ._numerics import VandermondeDiagnostics, grid_is_certified, solve_vandermonde +from ._numerics import build_n_matrix from .conversion.edges import create_edge_tree from .validation import validate_tree_model @@ -154,6 +155,10 @@ def __init__( if self._n_features_in_tree > 0: try: self._init_summary_polynomials() + except RepresentationLimitError: + # a structural float64 limit, never the trivial single-feature + # case the broader ValueError handler below absorbs + raise except ValueError: if self._n_features_in_tree == 1: self._trivial_computation = True # for one feature the computation is trivial @@ -710,35 +715,17 @@ def _get_n_matrix(interpolated_poly: np.ndarray) -> np.ndarray: The N matrix. """ - depth = interpolated_poly.shape[0] - Ns = np.zeros((depth + 1, depth)) - diagnostics = VandermondeDiagnostics() - certified = grid_is_certified(interpolated_poly) - for i in range(1, depth + 1): - Ns[i, :i] = solve_vandermonde( - interpolated_poly[:i], - 1.0 / np.array([sp.special.binom(i - 1, k) for k in range(i)]), - certified=certified, - diagnostics=diagnostics, - ) - diagnostics.emit(stacklevel=5) - return Ns + return build_n_matrix( + interpolated_poly, + lambda i: 1.0 / np.array([sp.special.binom(i - 1, k) for k in range(i)]), + ) def _get_n_cii_matrix(self, interpolated_poly: np.ndarray, order: int) -> np.ndarray: """Computes the N matrix for the CII index.""" - depth = interpolated_poly.shape[0] - Ns = np.zeros((depth + 1, depth)) - diagnostics = VandermondeDiagnostics() - certified = grid_is_certified(interpolated_poly) - for i in range(1, depth + 1): - Ns[i, :i] = solve_vandermonde( - interpolated_poly[:i], - i * np.array([self._get_subset_weight_cii(j, order) for j in range(i)]), - certified=certified, - diagnostics=diagnostics, - ) - diagnostics.emit(stacklevel=5) - return Ns + return build_n_matrix( + interpolated_poly, + lambda i: i * np.array([self._get_subset_weight_cii(j, order) for j in range(i)]), + ) def _get_subset_weight_cii(self, t: int, order: int) -> float | None: """Computes the weight for a given subset size and interaction order. @@ -769,17 +756,7 @@ def _get_subset_weight_cii(self, t: int, order: int) -> float | None: @staticmethod def _get_n_id_matrix(D: np.ndarray) -> np.ndarray: """Computes N_id matrix.""" - depth = D.shape[0] - Ns_id = np.zeros((depth + 1, depth)) - diagnostics = VandermondeDiagnostics() - certified = grid_is_certified(D) - for i in range(1, depth + 1): - Ns_id[i, :i] = solve_vandermonde( - D[:i], np.ones(i), certified=certified, diagnostics=diagnostics - ) - # stacklevel traces emit <- builder <- _init_summary_polynomials <- __init__ <- user - diagnostics.emit(stacklevel=5) - return Ns_id + return build_n_matrix(D, np.ones) @staticmethod def _cache(interpolated_poly: FloatVector) -> FloatVector: diff --git a/src/shapiq/tree/validation.py b/src/shapiq/tree/validation.py index c5e96456d..c38d1b0f9 100644 --- a/src/shapiq/tree/validation.py +++ b/src/shapiq/tree/validation.py @@ -56,7 +56,9 @@ def validate_tree_model( Raises: TypeError: If the model type is not supported (raised from the underlying - ``NotImplementedError`` of :func:`~shapiq.tree.conversion.convert_tree_model`). + ``NotImplementedError`` of :func:`~shapiq.tree.conversion.convert_tree_model`), or + if the model yields no trees (an empty list, or a list whose elements are not + ``TreeModel`` instances). """ tree_model = [] # direct returns for base tree models and dict as model @@ -78,4 +80,9 @@ def validate_tree_model( msg = f"Model type {type(model)} is not supported." raise TypeError(msg) from e tree_model = result if isinstance(result, list) else [result] + if not tree_model: + # an empty list, or a list whose elements are not TreeModel instances, + # would otherwise propagate as an opaque IndexError at explain time + msg = "The model does not contain any trees: expected a non-empty (list of) TreeModel." + raise TypeError(msg) return tree_model diff --git a/src/shapiq/utils/__init__.py b/src/shapiq/utils/__init__.py index cc9747b64..6c1043305 100644 --- a/src/shapiq/utils/__init__.py +++ b/src/shapiq/utils/__init__.py @@ -1,7 +1,7 @@ """Low-level utility functions for working with subsets, coalitions, and modules.""" from .datasets import shuffle_data -from .errors import raise_deprecation_warning +from .errors import RepresentationLimitError, raise_deprecation_warning from .modules import check_import_module, safe_isinstance from .sets import ( count_interactions, @@ -32,5 +32,6 @@ # datasets "shuffle_data", # errors + "RepresentationLimitError", "raise_deprecation_warning", ] diff --git a/src/shapiq/utils/errors.py b/src/shapiq/utils/errors.py index 19c8650b8..06b0caacf 100644 --- a/src/shapiq/utils/errors.py +++ b/src/shapiq/utils/errors.py @@ -17,3 +17,17 @@ def raise_deprecation_warning( ) warn(message, DeprecationWarning, stacklevel=2) + + +class RepresentationLimitError(ValueError): + """The interpolation degree exceeds what the float64 TreeSHAP pipeline can evaluate. + + The exact interpolation coefficients are available at any depth, but the + downstream double-precision evaluation would cancel them beyond an acceptable + error. For order-1 SV/SII configurations, ``TreeExplainer`` catches this and + re-routes the affected tree to ``TreeSHAPIQ`` (whose interpolation degree is + bounded by the features in the tree rather than the tree depth); in every + other configuration — higher orders, the feature-bounded indices, or trees + that exceed the limit on both paths — the error propagates from the + constructor. + """ diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py index 17b318e66..58817026f 100644 --- a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py +++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py @@ -422,7 +422,7 @@ def test_xgboost_shap_error(xgb_clf_model, background_clf_data): # per-tree explainer list was populated so the mutation actually takes effect at explain time. per_tree_explainers = ( explainer_shapiq_rounded._lineartreeshap_explainers - or explainer_shapiq_rounded._treeshapiq_explainers + + explainer_shapiq_rounded._treeshapiq_explainers ) for tree_explainer in per_tree_explainers: tree_explainer._tree.thresholds = np.round(tree_explainer._tree.thresholds, 4) diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py index 89ed65829..b6c33fe77 100644 --- a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py +++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py @@ -33,6 +33,12 @@ def test_validate_model(dt_clf_model, dt_reg_model, rf_reg_model, rf_clf_model, with pytest.raises(TypeError): validate_tree_model("unsupported_model") + # treeless inputs raise instead of propagating an IndexError at explain time + with pytest.raises(TypeError, match="does not contain any trees"): + validate_tree_model([]) + with pytest.raises(TypeError, match="does not contain any trees"): + validate_tree_model([dt_clf_model, dt_reg_model]) # a list of raw sklearn models + @pytest.mark.external_libraries @pytest.mark.parametrize(("model_fixture", "model_class"), TREE_MODEL_FIXTURES) diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py index f2615c502..7a72ad02f 100644 --- a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py +++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py @@ -1,94 +1,374 @@ -"""Tests for the guarded Chebyshev-Vandermonde solves in ``shapiq.tree._numerics``.""" +"""Tests for the exact Vandermonde solves in ``shapiq.tree._numerics``.""" from __future__ import annotations -import warnings +from fractions import Fraction import numpy as np import pytest from numpy.polynomial.chebyshev import chebpts2 -from shapiq.tree._numerics import grid_is_certified, solve_vandermonde -from shapiq.tree.linear.explainer import LinearTreeSHAP, get_N_v2 +from shapiq.tree import RepresentationLimitError, TreeExplainer, TreeModel +from shapiq.tree._numerics import ( + _solve_vandermonde_cached, + clear_solver_cache, + solve_vandermonde, +) +from shapiq.tree.linear.explainer import LinearTreeSHAP, get_norm_weight +from shapiq.tree.treeshapiq import TreeSHAPIQ -@pytest.mark.parametrize("grid_size", [4, 8, 12, 20, 26]) +def _solve_exact_reference(points: np.ndarray, rhs: np.ndarray) -> np.ndarray: + """Oracle: exact rational Gaussian elimination for ``vander(points).T x = rhs``.""" + n = len(points) + p = [Fraction(t) for t in points.tolist()] + rows = [ + [p[j] ** (n - 1 - k) for j in range(n)] + [Fraction(v)] for k, v in enumerate(rhs.tolist()) + ] + for col in range(n): + pivot = next(r for r in range(col, n) if rows[r][col] != 0) + rows[col], rows[pivot] = rows[pivot], rows[col] + for r in range(n): + if r != col and rows[r][col] != 0: + factor = rows[r][col] / rows[col][col] + rows[r] = [a - factor * b for a, b in zip(rows[r], rows[col], strict=True)] + return np.array([float(rows[k][n] / rows[k][k]) for k in range(n)]) + + +@pytest.mark.parametrize("grid_size", [2, 4, 8, 12, 16, 20]) def test_matches_explicit_inverse_on_well_conditioned_grids(grid_size): - """The fast path agrees with the previous ``inv(V) @ rhs`` formulation.""" + """On well-conditioned systems the result agrees with the previous inverse. + + Larger grids are deliberately excluded here: above size ~20 the explicit + inverse starts drifting at the ~1e-7 level (at size 26 the peak prefix + condition number is already ~4e11), so the old formulation stops being a + valid reference; those sizes are covered by the exact-rational oracle test + below instead. + """ D = chebpts2(grid_size) rng = np.random.default_rng(0) for i in range(2, grid_size + 1): rhs = rng.normal(size=i) - V = np.vander(D[:i]).T - expected = np.linalg.inv(V).dot(rhs) - with warnings.catch_warnings(): - warnings.simplefilter("error") # no warning may fire on the fast path - result = solve_vandermonde(D[:i], rhs, certified=grid_is_certified(D)) - np.testing.assert_allclose(result, expected, rtol=1e-8, atol=1e-10) - - -def test_precision_warning_fires_at_interior_prefixes_of_large_grids(): - """Conditioning peaks at interior prefixes (i ~ n/2): for a size-30 grid the - warning must fire even though the degenerate prefix is shorter than the safe - grid size — this is the regression test for the mis-calibrated size gate.""" - D = chebpts2(28) - i = 14 # peak conditioning region; full-rank with cond ~4e12 (comfortably inside - # the precision-warning band: rank-deficient prefixes only appear from grid ~32) - with pytest.warns(RuntimeWarning, match="condition number"): - solve_vandermonde(D[:i], np.ones(i)) - - -def test_no_warning_for_ordinary_grids(): - D = chebpts2(12) - with warnings.catch_warnings(): - warnings.simplefilter("error") - for i in range(2, 13): - solve_vandermonde(D[:i], np.ones(i)) - - -def test_singular_system_returns_least_squares_fallback_with_warning(): - """A rank-deficient system must not crash; it returns a least-squares solution - and warns that the values are not reliable.""" - points = np.concatenate([chebpts2(30), chebpts2(30)[:30]]) # duplicated nodes - size = len(points) - with pytest.warns(RuntimeWarning, match="NOT reliable"): - result = solve_vandermonde(points, np.ones(size), certified=grid_is_certified(points)) - assert result.shape == (size,) - assert np.all(np.isfinite(result)) - - -def test_deep_tree_constructs_and_explains_with_warning(): - """End-to-end: LinearTreeSHAP on a deep tree no longer raises LinAlgError.""" - sklearn = pytest.importorskip("sklearn.tree") + expected = np.linalg.inv(np.vander(D[:i]).T).dot(rhs) + np.testing.assert_allclose(solve_vandermonde(D[:i], rhs), expected, rtol=1e-7, atol=1e-9) + + +@pytest.mark.parametrize("grid_size", [12, 28, 35, 45]) +def test_bitwise_equal_to_exact_rational_solution(grid_size): + """Every prefix solve equals the exact rational solution rounded to float64. + + Sizes 28-45 are exactly the regime where the previous ``inv``-based code lost + precision (peak prefix condition numbers above 1e12) or silently returned wrong + values (rank-deficient to machine precision from size ~32). + """ + D = chebpts2(grid_size) + rng = np.random.default_rng(1) + for i in (2, grid_size // 2, grid_size): # interior prefixes are the worst-conditioned + rhs = rng.normal(size=i) + exact = _solve_exact_reference(D[:i], rhs) + result = solve_vandermonde(D[:i], rhs) + np.testing.assert_array_equal(result, exact) + + +def test_replaced_inverse_construction_was_wrong_in_the_failure_band(): + """Regression guard: the construction this module replaced fails where the solver is exact. + + At grid size 32 with the library's own right-hand side the explicit inverse is + off by an absolute error of ~2e2 — entries of magnitude ~3e7 that downstream + inner products cancel to an O(prediction) result, so an absolute error of that + size destroys the result entirely. This test fails if the exact solver is ever + swapped back for the inverse. + """ + D = chebpts2(32) + rhs = 1.0 / get_norm_weight(31) + exact = _solve_exact_reference(D, rhs) + np.testing.assert_array_equal(solve_vandermonde(D, rhs), exact) + try: + old = np.linalg.inv(np.vander(D).T) @ rhs + assert np.max(np.abs(old - exact)) > 1.0 # measured ~2e2; anything >1 is fatal downstream + except np.linalg.LinAlgError: + pass # an outright failure is equally disqualifying for the old construction + + +def test_issue_545_scenario_raises_instead_of_returning_wrong_values(): + """The end-to-end scenario of issue #545 is refused instead of silently wrong. + + A decision tree fit on one-hot-style data (every sample distinguished by its + own indicator feature) is forced into a depth-39 chain over 39 features, the + structure that previously returned silently wrong Shapley values (interpolation + degree in the ~32-59 failure band, no warning). Both the LinearTreeSHAP degree + (tree depth) and the TreeSHAPIQ degree (features in the tree) exceed the + representation limit, so there is no re-route and construction must raise. + """ + sklearn_tree = pytest.importorskip("sklearn.tree") + x_data = np.eye(40) + y_target = np.arange(40, dtype=float) + model = sklearn_tree.DecisionTreeRegressor(random_state=0).fit(x_data, y_target) + assert model.get_depth() >= 30 # in the previously silent failure band + assert len({int(f) for f in model.tree_.feature if f >= 0}) >= 30 # no re-route possible + with pytest.raises(RepresentationLimitError): + TreeExplainer(model=model, index="SV", max_order=1) + + +def test_single_point_system(): + """The size-1 system is the identity.""" + np.testing.assert_array_equal(solve_vandermonde(np.array([0.3]), np.array([2.5])), [2.5]) + + +def test_non_chebyshev_clustered_grid_is_solved_exactly(): + """Custom interpolation grids (LinearTreeSHAP(base_func=...)) are handled exactly.""" + grid = np.linspace(0.0, 1.0, 18) ** 8 + 0.01 # clustered near zero + rng = np.random.default_rng(2) + rhs = rng.normal(size=18) + np.testing.assert_array_equal(solve_vandermonde(grid, rhs), _solve_exact_reference(grid, rhs)) + + +def test_coincident_nodes_raise(): + """Genuinely singular systems (repeated nodes) raise instead of returning garbage.""" + with pytest.raises(ValueError, match="pairwise distinct"): + solve_vandermonde(np.array([0.1, 0.5, 0.5]), np.ones(3)) + + +def test_mismatched_input_shapes_raise(): + """Mismatched or non-1-D inputs are rejected up front with a clear message.""" + with pytest.raises(ValueError, match="equal length"): + solve_vandermonde(np.array([0.1, 0.3, 0.5]), np.ones(2)) + with pytest.raises(ValueError, match="1-D"): + solve_vandermonde(np.ones((2, 2)), np.ones((2, 2))) + + +def test_non_finite_inputs_raise(): + """NaN or inf nodes/right-hand sides are rejected with a clear message.""" + with pytest.raises(ValueError, match="must be finite"): + solve_vandermonde(np.array([0.1, np.nan, 0.5]), np.ones(3)) + with pytest.raises(ValueError, match="must be finite"): + solve_vandermonde(np.array([0.1, 0.3, 0.5]), np.array([1.0, np.inf, 1.0])) + + +def test_nearly_coincident_nodes_raise_a_documented_error(): + """Nodes packed at consecutive float64 ulps push the exact solution beyond the + float64 range; the solver raises the documented ValueError instead of an opaque + OverflowError.""" + nodes = [0.5] + for _ in range(21): + nodes.append(float(np.nextafter(nodes[-1], 1.0))) + with pytest.raises(ValueError, match="too close together"): + solve_vandermonde(np.array(nodes), np.ones(len(nodes))) + + +def test_nodes_below_the_first_rung_resolution_climb_the_ladder(): + """Distinct nodes that collapse at 128 fractional bits are not misreported. + + Tiny-magnitude nodes can be genuinely distinct in float64 yet land on the + same scaled integer at the first precision rung; the solver must climb to a + finer rung and solve exactly rather than report them as coincident. + """ + a = 2.0**-100 + points = np.array([a, a + 2.0**-152]) # identical at 128 bits, distinct at 256 + rhs = np.array([1.0, 3.0]) + np.testing.assert_array_equal( + solve_vandermonde(points, rhs), _solve_exact_reference(points, rhs) + ) + + +def _deep_chain_tree(depth: int, n_features: int | None = None): + """A deterministic decision tree of exactly the requested depth. + Builds a left-leaning chain: every internal node splits on a fresh feature + (or cycles through ``n_features`` features when given) and has a leaf on the + right. The depth is exact by construction, so the deep-tree tests never + silently skip. + """ + n_nodes = 2 * depth + 1 + children_left = np.full(n_nodes, -1) + children_right = np.full(n_nodes, -1) + features = np.full(n_nodes, -2) + thresholds = np.full(n_nodes, np.nan) + values = np.zeros(n_nodes) + weights = np.ones(n_nodes) rng = np.random.default_rng(0) - X = rng.normal(size=(20000, 10)) - y = rng.normal(size=20000) - tree = sklearn.DecisionTreeRegressor(max_depth=60, min_samples_leaf=1, random_state=0).fit(X, y) - if tree.get_depth() < 50: - pytest.skip("random data did not produce a deep enough tree") - with pytest.warns(RuntimeWarning): + # Internal nodes 0..depth-1 form the left chain; node i's right child is the + # leaf depth+i, and the last internal node's left child is the leaf 2*depth. + for level in range(depth): + children_left[level] = level + 1 if level < depth - 1 else 2 * depth + children_right[level] = depth + level + features[level] = level if n_features is None else level % n_features + thresholds[level] = 0.5 + weights[level] = depth - level + 1 # number of leaves below the node + for leaf in range(depth, 2 * depth + 1): + values[leaf] = float(rng.normal()) + return TreeModel( + children_left=children_left, + children_right=children_right, + children_missing=children_left.copy(), + features=features, + thresholds=thresholds, + values=values, + node_sample_weight=weights, + ) + + +@pytest.mark.parametrize(("depth", "tolerance"), [(20, 1e-6), (24, 1e-4), (28, 1e-2)]) +def test_deep_tree_explanations_satisfy_completeness(depth, tolerance): + """Deep trees produce values satisfying the completeness axiom. + + The tolerances are the documented downstream float64 bound (``max|N| * 1e-13``): + with exact coefficients the only remaining error is the double-precision + evaluation of the (exponentially large) N entries, which scales with the + largest entry. + """ + tree = _deep_chain_tree(depth) + explainer = LinearTreeSHAP(tree) + x = np.zeros(depth) # walks the full left chain through every split + explanation = explainer.explain_function(x) + prediction = tree.predict_one(x) + gap = float(np.sum(explanation.values) - prediction) + assert abs(gap) < tolerance * max(1.0, abs(prediction)) + + +@pytest.mark.parametrize(("depth", "is_supported"), [(29, True), (30, False)]) +def test_representation_limit_boundary_is_pinned(depth, is_supported): + """The largest supported interpolation degree on the default grids is exactly 29. + + This pins the "~29" quoted in the error message, the CHANGELOG, and the + docstrings to the actual refuse/accept boundary, so a drift in the grids or + the limit cannot silently decouple the prose from the behavior. + """ + tree = _deep_chain_tree(depth) + if is_supported: explainer = LinearTreeSHAP(tree) - values = explainer.explain_function(X[0]) - assert np.all(np.isfinite(values.values)) - - -def test_one_summary_warning_per_matrix_construction(): - """A deep grid touches many degenerate prefixes; the construction loop must - emit exactly one coalesced warning, not one per prefix.""" - with warnings.catch_warnings(record=True) as caught: - warnings.simplefilter("always") - get_N_v2(chebpts2(30)) - runtime_warnings = [w for w in caught if issubclass(w.category, RuntimeWarning)] - assert len(runtime_warnings) == 1 - - -def test_custom_clustered_grid_is_not_silently_fast_pathed(): - """Certification measures the actual grid: a small but badly conditioned custom - grid (legal via LinearTreeSHAP's base_func) must not take the unchecked path.""" - grid = np.linspace(0.0, 1.0, 26) ** 8 # clustered near zero, size <= 26 - assert not grid_is_certified(grid) - assert grid_is_certified(chebpts2(26)) - with pytest.warns(RuntimeWarning): - for i in range(2, len(grid) + 1): - solve_vandermonde(grid[:i], np.ones(i), certified=grid_is_certified(grid)) + x = np.zeros(depth) + gap = float(np.sum(explainer.explain_function(x).values)) - tree.predict_one(x) + assert abs(gap) < 1e-2 * max(1.0, abs(tree.predict_one(x))) + else: + with pytest.raises(RepresentationLimitError): + LinearTreeSHAP(tree) + + +@pytest.mark.parametrize(("depth", "is_supported"), [(25, True), (26, False)]) +def test_treeshapiq_representation_limit_boundary_is_pinned(depth, is_supported): + """TreeSHAPIQ's boundary is lower than LinearTreeSHAP's 29/30. + + The identity N matrix (rhs of ones, built for every index) crosses the + magnitude limit at degree 26 already, so it is the binding constraint for + all TreeSHAPIQ indices; this pins the "~25" quoted in the error message. + """ + tree = _deep_chain_tree(depth) + if is_supported: + TreeSHAPIQ(model=tree, max_order=2, index="SII") + else: + with pytest.raises(RepresentationLimitError): + TreeSHAPIQ(model=tree, max_order=2, index="SII") + + +@pytest.mark.parametrize("depth", [35, 60]) +def test_too_deep_trees_are_refused_with_a_clear_error(depth): + """Beyond the float64 representation limit the explainer refuses loudly. + + Depth 35 previously returned silently wrong values (the ~32-59 band where the + explicit inverse is rank-deficient to machine precision); depth 60 crashed + with an unexplained ``LinAlgError``. The coefficients are now exact at any + depth, but their magnitude exceeds what the downstream float64 pipeline can + evaluate without destroying the result, so construction raises an + explanatory error instead. + """ + tree = _deep_chain_tree(depth) + with pytest.raises(RepresentationLimitError, match=r"interpolation degree.*too large"): + LinearTreeSHAP(tree) + + +def test_clear_solver_cache_forces_recomputation(): + """clear_solver_cache drops memoized solutions; a fresh solve still matches the oracle.""" + D = chebpts2(10) + rhs = np.ones(10) + solve_vandermonde(D, rhs) + assert _solve_vandermonde_cached.cache_info().currsize > 0 + clear_solver_cache() + assert _solve_vandermonde_cached.cache_info().currsize == 0 + np.testing.assert_array_equal(solve_vandermonde(D, rhs), _solve_exact_reference(D, rhs)) + + +def test_solver_results_are_cached_and_isolated(): + """Identical (points, rhs) solves hit the cache, and callers cannot corrupt it.""" + D = chebpts2(14) + rhs = np.ones(14) + first = solve_vandermonde(D, rhs) + first[0] = 12345.0 # mutate the returned copy + second = solve_vandermonde(D, rhs) + assert second[0] != 12345.0 # the cached value must be unaffected + np.testing.assert_array_equal(second, _solve_exact_reference(D, rhs)) + + +def test_deep_narrow_tree_routes_to_treeshapiq(): + """A deep tree over few features is re-routed instead of refused. + + LinearTreeSHAP's interpolation degree is the tree depth (above the limit + here), but TreeSHAPIQ's is min(depth, features in the tree) = 5, so + ``TreeExplainer`` falls back to TreeSHAPIQ and still produces order-1 + Shapley values satisfying completeness. + """ + depth = 40 + with pytest.raises(RepresentationLimitError): + LinearTreeSHAP(_deep_chain_tree(depth, n_features=5)) + + tree = _deep_chain_tree(depth, n_features=5) + explainer = TreeExplainer(model=[tree], index="SV", max_order=1) + x = np.zeros(depth) + explanation = explainer.explain(x) + prediction = tree.predict_one(x) + gap = abs(float(np.sum(explanation.values)) - prediction) + assert gap < 1e-9 * max(1.0, abs(prediction)) + + +def test_deep_narrow_tree_is_rerouted_for_sii_as_well(): + """The re-route covers both order-1 fast-path indices, not just ``"SV"``. + + ``index="SII"`` with ``max_order=1`` also takes the LinearTreeSHAP fast path; + the explainer validation normalizes it to ``"SV"`` (with a warning), which is + exactly the label the re-route assigns, so the re-routed result satisfies the + same completeness. + """ + depth = 40 + tree = _deep_chain_tree(depth, n_features=5) + with pytest.warns(UserWarning, match="SII generalizes 'SV'"): + explainer = TreeExplainer(model=[tree], index="SII", max_order=1) + assert len(explainer._treeshapiq_explainers) == 1 + x = np.zeros(depth) + explanation = explainer.explain(x) + prediction = tree.predict_one(x) + gap = abs(float(np.sum(explanation.values)) - prediction) + assert gap < 1e-9 * max(1.0, abs(prediction)) + + +def test_treeshapiq_n_matrices_are_gated_too(): + """The representation limit also fires inside TreeSHAPIQ's own N matrices. + + For higher-order indices there is no re-route: a tree whose feature-bounded + interpolation degree exceeds the limit is refused at construction instead of + returning silently wrong interaction values. + """ + tree = _deep_chain_tree(35) # 35 distinct features -> TreeSHAPIQ degree 35 + with pytest.raises(RepresentationLimitError): + TreeSHAPIQ(model=tree, max_order=2, index="k-SII") + + +def test_mixed_ensemble_keeps_the_fast_path_for_shallow_trees(): + """Only the over-deep tree leaves the LinearTreeSHAP fast path in an ensemble. + + The shallow tree stays on LinearTreeSHAP, the deep-narrow tree is re-routed to + TreeSHAPIQ, and the mixed aggregation still satisfies completeness for the + two-tree ensemble. + """ + assert issubclass(RepresentationLimitError, ValueError) + depth = 40 + deep_tree = _deep_chain_tree(depth, n_features=5) + shallow_tree = _deep_chain_tree(8) + explainer = TreeExplainer(model=[shallow_tree, deep_tree], index="SV", max_order=1) + assert len(explainer._lineartreeshap_explainers) == 1 + assert len(explainer._treeshapiq_explainers) == 1 + + x = np.zeros(depth) + explanation = explainer.explain(x) + prediction = shallow_tree.predict_one(x) + deep_tree.predict_one(x) + gap = abs(float(np.sum(explanation.values)) - prediction) + assert gap < 1e-9 * max(1.0, abs(prediction)) From 73618de6bd7488feafa7ecd6398c8ae5e966aebc Mon Sep 17 00:00:00 2001 From: YUANLUO WU Date: Thu, 11 Jun 2026 12:02:29 +0200 Subject: [PATCH 3/3] chore: apply the ruff autofix in test_public_api (split -> rsplit) --- tests/shapiq/tests_unit/test_public_api.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/shapiq/tests_unit/test_public_api.py b/tests/shapiq/tests_unit/test_public_api.py index f7103e13a..8f6cd378d 100644 --- a/tests/shapiq/tests_unit/test_public_api.py +++ b/tests/shapiq/tests_unit/test_public_api.py @@ -78,7 +78,7 @@ def test_all_concrete_subclasses_in_all(module_path: str, base_path: str, label: exported = set(module.__all__) missing = concrete - exported - pkg_init = f"src/shapiq/{module_path.split('.')[-1]}/__init__.py" + pkg_init = f"src/shapiq/{module_path.rsplit('.', maxsplit=1)[-1]}/__init__.py" assert not missing, ( f"Concrete {label} subclasses not listed in {module_path}.__all__: " f"{missing}. Add them to {pkg_init}."