Implement LeverageSHAP approximator

.gitignore

@@ -181,6 +181,7 @@ illustration_sources/
 src/shapiq_games/datasets/data/tabarena_*.csv
 
 shapiq/games/benchmark/precomputed/
+benchmark/results/*
 precomputed.zip
 game_storage/*
 

docs/source/references.bib

@@ -39,6 +39,13 @@ @inproceedings{Witter.2025
   year      = {2025},
   url       = {https://openreview.net/forum?id=Qabko39AS5}
 }
+@inproceedings{Musco.2025,
+  title     = {Provably Accurate Shapley Value Estimation via Leverage Score Sampling},
+  author    = {Christopher Musco and R. Teal Witter},
+  booktitle = {Proceedings of the International Conference on Learning Representations {(ICLR)}},
+  year      = {2025},
+  url       = {https://arxiv.org/abs/2410.01917}
+}
 @article{Arya.2020,
   title   = {{AI} Explainability 360: An Extensible Toolkit for Understanding Data and Machine Learning Models},
   author  = {Vijay Arya and Rachel K. E. Bellamy and Pin{-}Yu Chen and Amit Dhurandhar and Michael Hind and Samuel C. Hoffman and Stephanie Houde and Q. Vera Liao and Ronny Luss and Aleksandra Mojsilovic and Sami Mourad and Pablo Pedemonte and Ramya Raghavendra and John T. Richards and Prasanna Sattigeri and Karthikeyan Shanmugam and Moninder Singh and Kush R. Varshney and Dennis Wei and Yunfeng Zhang},

src/shapiq/__init__.py

@@ -20,6 +20,7 @@
     InconsistentKernelSHAPIQ,
     KernelSHAP,
     KernelSHAPIQ,
+    LeverageSHAP,
     OddSHAP,
     OwenSamplingSV,
     PermutationSamplingSII,
@@ -103,6 +104,7 @@
     "StratifiedSamplingSV",
     "OwenSamplingSV",
     "KernelSHAP",
+    "LeverageSHAP",
     "RegressionFSII",
     "RegressionFBII",
     "KernelSHAPIQ",

src/shapiq/approximator/__init__.py

@@ -16,6 +16,7 @@
     InconsistentKernelSHAPIQ,
     KernelSHAP,
     KernelSHAPIQ,
+    LeverageSHAP,
     OddSHAP,
     RegressionFBII,
     RegressionFSII,
@@ -58,6 +59,7 @@ def __init__(self, *_args: object, **_kwargs: object) -> None:
     UnbiasedKernelSHAP,
     PermutationSamplingSV,
     KernelSHAP,
+    LeverageSHAP,
     kADDSHAP,
     SPEX,
     RegressionMSR,
@@ -128,6 +130,7 @@ def __init__(self, *_args: object, **_kwargs: object) -> None:
     "StratifiedSamplingSV",
     "OwenSamplingSV",
     "KernelSHAP",
+    "LeverageSHAP",
     "RegressionFSII",
     "RegressionFBII",
     "KernelSHAPIQ",

src/shapiq/approximator/regression/__init__.py

@@ -5,6 +5,7 @@
 from .kadd_shap import kADDSHAP
 from .kernelshap import KernelSHAP
 from .kernelshapiq import InconsistentKernelSHAPIQ, KernelSHAPIQ
+from .leverageshap import LeverageSHAP
 from .oddshap import OddSHAP
 
 __all__ = [
@@ -13,6 +14,7 @@
     "KernelSHAP",
     "KernelSHAPIQ",
     "InconsistentKernelSHAPIQ",
+    "LeverageSHAP",
     "Regression",
     "RegressionFBII",
     "OddSHAP",

src/shapiq/approximator/regression/base.py

@@ -289,6 +289,7 @@ def kernel_shap_iq_routine(
                     X=regression_matrix,
                     y=residual_game_values[interaction_size],
                     kernel_weights=regression_weights,
+                    use_svd=False,
                 )
             else:
                 # for order > 2 use ground truth weights for sizes < interaction_size and > n -
@@ -314,6 +315,7 @@ def kernel_shap_iq_routine(
                     X=regression_matrix,
                     y=game_values_plus,
                     kernel_weights=regression_weights,
+                    use_svd=False,
                 )
 
                 sii_values_current_size = (
@@ -372,6 +374,7 @@ def regression_routine(
             X=regression_matrix,
             y=regression_response,
             kernel_weights=regression_weights,
+            use_svd=False,
         )
 
         if index_approximation in ["kADD-SHAP", "FBII"]:
@@ -638,31 +641,44 @@ def _get_regression_matrices(
     return regression_matrix, regression_weights
 
 
-def solve_regression(X: np.ndarray, y: np.ndarray, kernel_weights: FloatVector) -> np.ndarray:
-    """Solves the Shapley regression problem.
+def solve_regression(
+    X: np.ndarray,
+    y: np.ndarray,
+    kernel_weights: FloatVector,
+    *,
+    use_svd: bool = False,
+) -> np.ndarray:
+    """Solves the Shapley regression problem using weighted least squares (WLS).
 
-    Solves the regression problem using the weighted least squares method. Returns all approximated
-    interactions.
+    By default, this attempts a fast solution using the normal equations. If the
+    Gram matrix is singular or ill-conditioned, it falls back to a robust
+    Singular Value Decomposition (SVD) solver.
 
     Args:
         X: The regression matrix of shape ``[n_coalitions, n_interactions]``.
         y: The response vector for each coalition of shape ``[n_coalitions]``.
-        kernel_weights: The weights for the regression problem for each coalition.
+        kernel_weights: The weights for the regression problem of shape ``[n_coalitions]``.
+        use_svd: If ``True``, skips the fast normal equation solver and directly uses
+            the robust SVD-based least squares solver (``np.linalg.lstsq``). Useful
+            for cases with extreme weight initializations or known rank-deficiencies.
 
     Returns:
-        The solution to the regression problem.
-
+        The approximated interaction values of shape ``[n_interactions]``.
     """
+    # Explicit override: go straight to the robust, SVD-backed solver
+    if use_svd:
+        W_sqrt = np.sqrt(kernel_weights)
+        return np.linalg.lstsq(W_sqrt[:, np.newaxis] * X, W_sqrt * y, rcond=None)[0]
+
+    # Standard fast path (try the fast way, catch the error if it fails)
     try:
-        # try solving via solve function
         WX = kernel_weights[:, np.newaxis] * X
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=RuntimeWarning)
-            phi = np.linalg.solve(X.T @ WX, WX.T @ y)
-    except np.linalg.LinAlgError:
-        # solve WLSQ via lstsq function and throw warning
+            # Solves (X^T * W * X) * phi = X^T * W * y
+            return np.linalg.solve(X.T @ WX, WX.T @ y)
+
+    except (np.linalg.LinAlgError, ValueError):
+        # Fallback: Gram matrix is singular. Use robust SVD approach.
         W_sqrt = np.sqrt(kernel_weights)
-        X = W_sqrt[:, np.newaxis] * X
-        y = W_sqrt * y
-        phi = np.linalg.lstsq(X, y, rcond=None)[0]
-    return phi
+        return np.linalg.lstsq(W_sqrt[:, np.newaxis] * X, W_sqrt * y, rcond=None)[0]