Technical Summary: Fast Approximate Estimation of Conditional Shapley Values with Linear Explainers

Problem Statement

The paper addresses the computational bottleneck associated with estimating conditional Shapley values for linear regression models when covariates are dependent. While Shapley values provide a game-theoretic breakdown of feature contributions to a specific prediction, calculating them conditionally requires estimating the expected model output for all possible subsets (coalitions) of features.

For a model with $p$ predictors, there are $2^p$ possible coalitions. Existing methods, such as the Sequential and Iterative approaches in the shapr R-package, face significant challenges:

1. Sequential Estimation: Requires fitting $2^p$ separate regression models, which becomes computationally intractable as $p$ increases.
2. Iterative Estimation: Samples coalitions until a convergence criterion is met. While faster when few coalitions suffice, it can still require a large fraction of all $2^p$ coalitions to achieve convergence, leading to long runtimes (hours) comparable to the sequential method in difficult cases.
3. Feature Dependence: Unlike marginal Shapley values, conditional values must account for the dependence structure among covariates, necessitating complex conditional expectation calculations.

The authors aim to develop methods that can estimate the coefficients for all $2^p$ submodels simultaneously and rapidly, leveraging the sparsity of the underlying mathematical structures.

Methodology

The proposed approach utilizes Constrained Gaussian Markov Random Field (GMRF) theory and sparse matrix algebra to estimate regression coefficients for all submodels jointly rather than sequentially. The authors derive three distinct algorithms:

1. The Diagonal Correction Method (Approximate)

This method approximates the constrained regression by modifying the precision matrix ($Q = X^T X$) of the full model.

• Mechanism: To enforce a constraint where a specific coefficient is zero (i.e., the feature is excluded from a submodel), a large scalar value $\kappa$ is added to the corresponding diagonal element of the precision matrix.
• Theoretical Basis: Based on the Woodbury matrix identity, adding $\kappa A^T A$ (where $A$ is a constraint matrix) to $Q$ approximates the conditional covariance matrix as $\kappa \to \infty$.
• Implementation: A joint constraint matrix $A$ is constructed for all $2^p$ models. The method solves the linear system $(Q_F + \kappa A^T A)^{-1}m$, where $Q_F$ is a block-diagonal matrix containing the full model's precision matrix repeated $2^p$ times.
• Tuning: Requires a tuning parameter $\kappa$, suggested to be $10^5$ times the largest eigenvalue of the original precision matrix.

2. The Projection Method (Approximate)

This method projects the parameter space to remove constrained parameters.

• Mechanism: A projection matrix $Z = I - A^T A$ is applied to the precision matrix. The constrained model is treated as an intrinsic GMRF. To ensure numerical stability (positive definiteness), a small scalar $\epsilon$ is added to the diagonal of the projected precision matrix ($Z Q Z + \epsilon I$).
• Theoretical Basis: As $\epsilon \to 0^+$, the solution converges to the exact constrained estimates.
• Implementation: Solves the system $Z(Z Q_F Z + \epsilon I)^{-1} Z m$.
• Tuning: Requires a small tuning parameter $\epsilon$ (e.g., $10^{-5}$), chosen to be small relative to the smallest eigenvalue of $Q$ but large enough for stable Cholesky decomposition.

3. The Exact Transformation Method (Exact)

This method provides an exact solution without approximation errors or tuning parameters.

• Mechanism: Instead of working in the full parameter space and forcing coefficients to zero, the method transforms the problem into a lower-dimensional space containing only the non-zero parameters.
• Implementation: A mapping matrix $E$ is constructed to select only the active parameters for each submodel. The system is solved as $(E Q_F E^T)^{-1} E m$.
• Advantage: Eliminates the need for large penalty terms ($\kappa$) or small regularization terms ($\epsilon$), avoiding potential numerical instability or approximation bias.

All three methods utilize sparse matrix algebra (specifically Cholesky decomposition) to handle the joint estimation of thousands of models efficiently. The joint precision matrix is highly sparse because the models are independent (block-diagonal structure), and off-diagonal zeros indicate conditional independence.

Key Contributions

1. Three New Algorithms: The paper introduces two approximate methods (Diagonal Correction, Projection) and one exact method (Exact Transformation) for the joint estimation of all submodels in linear regression.
2. Convergence Proofs: The authors provide mathematical proofs demonstrating that the approximate methods converge to the true conditional estimates as their respective tuning parameters ($\kappa \to \infty$ and $\epsilon \to 0^+$).
3. Computational Efficiency: By exploiting sparsity and joint estimation, the methods reduce computation time from hours (typical of shapr) to seconds or minutes, even when enumerating all $2^p$ coalitions.
4. Empirical Validation: Extensive numerical studies compare the new methods against the shapr package (both Sequential and Iterative versions) across three datasets: Adult Income, a Simulated Gaussian dataset, and the WHO Life Expectancy dataset.

Results

The numerical case studies yield the following findings:

• Accuracy: The Shapley value estimates produced by the new methods are virtually identical to the Sequential method in shapr (used as the ground truth for full enumeration).
  • In the Adult dataset ($p=11$), the median absolute error (MedAE) between the new methods and the Sequential method was on the order of $10^{-9}$.
  • In the Simulated dataset ($p=21$), MedAE was on the order of $10^{-7}$ to $10^{-5}$.
  • In the WHO dataset ($p=16$), MedAE was on the order of $10^{-6}$.
• Speed:
  • Adult Dataset: The new methods ran in 2.5 to 9 seconds, whereas shapr (Sequential and Iterative) took 17 to 19 minutes.
  • Simulated Dataset: The new methods took 1.4 to 2.5 minutes to estimate all $2^{21} \approx 2.1$ million models. The shapr Iterative method was faster (6.48 seconds) but only estimated 200 coalitions, not the full set.
  • WHO Dataset: The new methods took 1.1 to 1.2 minutes. The shapr Sequential method took 131 minutes (2.2 hours), and the Iterative method took 16 minutes (parallel) or 36 minutes (serial).
• Scalability: The methods successfully estimated over two million submodels in the simulated study within minutes. The primary limitation identified is memory constraints (Curse of Dimensionality) when $p$ exceeds ~21 on standard hardware, though the authors note this can be mitigated with more RAM or specialized sparse matrix packages.

Significance and Claims

The paper claims that the proposed methods offer a substantial improvement in computational efficiency for conditional Shapley value estimation without sacrificing accuracy.

• Superiority in Full Enumeration: When the shapr Iterative method requires a large fraction of coalitions to converge (as seen in the Adult and WHO datasets), the new methods are "substantially faster" while maintaining "equal or better accuracy."
• Comprehensive Enumeration: Unlike the iterative approach which samples coalitions, the new methods provide a full enumeration of all $2^p$ submodels, ensuring no information is lost due to sampling variance.
• Recommendation: The authors recommend the Exact Transformation Method for practical use. While all three methods yield comparable results, the Exact method has no tuning parameters, thereby removing the need for hyperparameter selection and potential sources of approximation error.
• Future Potential: The authors suggest that the methodology can be extended to handle splines (e.g., from the rms package) and that parallelization (e.g., via OpenMP) could further reduce runtimes, though the current implementation relies on serial execution for fair comparison with the parallelized shapr package.

The paper concludes that these methods make the exact calculation of conditional Shapley values for linear models with mixed covariate types feasible in scenarios where it was previously computationally prohibitive.