arXiv CSopen access

Probabilistic Linear Explanations

Probabilistic Linear Explanations⋆ Frédéric Korichea,∗,1 , Jean-Marie Lagnieza,2 and Chi Trana,2 a Computer science Research Institute of Lens (CRIL), UMR CNRS 8188, University of Artois, Rue Jean Souvraz SP 18, Lens, F-62307, France

ABSTRACT

Keywords: Explainable Artificial Intelligence Probabilistic Explanations Sparse Linear Functions Iterative Hard Thresholding Mixed Integer Programming

Formal explainability provides mathematically grounded justifications for individual predictions. However, abductive explanations often exceed human cognitive limits by involving too many features, while probabilistic relaxations have remained largely limited to categorical classification. We present a unified framework for probabilistic explainability based on sparse, anchored linear models, applicable to both binary classification and continuous regression. By mapping instances to the Boolean hypercube, our linear explanations strictly generalize subset-based approaches: they capture both the magnitude and direction of feature contributions while enforcing a prescribed sparsity budget 𝑘. We show that minimizing the relevance error for such explanations is NPPP -hard when the underlying model is a neural network, and we relate this intractable objective to a tractable surrogate—the fidelity error. For a parameterized family of local distributions, the relevance error of any 𝑘-sparse explanation is bounded by its fidelity error up to a multiplicative factor that remains small locally. We address the resulting empirical problem using two complementary approaches: a Mixed Integer Programming (MIP) formulation that yields provably optimal empirical solutions while maintaining polynomial sample complexity, and a polynomial-time Iterative Hard Thresholding (IHT) algorithm with provable approximation guarantees. Empirical evaluations show that, unlike state-of-the-art baselines such as LIME and MAPLE, our explanations satisfy both the anchoring and sparsity constraints by construction, while consistently achieving lower relevance error.

arXiv:2609.19077v1 [cs.LG] 16 Sep 2026

ARTICLE INFO

1. Introduction

As machine learning models increasingly influence high-stakes domains—such as criminal justice, medical diagnosis, and social scoring—the need for ethics, fairness, and safety has become more pressing than ever. Explainable Artificial Intelligence (XAI) confronts this challenge by providing methods that enable users to interpret model behavior without requiring deep technical expertise (Burkart and Huber, 2021; Molnar, 2025). More recently, the subfield of formal explainability has emerged to offer rigorous mathematical guarantees regarding explanation quality, size, and semantics. The primary aim is to establish robust theoretical foundations for explaining predictions, thereby systematically fostering trust and confidence in model capabilities (Ignatiev, 2020; Marques-Silva and Ignatiev, 2022). When a machine learning model makes a consequential decision, the most natural question is, “Why?” Formal explainability seeks to translate this question into a clear, interpretable answer with mathematical guarantees. For a given data instance 𝒙 and classifier 𝑓 , the aim is to identify a rule that justifies the output 𝑓 (𝒙). Such a rule can be represented by a feature subset 𝑆, so that changing features outside 𝑆 does not alter the outcome 𝑓 (𝒙). The restriction of 𝒙 to 𝑆, denoted 𝒙𝑆 , thus contains all necessary information to determine 𝑓 (𝒙). Accordingly, 𝑆 is often called a (weak) abductive explanation (Cooper and Marques-Silva, 2023; Ignatiev, Narodytska and Marques-Silva, 2019) or a sufficient reason (Barceló, Monet, Pérez and Subercaseaux, 2020; Darwiche and Hirth, 2020). Despite their logical rigor, abductive explanations are often too complex for practical use (Ignatiev, Cooper, Siala, Hebrard and Marques-Silva, 2020; Ignatiev et al., 2019). This limitation fundamentally stems from human cognitive constraints. Classic psychological studies show that while people can retain about seven distinct pieces of information in short-term memory (Miller, 1956), our ability to process interacting variables is even more limited—typically capped at four or five elements before cognitive overload (Halford, Wilson and Phillips, 1998; Johnson-Laird, 2010). As a result, empirical XAI research emphasizes the need for concise explanations to achieve true interpretability (Lage, Chen, He, Narayanan, Kim, Gershman and Doshi-Velez, 2019; Narayanan, Chen, He, Kim, Gershman and Doshi-Velez, 2018). ⋆

This manuscript is an extended version of a paper presented at the 41st Conference on Uncertainty in Artificial Intelligence (UAI 2025), which received a Best Paper Award. ∗ Corresponding author [email protected] (F. Koriche); [email protected] (J. Lagniez); [email protected] (C. Tran) ORCID (s):

F. Koriche et al.: Preprint submitted to Elsevier

Page 1 of 35

Probabilistic Linear Explanations

To overcome these cognitive bottlenecks, recent research has shifted toward probabilistic explanations, which seek to balance certainty with interpretability (Blanc, Lange and Tan, 2021; Izza, Huang, Ignatiev, Narodytska, Cooper and Marques-Silva, 2023; Wäldchen, MacDonald, Hauch and Kutyniok, 2021). In this framework, a probabilistic explanation for a classifier 𝑓 and instance 𝒙 is a small feature subset 𝑆 such that 𝒙𝑆 determines 𝑓 (𝒙) with high probability. The effectiveness of 𝑆 is measured by its relevance error: the probability that 𝑓 assigns a different label to a random instance 𝒛 than to 𝒙, even when 𝒛 and 𝒙 are indistinguishable on 𝑆. Evaluated with respect to a distribution — such as the uniform distribution or a local neighborhood—finding a probabilistic explanation becomes a constrained stochastic optimization problem. For instance, seeking an explanation with the lowest relevance error, subject to a user-defined size bound 𝑘, leads to the following formulation: minimize

ℙ𝒛∼ [𝑓 (𝒛) ≠ 𝑓 (𝒙) ∣ 𝒛𝑆 = 𝒙𝑆 ]

subject to

|𝑆| ≤ 𝑘

(P1)

Example 1. For illustration, consider a binary classifier 𝑓 that evaluates loan applications, where the applicant’s profile 𝒙 is represented by a set of discrete, Boolean propositions. Suppose applicant 𝒙 is approved (𝑓 (𝒙) = +1), even though their credit history is shorter than 10 years. A strict abductive explanation might require fixing several features—such as [Account Age ≥ 5 yrs], [Citizen = Yes], [Income ≥ 65k], [Debt-to-Income Ratio (DTI) ≤ 30%], [Credit History < 10 yrs], and [Employed ≥ 2 yrs]—to mathematically guarantee that 𝑓 outputs +1. However, presenting such a lengthy rule risks the very cognitive overload that XAI seeks to prevent. By instead framing the explanation as the stochastic optimization problem in (P1) and imposing a cognitive size limit of 𝑘 = 2, we might identify a smaller subset 𝑆 = {[Income ≥ 65k], [DTI ≤ 30%]}. In this case, the relevance error is the probability that a random applicant 𝒛 ∼  matching 𝒙 on these two specific conditions would be denied the loan (𝑓 (𝒛) = −1). If this error is 0.02, then among all applicants agreeing with 𝒙 on these two conditions, only 2% receive a different decision. The user receives a clear, digestible answer, while the probabilistic framework rigorously quantifies the rare edge cases where the explanation does not hold. While probabilistic explanations successfully address cognitive limitations, their use has largely been restricted to categorical classification. This prompts a key question: how can these guarantees be extended to a broader class of predictive models? In this paper, we present a unified framework that represents data instances as vectors in the Boolean hypercube {−1, +1}𝑑 . Continuous or categorical attributes are mapped to Boolean features in advance, following standard practice in model-agnostic explainability (Ribeiro, Singh and Guestrin, 2016). Our approach encompasses both binary classification models—where 𝑓 takes values in {−1, +1}—and continuous regression models, where 𝑓 spans an interval such as [−1, +1].1 Importantly, we impose no assumptions on the internal structure of 𝑓 : whether 𝑓 is a tree ensemble, support vector machine, or deep neural network, it is treated entirely as a black box, accessible only through value queries. When explaining a continuous prediction 𝑓 (𝒙), selecting a feature subset 𝑆 alone often fails to capture both the magnitude and direction of each feature’s influence. 𝑆 acts as a binary mask—a support of a weight vector that indicates only which features are present, without conveying the strength or sign of their contributions. To address this shortcoming, we define explanations as sparse linear functions 𝒘 ∈ ℝ𝑑 that satisfy the anchoring hyperplane condition 𝒘 ⋅ 𝒙 = 𝑓 (𝒙). This constraint, related to the local accuracy property (Lundberg and Lee, 2017), ensures that any linear explanation 𝒘 remains perfectly consistent with the black-box model at 𝒙. The sparsity of 𝒘—its number of nonzero coefficients, ‖𝒘‖0 —reflects its conciseness. By allowing weights to take continuous values, 𝒘 captures both the direction and magnitude of each feature’s contribution, thereby improving interpretability and expressiveness. Crucially, sparse linear explanations provide a strict generalization of subset-based explanations. Since 𝒘 is zero outside its support 𝑆, any instance 𝒛 that matches 𝒙 on 𝑆 will also satisfy 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙, so the set of instances described by 𝒘 always contains those described by 𝑆. If the coefficients of 𝒘 are strictly non-compensatory (that is, no combination of feature weights cancels another), the two sets coincide and we recover the original subset explanation. Otherwise, 𝒘 describes a broader set of instances and encodes more expressive rules: instead of requiring every feature in 𝑆 to be fixed, it only constrains a weighted combination of features, allowing the explanation to capture compensatory effects between conditions. To assess how well a sparse linear model approximates the black-box model, we replace the conditional zero-one loss in (P1) with a normalized conditional quadratic loss: 𝓁(𝑦, 𝑦′ ) = 14 (𝑦 − 𝑦′ )2 . In our unified framework, this loss is 1 Standardizing the codomain to [−1, +1] simplifies notation and highlights the mathematical parallel with the zero-one loss in binary classification. All theoretical results extend to any bounded continuous range [−𝑐, +𝑐] for constant 𝑐 > 0 by adjusting the normalization factor.

F. Koriche et al.: Preprint submitted to Elsevier

Page 2 of 35

Probabilistic Linear Explanations

applied consistently to both continuous regression and binary classification tasks. The factor of 1∕4 ensures that for binary labels {−1, +1}, the quadratic loss matches the zero-one classification error exactly. The relevance error for an explanation 𝒘 is defined as the conditional expectation of 𝓁(𝑓 (𝒛), 𝑓 (𝒙)), given 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙. With respect to the chosen distribution , the central optimization problem addressed in this paper is: minimize

𝔼𝒛∼ [𝓁(𝑓 (𝒛), 𝑓 (𝒙)) ∣ 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙]

subject to

𝒘 ⋅ 𝒙 = 𝑓 (𝒙)

and

‖𝒘‖0 ≤ 𝑘

(P2)

Example 2. Consider now the same applicant in a continuous regression task: assigning a personalized loan interest rate, where lower values are preferable. Suppose 𝒙 receives a rate of 8% (𝑓 (𝒙) = 0.08). Reporting the feature subset 𝑆 = {[DTI ≤ 30%], [Credit History ≥ 10 yrs]} shows which variables were important, but not how they influenced the outcome. Solving (P2) with 𝑘 = 2 yields 𝑤[DTI≤30%] = −0.03 and 𝑤[Credit History≥10 yrs] = −0.11. Since the contribution of a feature is 𝑤𝑖 𝑥𝑖 , and 𝒙 satisfies the first condition (𝑥𝑖 = +1) but not the second (𝑥𝑖 = −1), the rate decomposes as −0.03 + 0.11 = 0.08. The user thus learns that the low debt ratio contributes −3 points to the rate while the short credit history contributes +11, the two effects partially compensating—an insight that no feature subset can convey. While (P2) provides a unified framework for generating linear explanations, solving it exactly presents severe computational hurdles. In Section 3, we demonstrate that for neural networks, (P2) is NPPP -hard—placing it well beyond the capabilities of modern exact solvers. However, this hardness does not preclude the existence of algorithms with approximation guarantees relative to a tractable “surrogate” error. We focus on the fidelity error, a widely used metric in model-agnostic explainability (Lakkaraju, Kamar, Caruana and Leskovec, 2019; Li, Nagarajan, Plumb and Talwalkar, 2021; Ribeiro et al., 2016), which measures the normalized quadratic loss incurred by 𝒘 on random instances labeled by 𝑓 . By substituting the conditional objective in (P2) with the fidelity error, we obtain: minimize

𝔼𝒛∼ [𝓁(𝒘 ⋅ 𝒛, 𝑓 (𝒛))]

subject to

𝒘 ⋅ 𝒙 = 𝑓 (𝒙)

and

‖𝒘‖0 ≤ 𝑘

(P3)

The connection between (P2) and (P3) relies on the choice of probability distribution . Drawing on distance-based models for discrete spaces, such as the Mallows model (Mallows, 1957) and its extensions (Fligner and Verducci, 1993; Marden, 1996), we focus on the following family of distributions: ) ( (𝒛) ∝ exp − 𝜎2 ‖𝒛 − 𝒙‖1 Here, 21 ‖𝒛 − 𝒙‖1 is the Hamming distance (the number of features on which 𝒙 and 𝒛 differ), and 𝜎 is a concentration parameter controlling how tightly instances 𝒛 cluster around 𝒙. Notably, setting 𝜎 = 0 recovers the uniform distribution over the instance space. Throughout, 𝜎 is a user-supplied parameter that specifies how local the explanation is meant to be; Section 6 reports how our results vary with it. In Section 4, we prove that for this parameterized family, the relevance error of any feasible solution to (P2) is at most (1 + 𝑒−𝜎 )𝑘 times its fidelity error. Furthermore, since (P3) involves a non-conditional expected loss function, it is highly amenable to empirical sampling. Thus, by approximating the expected fidelity error over 𝑚 samples, we arrive at the following practical, data-driven optimization problem: minimize

1 𝑚

𝑚 ∑

𝓁(𝒘 ⋅ 𝒛𝑖 , 𝑓 (𝒛𝑖 ))

(P4)

𝑖=1

subject to

𝒘 ⋅ 𝒙 = 𝑓 (𝒙)

and

‖𝒘‖0 ≤ 𝑘

This task generalizes the classic subset selection (Beale, Kendall and Mann, 1967; Hocking and Leslie, 1967) problem, also known as sparse regression (Natarajan, 1995), and is therefore NP-hard. It can nonetheless be solved exactly by Mixed Integer Programming (MIP), which remains tractable on instances of moderate size. We further show that a sample size 𝑚 polynomial in 𝑘, 1∕𝜀 and log 𝑑 suffices: with high probability, the resulting 𝑘-sparse explanation has relevance error at most (1+𝑒−𝜎 )𝑘 (𝑣∗ +𝜀), where 𝑣∗ is the optimal value of (P4). The quality of MIP-based explanations therefore improves as the concentration 𝜎 increases and, notably, the sample complexity does not depend on 𝜎 at all. As a strictly polynomial-time alternative, Section 5 adapts the Iterative Hard Thresholding (IHT) algorithm (Blumensath and Davies, 2008, 2009) to compute approximate solutions to (P4). Our variant projects, at each iteration, F. Koriche et al.: Preprint submitted to Elsevier

Page 3 of 35

Probabilistic Linear Explanations

onto the intersection of the 𝑘-sparsity constraint and the anchoring hyperplane—an operation that is both exact and computable in (𝑑 log 𝑑 + 𝑘2 ) time. For our family of distributions, IHT produces 𝑘-sparse explanations whose relevance error is, with high probability, at most (1 + 𝑒−𝜎 )𝑘 (𝑐 ⋅ 𝑣∗ + 𝜀), where 𝑐 > 1 is an explicit multiplicative factor governed by 𝑘 and 𝜎. However, the sample complexity for IHT grows with cosh4 (𝜎∕2)—that is, as 𝑒2𝜎 for large 𝜎. This discrepancy with the MIP bound is not merely analytical: while the MIP guarantee only requires uniform convergence of the empirical fidelity error, the IHT guarantee also depends on restricted strong convexity and smoothness constants, which deteriorate as  becomes more concentrated around 𝒙. As a result, the parameter 𝜎 governs a meaningful computational and statistical trade-off. In the localized regime (large 𝜎), the relevance bound is tight and MIP is preferred, as its sample complexity remains independent of 𝜎. In contrast, IHT faces a double penalty: its sample complexity grows as 𝑒2𝜎 and its approximation factor worsens as 𝑒𝜎 . In the intermediate regime, where 𝜎 is large enough for the relevance bound to be informative but small enough to keep both the sampling cost and the approximation factor moderate, IHT becomes a compelling alternative, yielding highly relevant explanations in polynomial time. Section 6 empirically benchmarks our MIP and IHT methods against the LIME (Ribeiro et al., 2016) and MAPLE (Plumb, Molitor and Talwalkar, 2018) explainers on both classification and regression tasks. Our results show that MAPLE consistently exceeds the sparsity budget 𝑘 and frequently violates the anchoring condition as well, placing its explanations outside the intended interpretability regime. While LIME satisfies the sparsity budget, it systematically deviates from anchoring and thus fails to reproduce the black-box prediction 𝑓 (𝒙) at the instance 𝒙 being explained. Because LIME operates over a strictly larger hypothesis space, it sometimes achieves a lower empirical fidelity error; our measurements indicate that this apparent advantage results directly from the anchoring deviation. Nevertheless, across variations in 𝑘, 𝜎, and 𝑚, both IHT and MIP consistently outperform LIME in terms of relevance—the objective our framework is designed to optimize. Between our two methods, relevance rarely differs by more than a few hundredths, so the choice becomes essentially computational: MIP certifies optimality within seconds on small and medium benchmarks but loses this certificate on the largest ones, whereas IHT is by far the fastest explainer and remains effective at every scale. In summary, our main contributions are as follows: • We present a unified formal framework for probabilistic explainability, bridging categorical classification and continuous regression through sparse, anchored linear models that strictly generalize subset-based explanations. We also prove that computing the optimal relevance error for neural network black-box models is NPPP -hard. • We establish a theoretical link between the intractable relevance error and a tractable surrogate—the fidelity error. By introducing a parameterized family of local distributions, we show that the relevance error of a 𝑘-sparse explanation is bounded by (1 + 𝑒−𝜎 )𝑘 times its fidelity error. • We formulate the empirical fidelity problem, encode its exact solution as a MIP, and approximate it with a strictly polynomial-time IHT variant. Our experiments validate the theoretical guarantees, demonstrating that, in contrast to MAPLE and LIME, our methods consistently satisfy both sparsity and anchoring constraints while achieving superior relevance. Furthermore, they offer complementary practical advantages: MIP reliably finds optimal solutions for medium-dimensional problems, while IHT scales efficiently to high-dimensional datasets. The remainder of the paper is organized as follows. Section 2 reviews related work. Section 3 introduces sparse anchored linear explanations and establishes the hardness of (P2). Section 4 relates the relevance error to the fidelity error and derives the MIP formulation, while Section 5 presents our IHT variant and its approximation guarantee. Section 6 reports the experimental evaluation, and Section 7 concludes.

2. Related Work To contextualize our contributions, we organize the landscape of local explainability into two primary categories: probabilistic explanations, which emphasize rigorous theoretical guarantees, and linear explanations, which focus on flexibility and the ability to handle continuous outcomes.

2.1. Probabilistic Explanations Probabilistic explanations, as a natural generalization of abductive explanations, have become increasingly prominent in formal explainability. Whereas traditional abductive reasoning seeks a feature subset that logically entails F. Koriche et al.: Preprint submitted to Elsevier

Page 4 of 35

Probabilistic Linear Explanations

a prediction, probabilistic approaches relax this requirement to address human cognitive limits, aiming for subsets that guarantee the prediction only with high probability. In categorical classification, a feature subset 𝑆 is called a 𝑘-size 𝜏-relevant probabilistic explanation for a reference instance 𝒙, classifier 𝑓 , and distribution  if |𝑆| ≤ 𝑘 and its zero-one relevance error is at most 1 − 𝜏. Recall that for binary labels, this zero-one error coincides exactly with the normalized quadratic relevance error of (P2), so the two formulations agree on classification tasks. Within this framework and the setup of (P1), one can either fix a cognitive size limit 𝑘 and seek 𝑆 that minimizes relevance error (Bounia and Koriche, 2023; Koriche, Lagniez, Mengel and Tran, 2024), or set a relevance threshold 𝜏 and find the smallest 𝑆 achieving relevance error at most 1 − 𝜏 (Wäldchen et al., 2021; Izza et al., 2023). Since exactly computing these probabilities typically requires knowledge of the model’s internal architecture, both approaches are usually studied in a model-specific (white-box) setting. Determining whether a 𝑘-size 𝜏-relevant explanation exists for a given 𝒙 is computationally hard: the problem is NPPP -hard for neural networks (Wäldchen et al., 2021) and NP-hard for decision trees (Arenas, Barceló, Orth and Subercaseaux, 2022). Even at 𝜏 = 1—the case of abductive explanations—the task remains p2 -hard for random forests (Audemard, Bellart, Bounia, Koriche, Lagniez and Marquis, 2022) and neural networks (Barceló et al., 2020). The outlook for fixed-parameter tractability and approximability is similarly negative. Finding cardinality-minimal 𝜏relevant explanations is 𝖶[2]-hard in 𝑘, even for monotone neural networks and 𝜏 = 1 (Ordyniak, Paesani and Szeider, 2023), and is NP-hard to approximate within any factor of 𝑑 1−𝛼 for any 𝛼 > 0, for both neural networks (Wäldchen et al., 2021) and decision trees (Kozachinskiy, 2023). To the best of our knowledge, the main positive result to date applies to linear functions 𝑓 , where a relaxed problem becomes polynomial-time approximable (Subercaseaux, Arenas and Meel, 2025). All of these results concern subset-based explanations. Whether the same barriers apply to linear explanations is not immediate: the constraint 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙 selects a far richer family of regions than the subcube induced by a subset, so neither hardness nor tractability transfers mechanically. We settle this question in Section 3. Shifting to black-box classifiers, Blanc et al. (2021) showed that if user-supplied instances 𝒙 are drawn uniformly at random, a 𝜏-relevant explanation 𝑆 of size 𝑘 can be constructed, with high probability, from the decision path 𝑇 (𝒙) of a depth-𝑘 decision tree 𝑇 , provided that 𝑇 disagrees with 𝑓 on at most a 1 − 𝜏 fraction of the instances. If 𝑘 is polynomial in the “average certificate complexity” of 𝑓 , the tree 𝑇 can be learned in polynomial time. In our terminology, their result links fidelity to relevance: a surrogate with low fidelity error yields an explanation with low relevance error. Our bound in Section 4 has the same form, but uses a sparse linear surrogate instead of a tree. Notably, their guarantee is average-case (over random 𝒙), while ours is pointwise (for every fixed 𝒙), which is the setting local explainability calls for. Consequently, the two guarantees are complementary, and neither subsumes the other. Among model-agnostic approaches, A NCHORS (Ribeiro, Singh and Guestrin, 2018) deserves special mention. This method, designed for black-box classifiers, outputs a feature subset (or logical decision rule) 𝑆 whose precision— the probability that 𝑓 agrees with 𝑓 (𝒙) on instances satisfying 𝑆—exceeds a user-specified threshold, estimated via sampling rather than computed from a white-box model. Precision, being a conditional quantity, aligns with relevance rather than fidelity in our terminology. Thus, A NCHORS sits at the intersection of probabilistic and model-agnostic approaches—incorporating elements of both, but only heuristically. Notably, there are no formal guarantees on the size of the returned subset: the anchor simply expands until the precision threshold is met, and achieving that threshold is trivial, since including all features of 𝒙 guarantees a precision of 1. The main computational challenge—finding a short anchor with high precision—is delegated to a greedy beam search algorithm, which lacks any approximation guarantee. Most importantly for our purposes, A NCHORS returns a subset rather than a linear function, placing it outside the class of explanations considered in our work. Accordingly, our experimental comparisons in Section 6 focus exclusively on explainers that produce linear functions. A related but distinct line of work concerns contrastive (or counterfactual) explanations, which identify a subset of features that, when modified, alter the prediction; we refer the reader to Guidotti (2024) for a survey. These methods answer “what would have to change?”, whereas abductive explanations and their probabilistic extensions answer “why this prediction?”. The two questions have recently been brought together by Bassan, Huang and Katz (2026), who show that sufficient and contrastive reasons, in both their local and global forms, can all be characterized as minimizers of a single probabilistic value function, and who tie the complexity of computing them to monotonicity, submodularity, and supermodularity properties of that function. Their unification is orthogonal to ours: it ranges over explanation types while remaining within feature subsets and categorical classification, whereas we unify classification and regression tasks within the class of sparse linear functions. Because no prior probabilistic framework has encompassed continuous regression, we focus exclusively on this abductive trajectory and do not consider contrastive explanations further.

F. Koriche et al.: Preprint submitted to Elsevier

Page 5 of 35

Probabilistic Linear Explanations

2.2. Linear Explanations To address the limitations of discrete subsets—specifically their inability to explain continuous regression tasks— and to abstract away the internal architectures of target models, a separate line of research has focused on modelagnostic linear explanations. These methods extract an interpretable linear surrogate from the local neighborhood of a data instance (Ribeiro et al., 2016; Lundberg and Lee, 2017; Plumb et al., 2018; Agarwal, Jabbari, Agarwal, Upadhyay, Wu and Lakkaraju, 2021; Zhao, Huang, Huang, Robu and Flynn, 2021). A common goal across these approaches is to minimize an unconstrained objective of the form: 1 𝑚

𝑚 ∑

𝜙𝒙 (𝒛𝑖 )(𝑓 (𝒛𝑖 ) − 𝒘 ⋅ 𝒛𝑖 )2 + 𝜓(𝒘)

(1)

𝑖=1

Here, {(𝒛𝑖 , 𝑓 (𝒛𝑖 ))}𝑚 is a set of labeled samples generated from a neighborhood distribution around 𝒙, the weighting 𝑖=1 function 𝜙𝒙 (𝒛𝑖 ) assesses the importance of 𝒛𝑖 , and the regularization term 𝜓(𝒘) penalizes the complexity of the linear model 𝒘. For example, in the widely used LIME method (Ribeiro et al., 2016), 𝜙𝒙 (𝒛𝑖 ) is defined as a normalized distance between 𝒛𝑖 and 𝒙. Similarly, in K ERNELSHAP (Lundberg and Lee, 2017), 𝜙𝒙 (𝒛𝑖 ) is derived from combinatorial Shapley weights to ensure game-theoretic properties. In the MAPLE method (Plumb et al., 2018), 𝜙𝒙 (𝒛𝑖 ) is the fraction of trees, in a random forest trained on 𝑓 , for which 𝒛𝑖 falls in the same leaf as 𝒙. Despite their widespread adoption, the reliance on the unconstrained, heuristic objective in (1) introduces significant vulnerabilities. A growing body of literature highlights the local instability of these methods, demonstrating that semantically indistinguishable inputs can yield vastly different linear explanations (Alvarez-Melis and Jaakkola, 2018). Furthermore, because the sampling distribution and weighting functions are not explicitly tied to formal performance bounds, these explainers can be manipulated by measurement biases or adversarial perturbations (Ghorbani, Abid and Zou, 2019; Slack, Hilgard, Jia, Singh and Lakkaraju, 2020), and their theoretical convergence heavily depends on the choice of hyperparameters (Garreau and von Luxburg, 2020). Our framework departs from this template in three fundamental ways. First, we rigorously formalize locality, moving beyond heuristic approaches. The family of 𝜎-parameterized distributions introduced in Section 1 serves a role similar to LIME’s exponential kernel, but crucially, it acts as the formal probability law defining the expected loss—and thus the relevance error—instead of an arbitrary heuristic weight in an empirical sum. This distinction enables us to state formal bounds. Second, we impose consistency at the reference instance as a hard constraint: 𝒘 ⋅ 𝒙 = 𝑓 (𝒙). As noted in Section 1, this anchoring constraint is not new: it is the local accuracy property of Lundberg and Lee (2017), or ∑ equivalently, the efficiency axiom 𝜙0 + 𝑑𝑖=1 𝜙𝑖 = 𝑓 (𝒙) of Shapley-based explainers. Our innovation lies in enforcing this anchoring condition jointly with a strict sparsity budget. Third, we enforce sparsity as the hard constraint ‖𝒘‖0 ≤ 𝑘, rather than through a soft penalty 𝜓 or a relaxed selection procedure. While LIME targets a sparse support via LASSO selection of 𝑘 features, it provides no guarantee on the resulting support; MAPLE has no sparsity mechanism, as its local linear model assigns a weight to every feature. Solving (P3) under both constraints yields explanations with theoretically sound relevance and fidelity bounds. These distinctions clarify why K ERNELSHAP, though it produces linear explanations that satisfy the anchoring constraint, is not used as a baseline in Section 6. Its minimization of (1) is instrumental rather than predictive: the Shapley kernel is chosen so that the solution of the weighted least-squares program matches the Shapley values, not to ensure local surrogate accuracy. Furthermore, its samples correspond to coalition maskings rather than to neighbors of 𝒙. Most importantly, Shapley values are inherently dense—every feature receives a value—so K ERNELSHAP cannot satisfy a sparsity budget ‖𝒘‖0 ≤ 𝑘. Truncating to the 𝑘 largest values breaks the local accuracy property (the anchoring constraint). Any explainer that cannot be made feasible for (P3) without sacrificing its defining property is not a meaningful comparator.

3. Problem Formulation and Complexity To rigorously study the trade-offs between explanation sparsity, exactness, and computational tractability, we must first cast the intuitions from Section 1 into a formal optimization problem. In this section, we define the mathematical framework for generating sparse linear explanations and formulate our primary optimization objective. We then analyze the computational complexity of solving this problem exactly, demonstrating why empirical approximations are practically necessary.

F. Koriche et al.: Preprint submitted to Elsevier

Page 6 of 35

Probabilistic Linear Explanations

3.1. Notation Sets, Vectors, and Matrices. Plain letters represent functions and scalars, while boldface letters represent vectors

and matrices. We denote the all-ones and all-zeros vectors as 𝟏 and 𝟎, respectively. The identity matrix is denoted 𝑰. For a positive integer 𝑑, we use [𝑑] as shorthand for the index set {1, … , 𝑑}. The standard basis vectors of ℝ𝑑 are denoted 𝒆1 , … , 𝒆𝑑 , where 𝑒𝑖𝑗 = 1 if 𝑖 = 𝑗 and 0 otherwise. For a vector 𝒘 ∈ ℝ𝑑 and a subset 𝑆 ⊆ [𝑑], the projection of 𝒘 to the coordinates indexed by 𝑆 is the vector 𝒘𝑆 ∈ ℝ|𝑆| . The support set of 𝒘 ∈ ℝ𝑑 , denoted as support(𝒘), is the set of coordinates 𝑗 ∈ [𝑑] for which 𝑤𝑗 ≠ 0.

Operations and Norms. The transpose of a matrix 𝒁 is denoted 𝒁 ⊤ . The scalar product of two vectors 𝒗 and 𝒘 is

denoted as 𝒗 ⋅ 𝒘, and their coordinate-wise (or Hadamard) product as 𝒗 ⊙ 𝒘. For a scalar 𝑝 ∈ [1, ∞), the 𝐿𝑝 norm of 𝒘 is denoted as ‖𝒘‖𝑝 . Limit cases include the maximum absolute value ‖𝒘‖∞ = max𝑑𝑗=1 |𝑤𝑗 |. Following common usage, we extend this notation to the sparsity measure ‖𝒘‖0 = |support(𝒘)|—a measure that counts the number of non-zero coordinates.

Balls, Hyperplanes, and Projections. For a scalar 𝑝 ∈ [0, ∞] and a radius 𝑟 ≥ 0, the 𝐿𝑝 ball (centered at 𝟎) is defined as 𝑝 (𝑟) = {𝒘 ∈ ℝ𝑑 ∶ ‖𝒘‖𝑝 ≤ 𝑟}. For a normal vector 𝒖 ∈ ℝ𝑑 and a scalar 𝑟 ∈ ℝ, the corresponding hyperplane is defined as (𝒖, 𝑟) = {𝒘 ∈ ℝ𝑑 ∶ 𝒖 ⋅ 𝒘 = 𝑟}. Finally, the Euclidean projection of a vector 𝒘 ∈ ℝ𝑑 onto a set  ⊆ ℝ𝑑 is given by:2 Π (𝒘) = arg min ‖𝒘 − 𝒖‖2 𝒖∈

3.2. Problem Formulation In this study, we focus on explanation tasks where data instances are represented by a set of interpretable literals. Returning to our financial example, consider a bank customer who wants to understand why her loan interest rate was set at 8% (i.e., 𝑓 (𝒙) = 0.08). Interpretable literals such as [Income ≥ 65k], [DTI ≤ 30%], and [Credit History ≥ 10 yrs] can be used to construct the explanation. Each literal is assigned the polarity present in the customer’s profile, resulting in a clear and succinct if-then rule over weighted features, for example: −0.03 [DTI ≤ 30%] + 0.11 [Credit History < 10 yrs] → Rate = 0.08 where each coefficient is the contribution 𝑤𝑗 𝑥𝑗 of the corresponding literal, and the two contributions sum to the predicted rate. Formally, let [𝑑] be the set of interpretable literals. Treating them as binary features, the prediction models in this study are pseudo-Boolean functions of the form 𝑓 ∶ {−1, +1}𝑑 → [−1, +1]. Here, 𝑓 is a (binary) classifier if its codomain is {−1, +1}, and a regressor if its codomain is [−1, +1]. Each input to 𝑓 is a data instance 𝒙 ∈ {−1, +1}𝑑 , where 𝑥𝑗 ∈ {+1, −1} indicates whether the 𝑗-th literal is present positively or negatively in 𝒙. By convention, the first literal is the constant 𝑥1 = 1, so that a bias (or intercept) term 𝑤1 is always available; it is counted within the sparsity measure ‖𝒘‖0 like any other coefficient, which avoids unnecessary mathematical friction in what follows. A linear explanation for 𝑓 (𝒙) is a vector 𝒘 ∈ ℝ𝑑 that satisfies the anchoring condition 𝒘 ⋅ 𝒙 = 𝑓 (𝒙). As in the previous example, such an explanation can be interpreted as an if-then rule over weighted literals: the head is 𝑓 (𝒙), and the body consists of the literals 𝑗 ∈ support(𝒘), each taken with the polarity 𝑥𝑗 it has in 𝒙 and the weight 𝑤𝑗 . An explanation 𝒘 is considered 𝑘-sparse if ‖𝒘‖0 ≤ 𝑘. Throughout, the two constraints we impose on explanations are sparsity, which reflects the cognitive limits discussed in Section 1, and anchoring, which enforces consistency with the model at 𝒙. We formalize the set of all such valid explanations as our primary hypothesis space: 𝒙,𝑘 = (𝒙, 𝑓 (𝒙)) ∩ 0 (𝑘)

(2)

As illustrated in Figure 1, the geometric space 𝒙,𝑘 of 𝑘-sparse linear explanations for 𝑓 (𝒙) is the intersection of two distinct objects: the anchoring hyperplane (𝒙, 𝑓 (𝒙)) and the 𝐿0 ball 0 (𝑘). While the hyperplane is a ( ) continuous convex space, the 𝐿0 ball is a non-convex union of the 𝑑𝑘 coordinate subspaces of dimension 𝑘. Their intersection is therefore a combinatorial union of affine subspaces of dimension 𝑘−1, one for each candidate support— isolated points when 𝑘 = 1, as in Figure 1a, and lines when 𝑘 = 2, as in Figure 1b. None of them is empty, since 2 When the minimizer is not unique, ties are broken arbitrarily.

F. Koriche et al.: Preprint submitted to Elsevier

Page 7 of 35

Probabilistic Linear Explanations 1 𝑤3 1 𝑤2

(𝒙, 𝑓 (𝒙))

𝒙

1 2

1 2

𝒙

𝑤1 −1

− 12

0 − 12

1 2

1 0

(𝒙, 𝑓 (𝒙))

1 2

𝑤1 −1

(a) For 𝒙 = (1, 1) with 𝑓 (𝒙) = 21 and 𝑘 = 1, the hyperplane (𝒙, 𝑓 (𝒙)) appears in blue, and the 𝐿0 ball 0 (1) corresponds to the two black coordinate axes. Their intersection consists of the two red points ( 12 , 0) and (0, 12 ): one per candidate support, each of dimension 𝑘 − 1 = 0.

1 2

𝑤2 1

1

(b) For 𝒙 = (1, 1, 1) with 𝑓 (𝒙) = 12 and 𝑘 = 2, the hyperplane (𝒙, 𝑓 (𝒙)) appears as the large blue triangle, while the 𝐿0 ball 0 (2) is the union of the three coordinate planes 𝑤1 = 0, 𝑤2 = 0 and 𝑤3 = 0. Their intersection consists of the three red lines, one per candidate support, each of dimension 𝑘 − 1 = 1.

Figure 1: A geometric illustration of 𝑘-sparse linear explanations 𝒘.

𝑥𝑗 ≠ 0 for every 𝑗 ∈ [𝑑], so that every support carries linear explanations for 𝑓 (𝒙). This decomposition locates ( ) a key source of computational hardness: selecting a support is a combinatorial problem over 𝑑𝑘 alternatives. The quality of probabilistic explanations is assessed in relation to a probability distribution  over {−1, +1}𝑑 , such that (𝒙) > 0. For instance,  could represent the uniform distribution  across {−1, +1}𝑑 or, more restrictively, a localized neighborhood distribution surrounding the instance 𝒙 that is being explained. To quantify this quality, our framework adopts a normalized quadratic loss: 𝓁(𝑦, 𝑦′ ) = 14 (𝑦−𝑦′ )2 . This loss is applied consistently across both continuous regression and binary classification tasks. While its use in regression is standard, it is equally principled in the binary classification setting, where 𝑓 outputs values in {−1, +1}. The 1∕4 scaling ensures that the normalized quadratic loss exactly matches the zero-one loss 𝓁(𝑦, 𝑦′ ) = 𝟙[𝑦 ≠ 𝑦′ ] for 𝑦, 𝑦′ ∈ {−1, +1}. The theoretical justification and convergence guarantees for the quadratic loss as a surrogate in classification are well established in statistical learning theory (Bartlett, Jordan and McAuliffe, 2006; Suykens and Vandewalle, 1999) and supported by empirical results in machine learning (Rifkin and Klautau, 2004; Hui and Belkin, 2021). Definition 1 (Relevance Error). For a prediction model 𝑓 ∶ {−1, +1}𝑑 → [−1, +1], a reference instance 𝒙 ∈ {−1, +1}𝑑 , and a distribution , the relevance error of a vector 𝒘 ∈ ℝ𝑑 is given by:3 𝖱𝑓 ,𝒙, (𝒘) = 𝔼𝒛∼ [𝓁(𝑓 (𝒛), 𝑓 (𝒙)) ∣ 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙]

(3)

In other words, the relevance error of 𝒘 measures the expected discrepancy between 𝑓 (𝒛) and 𝑓 (𝒙) for random instances 𝒛 that are indistinguishable from 𝒙 under the linear projection defined by 𝒘. Notably, due to the exact equivalence established above, when 𝑓 is a binary classifier, this expectation perfectly recovers the original probabilistic objective formulated in (P1). With these concepts established, the decision version of the stochastic optimization problem presented in (P2) is defined as follows. Definition 2 (SLE Problem). An instance of the SPARSE L INEAR E XPLANATION (SLE) problem consists of a predictive model 𝑓 ∶ {−1, +1}𝑑 → [−1, +1], a data instance 𝒙 ∈ {−1, +1}𝑑 , a probability distribution  over {−1, +1}𝑑 , a sparsity level 𝑘 ≥ 1, and a relevance threshold 𝜏 ∈ [0, 1]. The question is whether there exists a linear explanation 𝒘 ∈ ℝ𝑑 for 𝑓 (𝒙) such that ‖𝒘‖0 ≤ 𝑘 and 𝖱𝑓 ,𝒙, (𝒘) ≤ 1 − 𝜏. 3 The conditioning event always contains 𝒛 = 𝒙, so that ℙ

𝒛∼ [𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙] ≥ (𝒙). Assuming (𝒙) > 0 is thus enough for 𝖱𝑓 ,𝒙, (𝒘) to be

well defined, simultaneously for every 𝒘 ∈ ℝ𝑑 .

F. Koriche et al.: Preprint submitted to Elsevier

Page 8 of 35

Probabilistic Linear Explanations

When 𝑓 (𝒙) = 0, the null vector 𝒘 = 𝟎 meets both constraints: it is 𝑘-sparse, and the anchoring condition holds trivially. However, this explanation is degenerate, as its conditioning event covers the entire instance space and its body is empty, providing no meaningful insight into the model. Consequently, we set this degenerate case aside in the remainder of the paper, while the complexity results below apply to the unrestricted problem. The above definition treats 𝑓 and  as abstract objects, as required for the remainder of the paper. Notably, our algorithms interact with 𝑓 solely through value queries—by evaluating 𝑓 (𝒛) at chosen instances 𝒛—and never through its internal structure. In contrast, a complexity statement requires an input with a measurable description length. To address this, we fix such representations in the next subsection: 𝑓 is specified as a neural network, and  is provided in closed form.

3.3. Problem Complexity (𝑑 ) Solving SLE exactly presents two intertwined challenges. First, we must select the correct support from among possible candidates. Second, for each candidate, we need to compute a probability over the exponentially 𝑘 many instances that agree with 𝒙 on that support. The first is a combinatorial search, the second a counting task. The complexity class that encompasses this combination is NPPP : problems solvable in polynomial time by a nondeterministic machine with access to a counting oracle. In this subsection, our goal is to show that SLE is hard for this class. Consequently, no algorithm can be expected to solve it exactly at scale, and the empirical relaxations developed in later sections are not merely convenient—they are essential. A natural approach is to reduce from the subset version of the problem, whose hardness was established by Wäldchen et al. (2021). Recall that a subset explanation 𝑆 selects instances agreeing with 𝒙 on 𝑆, while a linear explanation 𝒘 selects all instances 𝒛 such that 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙. If these two families of sets always coincided, the reduction would be immediate. However, in general, they do not, and discrepancies arise in both directions. This gap is the central obstacle addressed in the remainder of the subsection, so it is helpful to first illustrate it concretely before proceeding to the formal details. Example 3. Consider 𝑑 = 3, the instance 𝒙 = (+1, −1, +1), and the weight vector 𝒘 = (1, 1, 0), so that 𝒘 ⋅ 𝒙 = 0. The instance 𝒛 = (−1, +1, +1) also satisfies 𝒘 ⋅ 𝒛 = 0, even though it differs from 𝒙 on two coordinates. Thus, the set selected by 𝒘 is strictly larger than the one selected by the subset {1, 2}, as it also includes instances obtained by flipping both of the first two coordinates. In general, a linear explanation selects a union of such subsets—one passing through 𝒙 and others that do not. This difference is important for the reduction: converting a subset explanation into a linear one is straightforward, since the coefficients can be chosen so that the union collapses to a single set. The reverse direction is more subtle, since a linear explanation may owe its quality to a set that does not pass through 𝒙, and such a set does not necessarily certify the subset problem. The proof proceeds in three movements. First, any linear explanation can be replaced by one of the subcubes it selects, without increasing the relevance error (Lemma 1); this reduces the linear problem to a subset-like problem, though at the cost of losing the guarantee that the subcube passes through 𝒙. Second, we make this loss harmless: by replacing selected variables with the parity of 𝑘 + 1 fresh copies, a budget of 𝑘 can never extract any information about them (Lemma 2 and Corollary 1), so only the coordinates of interest can be exploited. Third, for the resulting family of instances, a subcube meeting the relevance threshold exists if and only if an anchored one does—a property we call alignment (Definition 5)—and both are equivalent to the existence of a witness for the source problem (Lemma 3 and Corollary 2). The reduction then follows (Theorem 1). Before proceeding, we clarify the representation of the model and the source problem.

Model Representation. For a meaningful complexity statement, the input must have a finite, measurable description

length. Accordingly, we assume that 𝑓 is specified as a feedforward ReLU neural network, with all weights and biases in [−1, +1] and activation functions restricted to the identity 𝑢 ↦ 𝑢 and the rectifier 𝑢 ↦ max{0, 𝑢}. The description length of 𝑓 is measured by its number of gates. These networks can emulate Boolean circuits of comparable size and depth (Mukherjee and Basu, 2017; Parberry, 1996), so we may equivalently describe a binary classifier as a circuit built from the connectives ∧ (AND), ∨ (OR), ¬ (NOT), and ⊕ (XOR), counting its gates. While circuits are traditionally defined over {0, 1}𝑑 , we retain the symmetric cube {−1, +1}𝑑 throughout this paper, identifying the bit 𝑏 with 2𝑏 − 1. This bijection preserves both the set of instances and the uniform distribution, ensuring every result from (Wäldchen et al., 2021) applies directly in our encoding. Finally, we assume that (𝒛) can be evaluated in polynomial time. F. Koriche et al.: Preprint submitted to Elsevier

Page 9 of 35

Probabilistic Linear Explanations

Source Problem. All the hardness in this subsection ultimately comes from E-M AJ-SAT (Littman, Goldsmith and

Mundhenk, 1998), the canonical complete problem for NPPP . This problem is the probabilistic analog of SAT: instead of asking whether some assignment of the free variables satisfies the formula, it asks whether most of them do.

Definition 3 (E-M AJ-SAT Problem). An instance of E-M AJ-SAT is a Boolean formula 𝜑 in conjunctive normal form over 𝑛 variables, together with an integer 𝜅 ≤ 𝑛. The question is whether there exists an assignment 𝒖∗ of the first 𝜅 variables, called a witness, such that a majority of the assignments of the remaining 𝑛 − 𝜅 variables satisfy 𝜑. The existential choice of the witness underlies the problem’s NP-hardness, while the majority test accounts for its PP-hardness. These are precisely the two ingredients found in SLE: selecting a support and then performing counting within it.

Subcubes. The sets discussed above have a standard name. Given a subset 𝑆 ⊆ [𝑑] and an assignment 𝒂 ∈ {−1, +1}𝑆 ,

the subcube (𝑆, 𝒂) = {𝒛 ∶ 𝒛𝑆 = 𝒂} is the set of instances agreeing with 𝒂 on 𝑆, and |𝑆| is its codimension. A subcube is anchored when 𝒂 = 𝒙𝑆 , that is, when it passes through the instance being explained. It is convenient to assess the quality of a subcube in exactly the same way as for linear explanations, by extending the relevance error of (3) to any event 𝐸 with positive probability: 𝖱𝑓 ,𝒙, (𝐸) = 𝔼𝒛∼ [𝓁(𝑓 (𝒛), 𝑓 (𝒙)) ∣ 𝒛 ∈ 𝐸] With this convention, 𝖱𝑓 ,𝒙, (𝒘) represents the relevance error of the set {𝒛 ∶ 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙} selected by 𝒘; a 𝑘-sparse subset explanation corresponds to an anchored subcube of codimension at most 𝑘, and a single threshold 1 − 𝜏 applies in both cases. The first lemma formalizes the observation from the example: a linear explanation selects a union of subcubes of equal size, so its relevance error is the average of their errors. Since an average cannot be less than its minimum, at least one subcube must achieve a relevance error no greater than that of the linear explanation itself. Lemma 1 (Extraction). Let 𝒘 ∈ ℝ𝑑 with 𝑆 = support(𝒘). Under the uniform distribution  , the set selected by 𝒘 is the disjoint union of the subcubes (𝑆, 𝒂) over the assignments 𝒂 satisfying 𝒘𝑆 ⋅ 𝒂 = 𝒘𝑆 ⋅ 𝒙𝑆 . Moreover, at least one of them satisfies 𝖱𝑓 ,𝒙, ((𝑆, 𝒂)) ≤ 𝖱𝑓 ,𝒙, (𝒘) PROOF. The value 𝒘 ⋅ 𝒛 depends on 𝒛 only through 𝒛𝑆 , so the set selected by 𝒘 is indeed the union of the subcubes of codimension |𝑆| described above, and this union is disjoint since distinct assignments of 𝑆 define disjoint subcubes. Each of these subcubes contains 2𝑑−|𝑆| instances, hence carries the same conditional mass under  , so that 𝖱𝑓 ,𝒙, (𝒘) is the plain average of their relevance errors. A minimum is never larger than an average. The subcube produced by Lemma 1 is not necessarily anchored, which is the gap we must address. We bridge this gap by ensuring that unanchored subcubes are uninformative, using a simple device: if a variable is replaced by the parity of 𝑘 + 1 fresh copies, then fixing at most 𝑘 of these leaves at least one unfixed, so the parity remains an unbiased coin. As a result, no explanation with budget 𝑘 can extract any information about that variable, regardless of the values assigned. Definition 4 (Shielding). Let 𝑔 be a function of 𝑛 variables and let 𝐴 be a subset of them. The 𝑘-shielding of 𝑔 over ⨁ 𝐴 is obtained by replacing every variable 𝑦𝑖 with 𝑖 ∈ 𝐴 by the parity 𝑘+1 𝑗=1 𝑦𝑖,𝑗 of 𝑘 + 1 fresh variables. It adds 𝑘|𝐴| variables and 𝑂(𝑘|𝐴|) gates. Lemma 2 (Shielding). Fix any assignment of at most 𝑘 of the fresh variables, so that every block of 𝑘 + 1 copies ⨁ retains at least one free variable. Given this assignment, the vector of parities ( 𝑗 𝑦𝑖,𝑗 )𝑖∈𝐴 is uniformly distributed and independent of all the other variables. PROOF. Each 𝑖 ∈ 𝐴 is associated with a block of 𝑘+1 fresh variables. Since at most 𝑘 variables are fixed in total, every block retains at least one unfixed variable. The parity of each block, being the parity of a nonempty set of independent unbiased bits, is itself unbiased and remains independent across blocks as well as from the unshielded variables. F. Koriche et al.: Preprint submitted to Elsevier

Page 10 of 35

Probabilistic Linear Explanations

We use this device as follows: a candidate explanation does not benefit by allocating any part of its budget to fresh variables, so we may always assume it spends none on them. Corollary 1 (Shielded Coordinates are Useless). Let 𝑓̃ be the 𝑘-shielding of a function 𝑓 , and let 𝐵 denote the set of fresh variables it introduces. Then, for every subcube (𝑆, 𝒂) of codimension at most 𝑘, ( ) 𝖱𝑓̃,𝒙, ((𝑆, 𝒂)) = 𝖱𝑓̃,𝒙, (𝑆 ⧵ 𝐵, 𝒂𝑆⧵𝐵 ) In particular, for every subcube of codimension at most 𝑘 there is another one, of no larger codimension, that fixes no fresh variable and has the same relevance error. PROOF. By Lemma 2, conditioning on the coordinates of 𝑆 ∩ 𝐵 leaves the vector of parities uniform and independent of every other variable. The conditional law of the arguments of the shielded function, and therefore that of 𝑓̃(𝒛), is thus the same as under the conditioning by 𝑆 ⧵ 𝐵 alone, and the two conditional expectations coincide. We can now name the property that our family of instances must satisfy. It is a statement at the threshold rather than an equality between optima: an unanchored subcube may well be strictly better than every anchored one, and this is harmless as long as it does not cross the threshold separating positive from negative instances. Definition 5 (Alignment). An instance (𝑓 , 𝒙, 𝑘) is aligned at level 𝜏 if, whenever some subcube of codimension at most 𝑘 has relevance error at most 1 − 𝜏, some anchored subcube of codimension at most 𝑘 has relevance error at most 1 − 𝜏 as well. On an aligned instance, the anchored and the free versions of the problem have the same answer, which is what our two directions require: the forward direction builds a linear explanation from an anchored subcube, while the backward direction receives an arbitrary subcube from Lemma 1.

The Construction. Let (𝜑, 𝜅) be an instance of E-M AJ-SAT and let 𝜏 ∈ (0, 1). We build a classifier 𝑓 , a reference

instance 𝒙 and a budget 𝑘 = 𝜅 in two steps borrowed from Wäldchen et al. (2021), followed by one application of Definition 4. The first step is the duplication construction of Wäldchen et al. (2021, Lemma 3.5). It produces a circuit Φ over two blocks 𝒖 and 𝒗 of size 𝜅 each, a block 𝒓 holding the remaining variables of 𝜑, and one extra variable. Its reference instance carries opposite values on the two duplicated blocks, so that for each index 𝑖, exactly one of the two 𝑖-th coordinates encodes a given value (True or False) of the 𝑖-th variable of 𝜑. Let EQ denote the event that the two duplicated blocks agree; conditionally on EQ, fixing exactly one of these two coordinates to its value in 𝒙 assigns that value to the 𝑖-th variable of 𝜑. For a subcube  fixing coordinates of 𝒖 and 𝒗 only, let 𝗊() denote the probability that Φ evaluates to +1 given . Their analysis (Wäldchen et al., 2021, Lemma 3.3) gives ( [ ) [ ] ] 𝗊() = 21 + ℙ 𝜑 ∣ , EQ − 12 ℙ EQ ∣  (4) No assignment of the variables of 𝜑 appears here: a subcube may leave coordinates free, in which case the first conditional probability averages over them. Note also that ℙ[EQ ∣ ] > 0 whenever  leaves at least one coordinate of each duplicated pair free, or fixes both to equal values, so that the conditional probability ℙ[𝜑 ∣ , EQ] is well defined in these cases; when  fixes both coordinates of some pair to opposite values, ℙ[EQ ∣ ] = 0 and (4) reads 𝗊() = 21 , with the second term interpreted as zero. The second step is the threshold adjustment from Wäldchen et al. (2021, Lemmas 3.8–3.11). This transforms Φ into the final circuit 𝑓 by composing it with auxiliary blocks, ensuring that 𝑓 (𝒙) = +1 and that, for every subcube  fixing coordinates of 𝒖 and 𝒗 only, 𝖱𝑓 ,𝒙, () ≤ 1 − 𝜏

⟺

𝗊() > 12

(5)

The null vector is not a feasible explanation here, since 𝑓 (𝒙) = +1 ≠ 0. Finally, we 𝑘-shield every variable except those belonging to the two duplicated blocks: specifically, the block 𝒓, the extra variable from the first step, and all auxiliary blocks from the second step are shielded. By Lemma 2, this modification does not affect any of the probabilities established above, and the construction remains polynomial since F. Koriche et al.: Preprint submitted to Elsevier

Page 11 of 35

Probabilistic Linear Explanations

shielding multiplies the number of variables in a block by 𝑘 + 1 and adds a proportional number of gates. The resulting classifier, still denoted 𝑓 , is a 𝑘-shielding, and therefore falls under Corollary 1. Shielding also renders the intermediate step of (Wäldchen et al., 2021, Lemma 3.7) unnecessary: its sole purpose was to prevent a budget from being allocated to 𝒓, which is now directly guaranteed by Corollary 1. One further notion completes the picture. The subcube induced by an assignment 𝒖∗ to the first 𝜅 variables of 𝜑, denoted  ∗ (𝒖∗ ), is defined as the subcube that fixes, for each 𝑖 ≤ 𝜅, the 𝑖-th coordinate of 𝒖 or 𝒗 that encodes 𝑢∗𝑖 to its value in 𝒙. By construction, this subcube is anchored, and its codimension is 𝜅 = 𝑘. The next lemma forms the technical heart of the reduction: it establishes the connection between subcubes in the construction and witnesses of the source instance, working in both directions. Lemma 3 (Witnesses and Subcubes). For every instance produced by the construction above: 1. if some subcube of codimension at most 𝑘 has relevance error at most 1 − 𝜏, then 𝜑 admits a witness; 2. conversely, for every witness 𝒖∗ , the induced subcube  ∗ (𝒖∗ ) has relevance error at most 1 − 𝜏. PROOF. Clause 1. Let  be such a subcube. By Corollary 1, we may assume that it fixes only coordinates of 𝒖 and 𝒗, since replacing it by its unshielded part changes neither its relevance error nor the bound on its codimension. Then ) ( (5) gives 𝗊() > 21 , so that by (4) the product ℙ[𝜑 ∣ , EQ] − 12 ℙ[EQ ∣ ] is strictly positive. Neither factor can therefore vanish, and both must share the same sign. Only the sign of the second factor matters here, not its magnitude: since ℙ[EQ ∣ ] is a probability, it is strictly positive, and hence ℙ[𝜑 ∣ , EQ] > 12 as well. Conditionally on  and EQ, the two duplicated blocks agree, the coordinates pinned by  carry the values they encode, the remaining ones are uniform, and 𝒓 is independent of all of them. The value ℙ[𝜑 ∣ , EQ] is thus an average over assignments of the first 𝜅 variables compatible with , of the probability that 𝜑 is satisfied by a uniform assignment to the remaining variables. At least one such assignment must yield a value exceeding 12 , and is thus a witness. Clause 2. Let 𝒖∗ be a witness and  ∗ =  ∗ (𝒖∗ ). This fixes one coordinate in each of the 𝜅 duplicated pairs and leaves the other free, so ℙ[EQ ∣  ∗ ] = 2−𝜅 > 0, and conditioning further on EQ pins the first 𝜅 variables of 𝜑 to 𝒖∗ . Then ℙ[𝜑 ∣  ∗ , EQ] > 12 because 𝒖∗ is a witness, and (4) gives 𝗊( ∗ ) > 21 , so (5) concludes. Both properties we need now follow formally, the first one by chaining the two clauses. Corollary 2 (Alignment). Let 𝜏 ∈ (0, 1). The construction above maps, in polynomial time, any instance of E-M AJSAT to a triple (𝑓 , 𝒙, 𝑘) such that 1. (𝑓 , 𝒙, 𝑘) is aligned at level 𝜏; and 2. the E-M AJ-SAT instance is positive if and only if some anchored subcube of codimension at most 𝑘 has relevance error at most 1 − 𝜏. PROOF. For clause 1, suppose that some subcube of codimension at most 𝑘 has relevance error at most 1 − 𝜏. Clause 1 of Lemma 3 yields a witness 𝒖∗ , and clause 2 turns it into the subcube  ∗ (𝒖∗ ), which is anchored, of codimension at most 𝑘, and of relevance error at most 1 − 𝜏. This is precisely Definition 5. For clause 2, a positive instance has a witness, hence an anchored subcube below the threshold by clause 2 of Lemma 3. Conversely, an anchored subcube below the threshold is in particular a subcube, so clause 1 of Lemma 3 produces a witness. In particular, if the E-M AJ-SAT instance is negative, then no anchored subcube of codimension at most 𝑘 meets the threshold, and thus, by alignment, no subcube of any kind does. This contrapositive form is the one to keep in mind when reading the backward direction below. Theorem 1 (Hardness). For every fixed relevance threshold 𝜏 ∈ (0, 1), the SPARSE L INEAR E XPLANATION problem for ReLU neural networks is NPPP -hard. PROOF. Let 𝑓 , 𝒙, and 𝑘 be produced by Corollary 2 from an instance of E-M AJ-SAT, and consider the SLE instance (𝑓 , 𝒙,  , 𝑘, 𝜏). This instance is well defined, as  (𝒙) > 0, and its construction is polynomial in size. Suppose first that the E-M AJ-SAT instance is positive, and let (𝑆, 𝒙𝑆 ) be the anchored subcube given by clause 2 of Corollary 2. If 𝑆 is empty, set 𝒘 = 𝒆1 ; since 𝑥1 = 1 is the constant literal, 𝒘 ⋅ 𝒛 = 1 = 𝑓 (𝒙) for every instance 𝒛, F. Koriche et al.: Preprint submitted to Elsevier

Page 12 of 35

Probabilistic Linear Explanations

so 𝒘 selects the entire cube, which is precisely the subcube of empty support. Otherwise, set 𝑤𝑗 = 𝑥𝑗 ∕|𝑆| for 𝑗 ∈ 𝑆 and 𝑤𝑗 = 0 elsewhere; the nonzero coefficients are then equal in magnitude, so the condition 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙 reduces to ∑ 𝑗∈𝑆 𝑥𝑗 𝑧𝑗 = |𝑆|, and since each term is at most one, this forces 𝒛𝑆 = 𝒙𝑆 : the set selected by 𝒘 is exactly (𝑆, 𝒙𝑆 ). In both cases 𝒘 ⋅ 𝒙 = 1 = 𝑓 (𝒙) and ‖𝒘‖0 ≤ 𝑘, so 𝒘 is a feasible explanation whose relevance error matches that of the subcube, hence is at most 1 − 𝜏, and the SLE instance is positive. Conversely, suppose the SLE instance is positive, witnessed by a feasible 𝒘, and let 𝑆 = support(𝒘), so that |𝑆| ≤ 𝑘. Lemma 1 guarantees the existence of a subcube (𝑆, 𝒂) whose relevance error is at most that of 𝒘, and thus at most 1 − 𝜏. If 1 ∈ 𝑆, fixing the constant coordinate is vacuous, so this subcube coincides with (𝑆 ⧵ {1}, 𝒂𝑆⧵{1} ), of no larger codimension. Because the instance is aligned at level 𝜏, some anchored subcube of codimension at most 𝑘 also achieves this threshold, and clause 2 of Corollary 2 then ensures the existence of a witness. The two instances are therefore equivalent, and since E-M AJ-SAT is NPPP -complete (Littman et al., 1998), the result follows.

4. Dealing with PP-Hardness Among the two sources of intractability highlighted in Section 3, only one is out of the ordinary. The combinatorial search over candidate supports 𝑆 ⊆ [𝑑] of size at most 𝑘 is the familiar cost of enforcing sparsity, and an extensive literature exists for addressing it—ranging from exact encodings to greedy and thresholding algorithms. The second source is more problematic: even for a fixed candidate 𝒘, determining whether its relevance error is at most 1 − 𝜏 is PP-hard. This complexity places the problem beyond the scope of standard optimization techniques, which require that the objective be efficiently computable. This section addresses the second obstacle by introducing the fidelity error, an unconditional surrogate loss. Section 4.1 formally relates relevance and fidelity errors within the parameterized family of neighborhood distributions, while Section 4.2 demonstrates that the fidelity error can be efficiently approximated via sampling. This allows for a tractable Mixed Integer Programming (MIP) formulation with formal sample complexity guarantees.

4.1. From Relevance to Fidelity To circumvent the PP-hard evaluation of the conditional relevance error, we focus on the fidelity error. Commonly used in model-agnostic explainability (Lakkaraju et al., 2019; Li et al., 2021; Ribeiro et al., 2016), the fidelity error assesses the expected loss of the linear explanation over the entire distribution, omitting the conditioning event 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙. The anchoring constraint 𝒘 ⋅ 𝒙 = 𝑓 (𝒙) remains an essential part of the definition and will continue to play a central role below. By evaluating the loss unconditionally, this surrogate metric becomes much more tractable for both theoretical analysis and empirical approximation. Definition 6 (Fidelity Error). For a prediction model 𝑓 , a reference instance 𝒙, and a distribution , the fidelity error of a vector 𝒘 under the normalized quadratic loss 𝓁 is given by: [ ( )] 𝖥𝑓 ,𝒙, (𝒘) = 𝔼𝒛∼ 𝓁 𝒘 ⋅ 𝒛, 𝑓 (𝒛) (6) To bridge the theoretical gap between the conditional relevance error and the unconditional fidelity error, we rely on the parameterized family of localized distributions introduced in Section 1. Recall that for a given concentration parameter 𝜎 ≥ 0 and a central reference instance 𝒙, the probability of sampling an instance 𝒛 decays exponentially with its Hamming distance 12 ‖𝒙 − 𝒛‖1 . Formally, we define the neighborhood distribution 𝒙,𝜎 as: 𝒙,𝜎 (𝒛) =

1 − 𝜎2 ‖𝒙−𝒛‖1 𝑒 𝑍𝜎

where

𝑍𝜎 =

𝑑 ( ) ∑ 𝑑 −𝜎𝑗 𝑒 = (1 + 𝑒−𝜎 )𝑑 𝑗 𝑗=0

(7)

As previously noted, the parameter 𝜎 directly controls the locality of the explanation. When 𝜎 = 0, the decay vanishes and 𝒙,0 recovers the uniform distribution  over the entire instance space. Conversely, as 𝜎 → ∞, the distribution concentrates entirely on the center 𝒙. Importantly, for neighborhood distributions, the fidelity error provides an upper bound for the relevance error of any explanation in 𝒙,𝑘 = (𝒙, 𝑓 (𝒙)) ∩ 0 (𝑘), modulated by the sparsity level 𝑘 and the concentration 𝜎.

F. Koriche et al.: Preprint submitted to Elsevier

Page 13 of 35

Probabilistic Linear Explanations

Lemma 4 (Approximating Relevance via Fidelity). Let 𝑓 ∶ {−1, +1}𝑑 → [−1, +1] be a prediction model, let 𝒙 ∈ {−1, +1}𝑑 be a data instance, let 𝑘 ≥ 1 be a sparsity parameter, and let 𝜎 ≥ 0 be a concentration parameter. Then, for any explanation 𝒘 ∈ 𝒙,𝑘 , its relevance error satisfies 𝖱𝑓 ,𝒙,𝒙,𝜎 (𝒘) ≤ (1 + 𝑒−𝜎 )𝑘 𝖥𝑓 ,𝒙,𝒙,𝜎 (𝒘) PROOF. Let 𝑆 = support(𝒘), so that |𝑆| ≤ 𝑘 since 𝒘 ∈ 0 (𝑘), and let 𝐴 = {𝒛 ∈ {−1, +1}𝑑 ∶ 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙} be the conditioning event. Since 𝒘 is anchored (i.e., 𝒘 ∈ (𝒙, 𝑓 (𝒙))), 𝑓 (𝒙) = 𝒘 ⋅ 𝒙, hence 𝑓 (𝒙) = 𝒘 ⋅ 𝒛 for every 𝒛 ∈ 𝐴. Therefore [ ( ) ] 𝖱𝑓 ,𝒙,𝒙,𝜎 (𝒘) = 𝔼 𝓁 𝒘 ⋅ 𝒛, 𝑓 (𝒛) ∣ 𝐴 =

( ) 𝖥𝑓 ,𝒙,𝒙,𝜎 (𝒘) 1 ∑ 𝒙,𝜎 (𝒛) 𝓁 𝒘 ⋅ 𝒛, 𝑓 (𝒛) ≤ , ℙ[𝐴] 𝒛∈𝐴 ℙ[𝐴]

where the first equality uses the symmetry of 𝓁, and the inequality follows from the non-negativity of the summands. It remains to lower-bound ℙ[𝐴]. As 𝒘 vanishes outside 𝑆, the event 𝒛𝑆 = 𝒙𝑆 implies 𝐴. Moreover, ‖𝒙−𝒛‖1 decomposes coordinatewise, so 𝒙,𝜎 is a product distribution with marginals ℙ[𝑧𝑗 = 𝑥𝑗 ] = 1∕(1 + 𝑒−𝜎 ). Hence ℙ[𝐴] ≥ ℙ[𝒛𝑆 = 𝒙𝑆 ] = (1 + 𝑒−𝜎 )−|𝑆| ≥ (1 + 𝑒−𝜎 )−𝑘 . It should be stressed that this relationship is one-sided. Decomposing the fidelity error over the conditioning event ̄ 𝔼[𝓁 ∣ 𝐴], ̄ and the second term is not controlled by the relevance 𝐴 and its complement gives 𝖥 = ℙ[𝐴] 𝖱 + ℙ[𝐴] error: an explanation can be highly relevant yet incur a large fidelity error, if it behaves poorly outside the anchoring subcube. Minimizing fidelity is therefore a conservative surrogate—it certifies relevance, but may overlook relevant explanations. Lemma 4 quantifies the price of this conservatism, which the discussion below shows to be a small constant in the regimes of interest. Lemma 4 has important practical consequences for our framework. The bound between the fidelity and relevance errors becomes tighter as the concentration parameter 𝜎 increases, since the term 𝑒−𝜎 drops off rapidly. Although a large sparsity level 𝑘 could, in principle, amplify the exponential factor (1 + 𝑒−𝜎 )𝑘 , such high values of 𝑘 are at odds with the very goal of interpretability. As discussed in Section 1, human cognitive limits necessitate highly sparse explanations, generally restricting 𝑘 to four or five features. In the regimes of interest, the factor therefore remains a small constant: for 𝑘 = 5, it is about 4.8 at 𝜎 = 1, 1.9 at 𝜎 = 2, and 1.3 at 𝜎 = 3. This alignment between cognitive constraints and mathematical properties ensures that, for real-world XAI applications, the fidelity error provides a reliable proxy for the computationally intractable relevance error.

4.2. Empirical Approximation and MIP Formulation Interestingly, the fidelity error in (6) involves an unconditional expectation, which is highly amenable to sampling. To approximate this expectation, let {(𝒛𝑖 , 𝑓 (𝒛𝑖 ))}𝑚 be a sample set where each 𝒛𝑖 is drawn independently at random 𝑖=1 according to 𝒙,𝜎 , and its target value 𝑓 (𝒛𝑖 ) is obtained through query access to 𝑓 . The corresponding empirical fidelity error is given by: 𝑚 ) 1∑ ( 1 ̂ 𝖥𝑓 ,𝒙,𝑚 (𝒘) = 𝓁 𝒘 ⋅ 𝒛𝑖 , 𝑓 (𝒛𝑖 ) = ‖𝒁𝒘 − 𝒚‖22 𝑚 𝑖=1 4𝑚

(8)

where 𝒁 ∈ {−1, +1}𝑚×𝑑 collects the sampled instances as rows and 𝒚 = (𝑓 (𝒛1 ), … , 𝑓 (𝒛𝑚 )) collects their target values. Based on this objective function, (P3) takes the form of a variant of the well-studied problem known as sparse regression, also referred to as best subset selection, which dates back at least to (Beale et al., 1967; Hocking and Leslie, 1967). While this problem is non-convex and NP-hard (Natarajan, 1995), Bertsimas, King and Mazumder (2016) showed that Mixed Integer Programming (MIP) formulations make it practically solvable at scale, an approach that has since benefited from steady progress in branch-and-cut solvers. The following formulation is a variation of their parameter-free approach utilizing Specially Ordered Sets (SOS) (Bertsimas and Weismantel, 2005):

F. Koriche et al.: Preprint submitted to Elsevier

Page 14 of 35

Probabilistic Linear Explanations

minimize

𝑚 ) 1∑ ( 𝓁 𝒘 ⋅ 𝒛𝑖 , 𝑓 (𝒛𝑖 ) 𝑚 𝑖=1

subject to

𝒘 ⋅ 𝒙 = 𝑓 (𝒙) (MIP)

𝟏⋅𝒖≤𝑘 ‖(𝑤𝑗 , 1 − 𝑢𝑗 )‖0 ≤ 1, 𝑢𝑗 ∈ {0, 1},

for all 𝑗 ∈ [𝑑]

for all 𝑗 ∈ [𝑑]

𝑤𝑗 ∈ [−1, +1],

for all 𝑗 ∈ [𝑑]

The constraint ‖(𝑤𝑗 , 1 − 𝑢𝑗 )‖0 ≤ 1 formally encodes a Specially Ordered Set of type 1: if 𝑢𝑗 = 0 then 𝑤𝑗 is forced to zero, and if 𝑢𝑗 = 1 then 𝑤𝑗 is free. Since 𝑤𝑗 ∈ [−1, +1], this is equivalent to the Big-M constraints −𝑢𝑗 ≤ 𝑤𝑗 ≤ 𝑢𝑗 , which modern MIP solvers handle with high efficiency. The formulation is also guaranteed to be feasible: as 𝑘 ≥ 1 and |𝑓 (𝒙)| ≤ 1, the vector 𝒘 = 𝑓 (𝒙) 𝒆1 satisfies every constraint, and it will serve as the admissible starting point for the iterative method developed in Section 5. Compared to (P3), this formulation introduces the bounding constraint 𝑤𝑗 ∈ [−1, +1] for all 𝑗 ∈ [𝑑]. From the standpoint of interpretability, the restriction is mild: since the output of 𝑓 ranges over [−1, +1], a coefficient of magnitude greater than one would attribute to a single feature an effect exceeding the entire range of the model, which runs against the very purpose of feature attribution. Technically, however, the constraint is essential. Combined with the sparsity constraint, it guarantees that ‖𝒘‖1 ≤ ‖𝒘‖0 ≤ 𝑘, and it is this 𝐿1 bound—the radius 𝐵 of Theorem 2 below—that governs the sample complexity of the formulation; without it, no finite sample size would suffice. The following result shows that if the solver is supplied with a number of samples 𝑚 that is quartic in 𝑘 and logarithmic in 𝑑, then with high probability, the relevance error of every feasible explanation—and hence, in particular, that of the solution returned by the solver—is upper-bounded by a function of its empirical fidelity. We state it in terms of a generic 𝐿1 radius 𝐵, the box-constrained case of (MIP) corresponding to 𝐵 = 𝑘; Section 5 applies the same result with a larger radius. Theorem 2 (Approximating Relevance via Empirical Fidelity). Let 𝑓 ∶ {−1, +1}𝑑 → [−1, +1] be a prediction model, 𝒙 ∈ {−1, +1}𝑑 be a data instance, 𝜎 ≥ 0 be a concentration parameter, 𝑘 ≥ 1 be a sparsity parameter, and 𝐵 ≥ 1 be a magnitude parameter. Then, for any 𝛿 ∈ (0, 1] and any 𝜀 ∈ (0, 1], if 𝑚≥

( )) (𝐵 + 1)4 ( 16 ln(2𝑑) + ln 𝛿2 4 𝜀2

then with probability at least 1 − 𝛿 over the choice of an i.i.d. sample set of size 𝑚, every explanation 𝒘 ∈ 𝒙,𝑘 with ‖𝒘‖1 ≤ 𝐵 satisfies ( ) ̂𝑓 ,𝒙,𝑚 (𝒘) + 𝜀 𝖱𝑓 ,𝒙,𝒙,𝜎 (𝒘) ≤ (1 + 𝑒−𝜎 )𝑘 𝖥 PROOF. The result relies on standard uniform convergence bounds for linear predictors over 𝐿1 spaces (Shalev-Shwartz and Ben-David, 2014; Kakade, Sridharan and Tewari, 2008). The underlying setting is as follows: for a scalar 𝐵, consider a bounded input space  ⊆ {𝒛 ∈ ℝ𝑑 ∶ ‖𝒛‖∞ ≤ 1} and an output space  ⊆ ℝ, alongside the hypothesis class 1 (𝐵) = {𝒘 ∈ ℝ𝑑 ∶ ‖𝒘‖1 ≤ 𝐵} of 𝐿1 -bounded linear functions. In addition, given two scalars 𝜌 and 𝑐, let 𝓁 ∶  ×  → ℝ be a loss function such that, for all 𝑦 ∈ , the mapping 𝑎 ↦ 𝓁(𝑎, 𝑦) is 𝜌-Lipschitz, and the magnitude max𝑎∈[−𝐵,+𝐵] |𝓁(𝑎, 𝑦)| is upper bounded by 𝑐. Then, as established by Theorem 26.15 in Shalev-Shwartz and Ben-David (2014), for any distribution  over  × , with probability at least 1 − 𝛿 over the choice of an i.i.d. sample set of size 𝑚, every 𝒘 ∈ 1 (𝐵) satisfies: √ √ 𝑚 2 ln(2∕𝛿) 2 ln(2𝑑) 1∑ 𝓁(𝒘 ⋅ 𝒛𝑖 , 𝑦𝑖 ) + 2𝜌𝐵 +𝑐 𝔼(𝒛,𝑦)∼ [𝓁(𝒘 ⋅ 𝒛, 𝑦)] ≤ 𝑚 𝑖=1 𝑚 𝑚 In the setting of our framework,  = {−1, +1}𝑑 , which meets the above requirement since ‖𝒛‖∞ = 1 for all 𝒛 ∈ . Furthermore, the explanations of interest belong to our primary hypothesis space 𝒙,𝑘 and have an 𝐿1 norm of at most F. Koriche et al.: Preprint submitted to Elsevier

Page 15 of 35

Probabilistic Linear Explanations

𝐵; they therefore lie in 1 (𝐵), so the above bound applies to all of them with this radius, the anchoring and sparsity constraints defining 𝒙,𝑘 only shrinking the class and hence only helping. Note that any 𝒘 feasible for (MIP) satisfies ‖𝒘‖1 ≤ ‖𝒘‖0 ≤ 𝑘, since each non-zero coordinate is bounded by 1 in absolute value; the box-constrained case is thus recovered by setting 𝐵 = 𝑘. The radius is kept explicit because Section 5 applies this result to iterates that are not confined to the unit box. It remains to instantiate the two constants of the loss 𝓁(𝑎, 𝑦) = 41 (𝑎 − 𝑦)2 . For any 𝒛 ∈ {−1, +1}𝑑 and any such 𝒘, the prediction 𝒘 ⋅ 𝒛 is bounded in [−𝐵, 𝐵] and the true target 𝑓 (𝒛) is bounded in [−1, +1], so their absolute difference is at most 𝐵 + 1. The magnitude of the loss is therefore bounded by 𝑐 = (𝐵 + 1)2 ∕4. In addition, the derivative of 𝓁 with respect to the prediction is 12 (𝒘 ⋅ 𝒛 − 𝑓 (𝒛)), whose magnitude is at most (𝐵 + 1)∕2, so that 𝜌 = (𝐵 + 1)∕2. Based on these observations, and upper-bounding 2𝜌𝐵 = 𝐵(𝐵 + 1) by (𝐵 + 1)2 , it follows that with probability at least 1 − 𝛿, every such explanation 𝒘 satisfies: ̂𝑓 ,𝒙,𝑚 (𝒘) + 𝖥𝑓 ,𝒙,𝒙,𝜎 (𝒘) ≤ 𝖥

) √ (𝐵 + 1)2 (√ 2 ln(2𝑑) + 14 2 ln(2∕𝛿) √ 𝑚

Requiring the deviation term to be at most 𝜀, solving for 𝑚, and applying the algebraic inequality (𝑎 + 𝑏)2 ≤ 2𝑎2 + 2𝑏2 ̂𝑓 ,𝒙,𝑚 (𝒘) + 𝜀 for every such yields the lower bound on 𝑚 stated in Theorem 2. Under this condition, 𝖥𝑓 ,𝒙,𝒙,𝜎 (𝒘) ≤ 𝖥 −𝜎 𝑘 explanation, and multiplying both sides by (1 + 𝑒 ) before applying Lemma 4 concludes the proof. Two important remarks follow. First, for a fixed radius 𝐵, the required sample size is independent of the concentration parameter 𝜎; a single bound applies across the entire family of neighborhood distributions, from the uniform case to sharply localized scenarios. This is a genuine advantage of (MIP): its box-constrained feasible set ensures 𝐵 = 𝑘 for any 𝜎; Section 5 shows that, in contrast, the polynomial-time alternative incurs a cost for locality, with its radius increasing alongside 𝜎. Second, since these guarantees are worst-case results reliant on uniform convergence, the constants involved are quite conservative. For example, setting 𝑘 = 5, 𝑑 = 100, 𝜀 = 0.1, and 𝛿 = 0.05, yields a theoretical requirement of roughly 2.9 × 106 samples. In practice, much smaller sample sizes are sufficient, as demonstrated empirically in Section 6. Crucially, the dependence on dimensionality is only logarithmic, ensuring the bound degrades gracefully as 𝑑 increases.

5. Dealing with NP-Hardness While the MIP formulation presented in Section 4.2 characterizes the exact optimum of the empirical problem, its worst-case exponential runtime reflects the fundamental NP-hardness of the underlying sparse regression problem. To transition from this computational bottleneck to a strictly polynomial-time approximation, we must move beyond combinatorial search and exploit the geometric structure of the empirical objective. To this end, we rely on the Restricted Strong Convexity (RSC) and Restricted Strong Smoothness (RSS) properties (Negahban, Yu, Wainwright and Ravikumar, 2009). Widely utilized in the statistical learning literature for sparse recovery (Agarwal, Negahban and Wainwright, 2010; Shalev-Shwartz, Srebro and Zhang, 2010; Jalali, Johnson and Ravikumar, 2011; Jain, Tewari and Kar, 2014; Yuan, Li and Zhang, 2017), these conditions ensure that the empirical loss function behaves much like a strongly convex and smooth function, provided the optimization trajectory is restricted to sparse vectors. This section establishes that our parameterized family of neighborhood distributions 𝒙,𝜎 satisfies these RSC and RSS properties with high probability, provided the sample size scales appropriately with the concentration parameter 𝜎. This geometric regularity is what allows us to deploy a variant of the Iterative Hard Thresholding (IHT) algorithm (Blumensath and Davies, 2008, 2009), adapted to the anchoring constraint of our framework and equipped with convergence guarantees. The resulting method runs in polynomial time and offers a scalable alternative to the MIP formulation.

5.1. Restricted Strong Convexity and Smoothness In standard continuous optimization, strong convexity ensures that an objective function curves upward everywhere, allowing gradient-based methods to rapidly converge to a unique global minimum. In contrast, generating 𝑘-sparse linear explanations constrains the feasible space within the non-convex ball 0 (𝑘). When instances are sampled uniformly at random (𝜎 = 0), the data matrix satisfies, with high probability and for 𝑚 large enough, the classic √ Restricted Isometry Property (RIP) after the usual normalization (1∕ 𝑚)𝒁—a foundational concept in compressed F. Koriche et al.: Preprint submitted to Elsevier

Page 16 of 35

Probabilistic Linear Explanations

sensing that guarantees a matrix acts nearly as an isometry on sparse vectors (Candes and Tao, 2005). Unfortunately, when we localize explanations (𝜎 > 0), the sampled instances become biased toward the central instance 𝒙. The individual features remain mutually independent, but they are no longer centered, so that the second-moment matrix ceases to be near-isotropic: it acquires a rank-one component along the single direction 𝒙, which is what violates standard RIP assumptions. The RSC and RSS properties handle this issue by extending the isometry concept to more general loss functions—such as our empirical fidelity error—and, as we will see, the offending direction is precisely the one that the anchoring constraint removes from consideration. Geometrically, these properties guarantee that the loss function’s curvature remains well-behaved (neither excessively flat nor steep) as long as updates are limited to a small number of features. Definition 7 (Orthogonal RSC/RSS). Given a reference instance 𝒙 ∈ {−1, +1}𝑑 , an integer 𝑠 ≥ 1, and two scalars 𝛼, 𝛽 > 0, a matrix 𝒁 ∈ ℝ𝑚×𝑑 is said to satisfy the 𝛼-RSC and 𝛽-RSS properties of order 𝑠 orthogonal to 𝒙 if for any 𝒘 ∈ 0 (𝑠) such that 𝒘 ⋅ 𝒙 = 0, we have: 𝛼‖𝒘‖22 ≤ 𝑚1 ‖𝒁𝒘‖22 ≤ 𝛽‖𝒘‖22 The ratio 𝛽∕𝛼 is called the restricted condition number. Note that 0 (𝑠′ ) ⊆ 0 (𝑠) whenever 𝑠′ ≤ 𝑠, so these properties at a given order automatically hold at every smaller order. To ground these properties in our unified framework, recall that our goal is to minimize the empirical fidelity error (8) subject to the size limit ‖𝒘‖0 ≤ 𝑘 and the anchoring constraint 𝒘 ⋅ 𝒙 = 𝑓 (𝒙). The reader may notice a discrepancy between this anchoring constraint and the condition 𝒘⋅𝒙 = 0 in Definition 7. This distinction is intentional: in gradientbased optimization, we evaluate curvature along update directions. If 𝒘1 and 𝒘2 are two feasible linear explanations, their difference Δ𝒘 = 𝒘1 − 𝒘2 is a 2𝑘-sparse vector that satisfies Δ𝒘 ⋅ 𝒙 = 0. Thus, bounding the Hessian over these zero-anchored sparse directions directly dictates the stability of the algorithm. Because the empirical fidelity error is ̂ = (1∕2𝑚)𝒁 ⊤ 𝒁. The scalars of Definition 7 quadratic, its Hessian is the constant empirical second-moment matrix ∇2 𝖥 ̂ therefore translate into curvature bounds 𝛼∕2 and 𝛽∕2 for 𝖥 along sparse anchored directions, leaving the restricted ̂ for which these curvature condition number unchanged. Section 5.3 works instead with the rescaled objective 𝐿 = 2𝖥, bounds are exactly 𝛼 and 𝛽. The lemma below establishes these bounds, together with a second, auxiliary property. The rank-one shift separating 𝒁 from its centered counterpart is annihilated along anchored directions, but not elsewhere: the residual of the empirical fidelity error carries a constant component, whose interaction with sparse directions is governed by the empirical mean of the centered rows. Property (ii) states that this mean is uniformly small along such directions, at a tolerance 𝜔 left free at this stage. It plays no part in the RSC/RSS geometry itself and will be invoked in Section 5.3. Lemma 5 (Orthogonal RSC/RSS of Neighborhood Distributions). Given a concentration parameter 𝜎 ≥ 0, a confidence parameter 𝛿 ∈ (0, 1], an order 𝑠 ≥ 1, a target condition scalar 𝛾 > 1, and a tolerance 𝜔 > 0, write 𝛾−1 𝜈 = tanh(𝜎∕2), 𝜆 = sech2 (𝜎∕2), 𝜃 = 2(𝛾+1) , and let 𝒁 = 𝒁 − 𝜈𝟏𝒙⊤ be the centered data matrix. Then there exist absolute constants 𝑐0 , 𝑐1 > 0 such that if the sample size satisfies { ( )2 ( ( ) ( )) ( ) ( )} 2𝑠 2𝑐1 𝑒𝑑 4𝑑 4 𝜎 𝑚 ≥ max 𝑐0 𝛾+1 𝑠 ln + ln cosh , ln 𝛾−1 𝑠 𝛿 2 𝛿 𝜔2 then, with probability at least 1 − 𝛿, the following two properties hold simultaneously: (i) for every 𝒘 ∈ 0 (𝑠), 𝛼‖𝒘‖22 ≤ 𝑚1 ‖𝒁𝒘‖22 ≤ 𝛽‖𝒘‖22 , where 𝛼 = 𝜆(1 − 𝜃) and 𝛽 = 𝜆(1 + 𝜃), so that

𝛽∕𝛼 = 3𝛾+1 < 𝛾; 𝛾+3 |1 ⊤ | (ii) sup | 𝑚 𝟏 𝒁𝒖| ≤ 𝜔. | 𝒖∈ (𝑠),‖𝒖‖ =1 | 0

2

In particular, since 𝒁𝒘 = 𝒁𝒘 for every 𝒘 with 𝒘 ⋅ 𝒙 = 0, property (i) implies that 𝒁 satisfies the 𝛼-RSC and 𝛽-RSS properties of order 𝑠 orthogonal to 𝒙.

F. Koriche et al.: Preprint submitted to Elsevier

Page 17 of 35

Probabilistic Linear Explanations

∑ PROOF. Since the Hamming distance decomposes additively across features as 12 ‖𝒙 − 𝒛‖1 = 𝑑𝑗=1 𝟙[𝑥𝑗 ≠ 𝑧𝑗 ], the distribution 𝒙,𝜎 is a product distribution: the coordinates 𝑧1 , … , 𝑧𝑑 are mutually independent, with ℙ[𝑧𝑗 = 𝑥𝑗 ] = 1∕(1 + 𝑒−𝜎 ), hence ( ) ( −𝜎 ) 𝑒 1−𝑒−𝜎 𝔼[𝑧𝑗 ] = 𝑥𝑗 1+𝑒1−𝜎 − 𝑥𝑗 1+𝑒 = 𝑥𝑗 1+𝑒 Var(𝑧𝑗 ) = 1 − 𝜈 2 = 𝜆 −𝜎 −𝜎 = 𝜈𝑥𝑗 , Writing 𝝁 = 𝜈𝒙 for the mean vector, so that 𝒁 = 𝒁 − 𝟏𝝁⊤ , the covariance matrix of 𝒙,𝜎 is perfectly isotropic, while the second-moment matrix carries an additional rank-one term: 𝚺 = 𝜆𝑰,

𝑸 = 𝔼[𝒛𝒛⊤ ] = 𝚺 + 𝜈 2 𝒙𝒙⊤

Locality thus distorts the geometry in two distinct ways: it shrinks the isotropic component by a factor 𝜆, and it adds a bias along the single direction 𝒙. The rows of 𝒁 are i.i.d. copies of 𝒛 − 𝝁; they are centered, have covariance 𝚺, and each of their coordinates lies in an interval of length 2, hence is sub-Gaussian with parameter 1 by Hoeffding’s Lemma. Property (i). By standard restricted eigenvalue bounds for empirical covariance matrices of centered sub-Gaussian vectors (Wainwright, 2019, Theorem 6.5), for any tolerance 𝑡 ∈ (0, 1] the uniform deviation over all 𝑠-sparse unit vectors satisfies sup

| |1 | 𝑚 ‖𝒁𝒘‖22 − 𝒘⊤ 𝚺𝒘| ≤ 𝑡 |

𝒘∈0 (𝑠),‖𝒘‖2 =1 |

with probability at least 1 − 𝑐1 exp(−𝑐2 𝑚𝑡2 + 𝑠 ln(𝑒𝑑∕𝑠)), where 𝑐1 , 𝑐2 are absolute constants; the tolerance we select below is at most 𝜆∕2 ≤ 1∕2, so this quadratic regime indeed applies. Since 𝒘⊤ 𝚺𝒘 = 𝜆‖𝒘‖22 and both sides scale with

‖𝒘‖22 , the bound extends to every 𝒘 ∈ 0 (𝑠) in the form | 𝑚1 ‖𝒁𝒘‖22 − 𝜆‖𝒘‖22 | ≤ 𝑡‖𝒘‖22 . Setting 𝑡 = 𝜃𝜆 yields exactly the scalars 𝛼 = 𝜆(1 − 𝜃) and 𝛽 = 𝜆(1 + 𝜃) announced in the statement, whose ratio is 𝛽 3𝛾 + 1 1+𝜃 = = <𝛾 𝛼 1−𝜃 𝛾 +3 the last inequality being equivalent to 𝛾 2 > 1. The appearance of cosh4 (𝜎∕2) = 1∕𝜆2 in the sample size traces back to a mismatch of scales. The deviation bound above is stated in absolute terms, the rows of 𝒁 being sub-Gaussian with parameter 1 irrespective of 𝜎, whereas the curvature to be estimated is itself of order 𝜆. Securing a relative precision 𝜃 on that curvature thus requires an absolute precision 𝜃𝜆, and the sample size scales with the inverse square of this quantity.4 Formally, requiring 𝑐1 exp(−𝑐2 𝑚𝑡2 + 𝑠 ln(𝑒𝑑∕𝑠)) ≤ 𝛿∕2, taking logarithms, isolating 𝑚, and substituting 𝑡 = 𝜃𝜆 with 1∕𝜆 = cosh2 (𝜎∕2) gives ( ( ) ( ( )) )2 ( ( ) ( )) ( ) 2𝑐1 2𝑐1 𝛾+1 4 𝑒𝑑 4 𝜎 𝑚 ≥ 𝑐 1𝑡2 𝑠 ln 𝑒𝑑 + ln = 𝑠 ln + ln cosh 𝑠 𝛿 𝑐 𝛾−1 𝑠 𝛿 2 2

2

which is the first term of the maximum, with 𝑐0 = 4∕𝑐2 . ⊤

Property (ii). Let 𝜻 = 𝑚1 𝒁 𝟏 be the empirical mean of the centered rows, so that the supremum in (ii) equals the √ largest 𝓁2 norm of any 𝑠 coordinates of 𝜻, itself at most 𝑠 ‖𝜻‖∞ . Each coordinate is an average of 𝑚 i.i.d. centered variables that are √ sub-Gaussian with parameter 1, so Hoeffding’s inequality and a union bound over the 𝑑 coordinates give ‖𝜻‖∞ ≤ 2 ln(4𝑑∕𝛿)∕𝑚 with probability at least 1 − 𝛿∕2. The second term of the maximum is precisely the √ condition ensuring 2𝑠 ln(4𝑑∕𝛿)∕𝑚 ≤ 𝜔. A union bound over the two events concludes the proof. Two consequences of Lemma 5 deserve emphasis, particularly as they formalize the fundamental computational trade-off with the MIP approach of Section 4.2. First, while MIP’s sample complexity is entirely independent of 𝜎, the required number of samples here scales with cosh4 (𝜎∕2). This is a direct statistical consequence of the distribution’s locality: as 𝜎 → ∞, the variance 𝜆 of each individual feature vanishes. Estimating a geometric curvature that is itself flattening out requires proportionally more data to separate the signal from the noise. Fortunately, within the 4 This dependence is not tight: the sub-Gaussian parameter obtained from Hoeffding’s Lemma is loose for strongly biased coordinates, whose optimal proxy also decreases with 𝜎. Sharpening it would complicate the analysis without altering the qualitative picture.

F. Koriche et al.: Preprint submitted to Elsevier

Page 18 of 35

Probabilistic Linear Explanations

Algorithm 1: Iterative Hard Thresholding Explainer Input: explanation query (𝒙, 𝑓 (𝒙)), labeled sample set {(𝒛𝑖 , 𝑓 (𝒛𝑖 ))}𝑚 , sparsity level 𝑘, step-size 𝜂, iteration 𝑖=1 count 𝑇 𝒁 ← (𝒛1 , ⋯ , 𝒛𝑚 ) 𝒚 ← (𝑓 (𝒛1 ), ⋯ , 𝑓 (𝒛𝑚 )) 𝒘0 ← 𝑓 (𝒙) 𝒆1 for 𝑡 = 1, 2, … , 𝑇 do 𝒗𝑡 ← 𝒘𝑡−1 − 𝑚𝜂 𝒁 ⊤ (𝒁𝒘𝑡−1 − 𝒚) 𝒖𝑡 ← 𝒗𝑡 ⊙ 𝒙 𝒖∗𝑡 ← GSHP(𝒖𝑡 , 𝑘, 𝑓 (𝒙)) 𝒘𝑡 ← 𝒖∗𝑡 ⊙ 𝒙

cognitively meaningful range of concentration parameters discussed in Section 4, this overhead remains manageable— the multiplier cosh4 (𝜎∕2) is approximately 1.6 at 𝜎 = 1, 5.7 at 𝜎 = 2, and 31 at 𝜎 = 3—only becoming prohibitive for 𝜎 ≳ 5. Second, and crucially, both 𝛼 and 𝛽 shrink by the exact same factor 𝜆, meaning the restricted condition number 𝛽∕𝛼 remains strictly bounded by 𝛾, uniformly across all 𝜎. In short, locality costs samples, not conditioning— and as the next subsection demonstrates, it is the conditioning that ultimately dictates the convergence guarantees of gradient-based methods. As for property (ii), the tolerance 𝜔 enters the sample size only through √ the second term of the maximum, which is linear in 𝑠 and logarithmic in 𝑑. Section 5.3 calls for 𝜈𝜔 = Θ(𝛽∕ 𝑘) at order 𝑠 = 3𝑘, under which this second term scales as 𝑘2 sinh2 (𝜎) ln(4𝑑∕𝛿). It carries the same asymptotic dependence on 𝜎 as the first term, since sinh2 (𝜎) ∼ 4 cosh4 (𝜎∕2), but a quadratic rather than linear dependence on the order; it vanishes altogether at 𝜎 = 0, where the rows are already centered and the property is void. Both terms are in turn dominated by the sample size required for generalization in Theorem 3, so the auxiliary property does not drive the overall complexity.

5.2. The Iterative Hard Thresholding Explainer With restricted convexity and smoothness ensured by our exponentially localized distributions, we can efficiently approximate the optimal sparse explanation using a tailored variant of the Iterative Hard Thresholding (IHT) algorithm. The classic IHT method (Blumensath and Davies, 2008, 2009) operates by alternating between two steps: a standard gradient descent update to reduce the loss, followed by a nonlinear thresholding operation that retains only the 𝑘 largest-magnitude coefficients, projecting the weights back onto the 𝐿0 ball. However, standard IHT alone is insufficient for our explainability framework. A valid probabilistic linear explanation must not only be sparse but also exactly reconstruct the black-box prediction for the given instance— that is, it must satisfy the anchoring constraint 𝒘 ⋅ 𝒙 = 𝑓 (𝒙). To achieve this, our variant (described in Algorithm 1) incorporates a specialized projection step: after each gradient update, the continuous weights are projected onto the geometric intersection 𝒙,𝑘 = (𝒙, 𝑓 (𝒙)) ∩ 0 (𝑘) of the anchoring hyperplane and the 𝑘-sparse ball. The algorithm is initialized with the trivial anchored explanation 𝒘0 = 𝑓 (𝒙) 𝒆1 , which places the entire prediction on the intercept and attributes nothing to any feature. This 1-sparse vector is already feasible—anchored and inside the unit box—so that every iterate, from 𝒘0 onwards, is a valid 𝑘-sparse anchored explanation. As we show below, this seemingly cosmetic choice is what keeps the iterates bounded, and hence the sample complexity under control. The following result ensures that the projection step inside the loop operates in low polynomial time, circumventing the combinatorial explosion typically associated with sparse constrained optimization. Lemma 6 (Anchored Sparse Projection). Let (𝒙, 𝑓 (𝒙)) ∈ {−1, +1}𝑑 × [−1, +1] be a labeled data instance and 𝑘 ≥ 1 be a sparsity level. Then, the Euclidean projection of any vector 𝒗 ∈ ℝ𝑑 onto 𝒙,𝑘 = (𝒙, 𝑓 (𝒙)) ∩ 0 (𝑘) can be computed exactly in (𝑑 log 𝑑 + 𝑘2 ) time. PROOF. As outlined in the projection phase of Algorithm 1, the objective is to compute 𝒘𝑡 = Π𝒙,𝑘 (𝒗𝑡 ). To leverage existing sparse projection techniques, we perform a change of variables by taking the Hadamard product 𝒖𝑡 = 𝒗𝑡 ⊙ 𝒙. Let  denote the intersection of 0 (𝑘) with the standard diagonal hyperplane (𝟏, 𝑓 (𝒙)). F. Koriche et al.: Preprint submitted to Elsevier

Page 19 of 35

Probabilistic Linear Explanations

Because 𝒙 ∈ {−1, +1}𝑑 , we have 𝒙 ⊙ 𝒙 = 𝟏, meaning the Hadamard product is invertible via self-multiplication (i.e., 𝒖𝑡 ⊙ 𝒙 = 𝒗𝑡 ⊙ 𝒙 ⊙ 𝒙 = 𝒗𝑡 ). Consequently, for any vector 𝒘′ and its transformed counterpart 𝒖′ = 𝒘′ ⊙ 𝒙, we have the algebraic equivalence 𝒖′ ⋅ 𝟏 = 𝑓 (𝒙) if and only if 𝒘′ ⋅ 𝒙 = 𝑓 (𝒙). This implies that 𝒖′ ∈  if and only if 𝒘′ ∈ 𝒙,𝑘 . This structural equivalence, combined with the fact that Euclidean distance is preserved under this sign-flipping bijection (‖𝒖′ − 𝒖𝑡 ‖2 = ‖𝒘′ − 𝒗𝑡 ‖2 ), implies that the target projection can be factored as: ( ) Π𝒙,𝑘 (𝒗𝑡 ) = Π (𝒖𝑡 ) ⊙ 𝒙 Let 𝒖∗𝑡 be the projection of 𝒖𝑡 onto . Since  is exactly the intersection of the 𝑘-sparse set with a diagonal hyperplane, computing 𝒖∗𝑡 is an instance of the sparse hyperplane projection problem of Kyrillidis, Becker, Cevher and Koch (2013), whose Greedy Selector and Hyperplane Projector (GSHP) algorithm solves it exactly (their Theorem 2) in (𝑑 log 𝑑 + 𝑘2 ) time. By setting 𝒘𝑡 = 𝒖∗𝑡 ⊙ 𝒙, we obtain the exact projection of 𝒗𝑡 onto 𝒙,𝑘 within the stated time complexity limit, which concludes the proof.

5.3. Convergence Analysis Two design choices deserve a comment before we turn to the convergence analysis. First, the feasible set 𝒙,𝑘 deliberately omits the box constraint ‖𝒘‖∞ ≤ 1 that appears in the MIP formulation (MIP). The reason is that exactness of the projection—the hypothesis on which the whole convergence argument rests—is available for 0 (𝑘) intersected with a hyperplane, but not for the further intersection with a box: the correctness proof of GSHP relies on a telescoping decomposition of the objective that clipping destroys, and Kyrillidis et al. (2013) themselves show that the naive greedy selector already fails on the hyperplane variant. Second, and independently, the box remains present where it matters: we take as reference optimum { } ̂𝑓 ,𝒙,𝑚 (𝒘) ∶ 𝒘 ∈ 𝒙,𝑘 ∩ [−1, +1]𝑑 ̂ ∈ argmin 𝖥 𝒘 ̂ ∈ 𝒙,𝑘 , which holds; and comparing that is, exactly the optimum targeted by (MIP). The analysis only requires 𝒘 Algorithm 1 against the very same point as the MIP solver is what makes the two approaches directly commensurable. All the results below rest on a single instance of Lemma 5, which we establish as a standing assumption. Assumption 1. The centered matrix 𝒁 satisfies properties (i) and (ii) of Lemma 5 at order 𝑠 = 3𝑘 with target condition √ 17 scalar 𝛾 = 1.25—hence 𝜃 = 1∕18, 𝛼 = 18 3𝑘). 𝜆 and 𝛽 = 19 𝜆—and with a tolerance 𝜔 such that 𝜈 𝜔 ≤ 𝛽∕(16 18 Moreover, Algorithm 1 is run with the step-size 𝜂 = 1∕𝛽. The order 3𝑘 is the largest one the analysis requires: it arises in Step 1 of Lemma 7 below, where the supports of ̂ must be considered jointly. Every other sparse direction appearing in this subsection two consecutive iterates and of 𝒘 is a difference of two anchored 𝑘-sparse explanations, hence 2𝑘-sparse; by the monotonicity noted in Definition 7, Assumption 1 covers these as well, so we never invoke a second order. Lemma 5 instantiated at 𝑠 = 3𝑘 and 𝛾 = 1.25 shows that the assumption holds with probability at least 1−𝛿 as soon as 𝑚 meets the corresponding sample size bound. The resulting restricted condition number 𝛽∕𝛼 = 19∕17 ≈ 1.12 sits comfortably below the threshold 16∕9 ≈ 1.78 that Lemma 7 requires for contraction, and yields the factor 𝜌 = 51∕152 ≈ 0.34. Fixing 𝛾 in this way is what makes the constants of Theorem 3 explicit; any other admissible choice would only change their numerical values. Note also that the prescribed step-size is fully determined by the concentration parameter, since 𝜂 = 1∕𝛽 = 18 2 cosh (𝜎∕2). Unlike most IHT variants, whose step-size must be tuned or estimated from the data, Algorithm 1 19 therefore requires no calibration: 𝜎 is chosen by the user as part of the explanation query. One further point deserves attention. In classical sparse recovery, the target is a stationary point of the unconstrained ̂ minimizes the empirical fidelity objective, so that its gradient vanishes and the iterates converge to it exactly. Here, 𝒘 only within a constrained set, and there is no reason for its gradient to vanish: the model 𝑓 is not assumed to agree with any anchored 𝑘-sparse linear function on the sample. The natural quantity measuring this residual is the restricted gradient norm at the optimum, ̂‖ 𝜉 = max ‖ ‖∇𝑆 𝐿(𝒘) ‖2 , |𝑆|≤3𝑘

1 ̂𝑓 ,𝒙,𝑚 (𝒘) where 𝐿(𝒘) = 2𝑚 ‖𝒁𝒘 − 𝒚‖22 = 2 𝖥

(9)

and ∇𝑆 denotes the gradient restricted to the coordinates in 𝑆. The scaling of 𝐿 is chosen so that ∇𝐿(𝒘) = 1 ⊤ 𝒁 (𝒁𝒘 − 𝒚) is exactly the gradient appearing in the update rule of Algorithm 1; it differs from the normalized 𝑚 F. Koriche et al.: Preprint submitted to Elsevier

Page 20 of 35

Probabilistic Linear Explanations

̂ itself, but to a fidelity error of (8) by the constant factor 2. Accordingly, Algorithm 1 converges linearly not to 𝒘 ̂ whose radius is proportional to 𝜉. The degenerate case 𝜉 = 0—which occurs exactly when 𝒘 ̂ is neighborhood of 𝒘 also an unconstrained stationary point—recovers the classical exact linear convergence. Lemma 7 (Convergence of IHT). Under Assumption 1, the iterates of Algorithm 1 satisfy, for every 𝑡 ≥ 0, ̂ 2 ≤ 𝜌 𝑡 ‖𝒘0 − 𝒘‖ ̂ 2+ ‖𝒘𝑡 − 𝒘‖

2𝜉 , 𝛽(1 − 𝜌)

( ) 𝜌 = 2 1 − 𝛼𝛽 + 18 < 1

̂ and let 𝑆 = support(𝒘𝑡 ) ∪ support(𝒘𝑡−1 ) ∪ support(𝒘), ̂ so that |𝑆| ≤ 3𝑘 and PROOF. Fix 𝑡 ≥ 1, let 𝒉 = 𝒘𝑡−1 − 𝒘, ̂ are both anchored, 𝒉 ⋅ 𝒙 = 0. With 𝒗𝑡 = 𝒘𝑡−1 − 𝜂∇𝐿(𝒘𝑡−1 ), the update rule 𝒉 is supported in 𝑆. Since 𝒘𝑡−1 and 𝒘 reads 𝒘𝑡 = Π𝒙,𝑘 (𝒗𝑡 ). ̂ vanish, so the vectors 𝒘𝑡 − 𝒗𝑡 and 𝒘 ̂ − 𝒗𝑡 agree on Step 1: reduction to the support 𝑆. Outside 𝑆, both 𝒘𝑡 and 𝒘 ̂ ∈ 𝒙,𝑘 and 𝒘𝑡 is the exact projection of 𝒗𝑡 onto 𝒙,𝑘 (Lemma 6), we have ‖𝒘𝑡 − 𝒗𝑡 ‖2 ≤ ‖𝒘 ̂ − 𝒗𝑡 ‖2 , and 𝑆 𝑐 . Since 𝒘 ̂ − 𝒗𝑡 )𝑆 ‖2 . The triangle inequality then gives subtracting the common 𝑆 𝑐 contribution leaves ‖(𝒘𝑡 − 𝒗𝑡 )𝑆 ‖2 ≤ ‖(𝒘 ̂ 2 = ‖(𝒘𝑡 − 𝒘) ̂ 𝑆 ‖2 ≤ 2 ‖(𝒗𝑡 − 𝒘) ̂ 𝑆 ‖2 ‖𝒘𝑡 − 𝒘‖ ⊤

̂ = 1 𝒁 𝒁 and using ∇𝐿(𝒘𝑡−1 ) = ∇𝐿(𝒘) ̂ + 𝑚1 𝒁 ⊤ 𝒁𝒉 Step 2: decomposition of the gradient step. Writing 𝑯 𝑚 together with 𝒁𝒉 = 𝒁𝒉 (valid since 𝒉 ⋅ 𝒙 = 0), we obtain ( ) ( ) ( ) 1 ⊤ ̂ + 𝜈 1 𝟏⊤ 𝒁𝒉 𝒙 ⟹ 𝒗𝑡 − 𝒘 ̂ 𝒉 − 𝜂𝜈 1 𝟏⊤ 𝒁𝒉 𝒙 − 𝜂∇𝐿(𝒘) ̂ = 𝑰 − 𝜂𝑯 ̂ 𝒁 𝒁𝒉 = 𝑯𝒉 𝑚 𝑚 𝑚 The single rank-one leakage term is aligned with 𝒙; it is the only place where the non-centered part of the data matrix survives, and it vanishes identically at 𝜎 = 0. Step 3: bounding the three terms on 𝑆. Since 𝒉 is supported in 𝑆, property (i) at order 3𝑘 states that the eigenvalues ̂ ̂𝑆𝑆 lie in [0, 1 − 𝛼∕𝛽], whence ‖((𝑰 − 𝜂 𝑯)𝒉) ̂ 𝑆 ‖2 ≤ (1 − 𝛼∕𝛽)‖𝒉‖2 . of 𝑯𝑆𝑆 lie in [𝛼, 𝛽]; with 𝜂 = 1∕𝛽, those of 𝑰𝑆 − 𝜂 𝑯 √ √ 1 ⊤ For the leakage term, property (ii) gives | 𝑚 𝟏 𝒁𝒉| ≤ 𝜔‖𝒉‖2 , while ‖𝒙𝑆 ‖2 = |𝑆| ≤ 3𝑘; the assumption √ 1 𝜈𝜔 ≤ 𝛽∕(16 3𝑘) therefore bounds it by 16 ‖𝒉‖2 . The last term is at most 𝜂𝜉 by definition of 𝜉. Combining with Step 1, ) ( 𝛼 ̂ 2 ≤ 𝜌‖𝒘𝑡−1 − 𝒘‖ ̂ 2 + 2𝜉 + 18 ‖𝒘𝑡 − 𝒘‖ , 𝜌 = 2 1 − 𝛽 𝛽 The hypothesis 𝛽 < 16 𝛼, guaranteed by Assumption 1, is exactly 𝜌 < 1. Unrolling this recursion and summing the 9 geometric series yields the stated bound. Two further properties of the trajectory are needed to close the analysis. The first bounds the residual 𝜉 by the optimal empirical fidelity itself, which turns the additive floor of Lemma 7 into a multiplicative approximation factor. The second bounds the magnitude of the iterates, which is what allows the uniform convergence guarantee of Theorem 2 to be applied to them: the iterates are not confined to the unit box, so their 𝐿1 norm—the radius 𝐵 of that theorem—must be controlled. Lemma 8 (Trajectory Bounds and Residual Gradient). Under Assumption 1: √ √ √ ̂𝑓 ,𝒙,𝑚 (𝒘) ̂ ≤ 2 3𝑘; (i) 𝜉 ≤ 2 3𝑘 𝖥 √ ̂𝑓 ,𝒙,𝑚 (𝒘𝑡 ) ≤ 𝖥 ̂𝑓 ,𝒙,𝑚 (𝒘0 ) ≤ 1 and ‖𝒘𝑡 ‖1 ≤ 𝑘 + 4 𝑘∕𝛼. (ii) for every 𝑡 ≥ 0, 𝖥 ̂𝑓 ,𝒙,𝑚 (𝒘). ̂ − 𝒚, so that ∇𝐿(𝒘) ̂ = 𝑚1 𝒁 ⊤ 𝒓 and 𝑚1 ‖𝒓‖22 = 4 𝖥 ̂ Since every entry of 𝒁 lies in PROOF. (i) Let 𝒓 = 𝒁 𝒘

̂ ≤ 1 ‖𝒓‖1 ≤ ( 1 ‖𝒓‖22 )1∕2 by Cauchy–Schwarz; taking the {−1, +1}, each coordinate of the gradient satisfies |∇𝑗 𝐿(𝒘)| 𝑚 𝑚 ̂𝑓 ,𝒙,𝑚 (𝒘) ̂𝑓 ,𝒙,𝑚 (𝒘0 ) ≤ 1, ̂ ≤𝖥 𝓁2 norm over at most 3𝑘 coordinates gives the first inequality. The second follows from 𝖥 ̂ problem and |𝒘0 ⋅ 𝒛𝑖 − 𝑓 (𝒛𝑖 )| ≤ 2, so that each loss term 𝓁(𝒘0 ⋅ 𝒛𝑖 , 𝑓 (𝒛𝑖 )) is at most 1. since 𝒘0 is feasible for 𝒘’s F. Koriche et al.: Preprint submitted to Elsevier

Page 21 of 35

Probabilistic Linear Explanations

(ii) Monotonicity. Fix 𝑡 ≥ 1. The loss 𝐿 is quadratic and 𝒘𝑡 − 𝒘𝑡−1 is a difference of two anchored 𝑘-sparse explanations, hence 2𝑘-sparse and orthogonal to 𝒙; the 𝛽-RSS property of Assumption 1 therefore applies to it and gives the standard descent inequality 𝐿(𝒘𝑡 ) ≤ 𝐿(𝒘𝑡−1 ) + ⟨∇𝐿(𝒘𝑡−1 ), 𝒘𝑡 − 𝒘𝑡−1 ⟩ + 𝛽2 ‖𝒘𝑡 − 𝒘𝑡−1 ‖22 . Completing the ( ) square with 𝜂 = 1∕𝛽 rewrites the right-hand side as 𝐿(𝒘𝑡−1 ) + 𝛽2 ‖𝒘𝑡 − 𝒗𝑡 ‖22 − ‖𝒘𝑡−1 − 𝒗𝑡 ‖22 , and the bracket is non-positive because 𝒘𝑡−1 ∈ 𝒙,𝑘 while 𝒘𝑡 is the exact projection of 𝒗𝑡 onto 𝒙,𝑘 . Hence 𝐿(𝒘𝑡 ) ≤ 𝐿(𝒘𝑡−1 ), and the ̂𝑓 ,𝒙,𝑚 = 𝐿∕2. claim follows by induction since 𝖥 ̂ is likewise 2𝑘-sparse and orthogonal to 𝒙, so the 𝛼-RSC property and the Magnitude. The difference 𝒉 = 𝒘𝑡 − 𝒘 elementary inequality ‖𝒂 − 𝒃‖22 ≤ 2‖𝒂‖22 + 2‖𝒃‖22 applied to the two residuals give ( ) ̂𝑓 ,𝒙,𝑚 (𝒘𝑡 ) + 𝖥 ̂𝑓 ,𝒙,𝑚 (𝒘) ̂ ≤ 16 𝛼‖𝒉‖22 ≤ 𝑚1 ‖𝒁𝒉‖22 ≤ 8 𝖥 ̂ 2 + ‖𝒉‖2 ≤ Therefore ‖𝒘𝑡 ‖2 ≤ ‖𝒘‖

√

√ √ 𝑘 + 4∕ 𝛼, and since 𝒘𝑡 is 𝑘-sparse, ‖𝒘𝑡 ‖1 ≤ 𝑘‖𝒘𝑡 ‖2 .

√ Because 1∕𝛼 = cosh2 (𝜎∕2)∕(1 − 𝜃), part (ii) bounds the 𝐿1 norm of every iterate by 𝐵 = 𝑘 + ( 𝑘 cosh(𝜎∕2)). This is the radius that enters Theorem 2, the box constraint of (MIP) corresponding to the particular case 𝐵 = 𝑘. Applying that result with the larger radius leaves the generalization guarantee intact, and extends it to the “unclipped” IHT iterates. We are now in a position to state the main result of this section. Theorem 3 (End-to-End Guarantee). Let 𝜎 ≥ 0 and 𝑘 ≥ 1, √and suppose Algorithm 1 is run with the step-size 𝜂 = 1∕𝛽 prescribed by Assumption 1. Let 𝜉, 𝜌 and 𝐵 = 𝑘 + 4 𝑘∕𝛼 be as above. For any tolerance 𝜀 ∈ (0, 1] and confidence level 𝛿 ∈ (0, 1], if the sample size satisfies both the condition of Lemma 5 at order 3𝑘 with confidence 𝛿∕2 and 𝑚≥

( )) (𝐵 + 1)4 ( 16 ln(2𝑑) + ln 𝛿4 𝜀2

and if the iteration count satisfies (√ ) 𝑘+1 1 ln , 𝑇 ≥ ln(1∕𝜌) 𝜗

𝜗 = min

{

𝜀 , 2𝜉

√

𝜀 2𝛽

}

with the convention 𝜀∕0 = +∞, then, with probability at least 1 − 𝛿, the returned explanation 𝒘𝑇 satisfies (( ) ) ̂𝑓 ,𝒙,𝑚 (𝒘) ̂ +𝜀 𝖱𝑓 ,𝒙,𝒙,𝜎 (𝒘𝑇 ) ≤ (1 + 𝑒−𝜎 )𝑘 1 + 81𝑘 cosh2 (𝜎∕2) 𝖥 PROOF. By Lemma 5 at order 3𝑘 with confidence 𝛿∕2, the first condition on 𝑚 ensures that Assumption 1 holds with probability at least 1 − 𝛿∕2; we work on that event throughout and account for the remaining 𝛿∕2 below. Note that 𝛼, 𝛽, 𝐵 and 𝜂 are deterministic functions of 𝜎 and 𝑘 alone, so no circularity arises. We now combine the ingredients in turn. First, by Lemma 4, 𝖱𝑓 ,𝒙,𝒙,𝜎 (𝒘𝑇 ) ≤ (1 + 𝑒−𝜎 )𝑘 𝖥𝑓 ,𝒙,𝒙,𝜎 (𝒘𝑇 )

(10)

Second, every iterate satisfies ‖𝒘𝑇 ‖1 ≤ 𝐵 by Lemma 8(ii), so Theorem 2 applied with radius 𝐵, tolerance 𝜀∕2 and confidence 𝛿∕2—which is exactly the second displayed condition on 𝑚—gives, with probability at least 1 − 𝛿∕2, ̂𝑓 ,𝒙,𝑚 (𝒘𝑇 ) + 𝜀 𝖥𝑓 ,𝒙,𝒙,𝜎 (𝒘𝑇 ) ≤ 𝖥 2

(11)

̂𝑓 ,𝒙,𝑚 = 𝐿∕2 is 𝛽 -smooth ̂ 2 . The difference 𝒘𝑇 − 𝒘 ̂ is 2𝑘-sparse and orthogonal to 𝒙, so 𝖥 Third, write Δ𝑇 = ‖𝒘𝑇 − 𝒘‖ 2 ̂ gives along it and a second-order expansion around 𝒘 ̂𝑓 ,𝒙,𝑚 (𝒘𝑇 ) − 𝖥 ̂𝑓 ,𝒙,𝑚 (𝒘) ̂ ≤ 𝖥

⟨1 2

⟩ ̂ 𝒘𝑇 − 𝒘 ̂ + 𝛽4 Δ2𝑇 ≤ 12 𝜉Δ𝑇 + 𝛽4 Δ2𝑇 ∇𝐿(𝒘),

F. Koriche et al.: Preprint submitted to Elsevier

Page 22 of 35

Probabilistic Linear Explanations

̂ is only a constrained optimum. By Lemma 7, Δ𝑇 ≤ 𝜗𝑇 + Δ∞ the first-order term being nonzero precisely because 𝒘 ̂ 2 and Δ∞ = 2𝜉∕(𝛽(1−𝜌)). Applying (𝑎+𝑏)2 ≤ 2𝑎2 +2𝑏2 and separating the two contributions, with 𝜗𝑇 = 𝜌 𝑇 ‖𝒘0 − 𝒘‖ ̂𝑓 ,𝒙,𝑚 (𝒘𝑇 ) − 𝖥 ̂𝑓 ,𝒙,𝑚 (𝒘) ̂ ≤ 1 𝜉Δ∞ + 𝛽 Δ2∞ + 1 𝜉𝜗𝑇 + 𝛽 𝜗2𝑇 𝖥 2 2 2 2 ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟ ≤ 3𝜉 2 ∕(𝛽(1−𝜌)2 )

≤ 𝜀∕2 when 𝜗𝑇 ≤𝜗

√ ̂ 2 ≤ 𝑘, the stated iteration count guarantees 𝜗𝑇 ≤ 𝜗; and Lemma 8(i) bounds the first Since ‖𝒘0 ‖2 ≤ 1 and ‖𝒘‖ 51 19 ̂𝑓 ,𝒙,𝑚 (𝒘)∕(𝛽(1 ̂ 𝜆 and 𝜌 = 152 fixed by Assumption 1, we can bracket by 36𝑘 𝖥 − 𝜌)2 ). Inserting the values 𝛽 = 18

upper-bound the multiplicative factor 36𝑘∕(𝛽(1 − 𝜌)2 ) by 81𝑘∕𝜆 = 81𝑘 cosh2 (𝜎∕2). Chaining with (11) and (10), and taking a union bound over the two events, concludes the proof. The iteration count of Theorem 3 involves the residual √ 𝜉, which is not known before running the algorithm. It can nonetheless be made explicit: Lemma 8(i) gives 𝜉 ≤ 2 3𝑘, so that running Algorithm 1 for ) { (√ √ } 𝑘+1 𝜀 1 𝜀 , 𝜗0 = min ln 𝑇 ≥ √ , ln(1∕𝜌) 𝜗0 2𝛽 4 3𝑘

iterations always suffices. This count depends only on 𝑘, 𝜎 and 𝜀, all of which are available upfront, and each iteration costs (𝑚𝑑 + 𝑑 log 𝑑 + 𝑘2 ) time—dominated by the gradient step for any nontrivial sample size. The guarantee has a simple reading: Algorithm 1 recovers the optimum of the MIP formulation (MIP) up to a multiplicative factor, and recovers it exactly in the limit where the model is genuinely explainable by an anchored 𝑘-sparse linear function on the neighborhood of 𝒙. A multiplicative guarantee of this kind is the expected shape for a polynomial-time method applied to an NP-hard problem; what the analysis adds is that the factor is explicit, and that the number of iterations required is only logarithmic in 1∕𝜀. The multiplier 81𝑘 cosh2 (𝜎∕2) in Theorem 3 is a strict, worst-case analytical bound driven by two structural parameters. The dependence on 𝑘 originates from a single union bound in Lemma 8(i), which pessimistically assumes that all 3𝑘 coordinates of the restricted gradient are simultaneously extremal and that no cancellation occurs between the residual and the Boolean columns of 𝒁. The dependence on cosh2 (𝜎∕2) reflects the vanishing curvature of the data matrix (the 1∕𝜆 term in our RSC bound) as the distribution localizes. In practice, such adversarial alignment is exceedingly rare. As demonstrated in Section 6, the empirical gap between the IHT iterate 𝒘𝑇 and the exact MIP ̂ is orders of magnitude smaller than this worst-case factor suggests. The primary value of the theorem is optimum 𝒘 therefore qualitative and structural: it guarantees that the approximation ratio is finite, explicit, and vanishes as the optimum empirical fidelity approaches zero. Overall, Theorem 3 captures the core operational trade-off discussed in Section 1. Locality enters the sample complexity twice: through the restricted eigenvalue condition of Lemma 5, in cosh4 (𝜎∕2), and through the magnitude of the iterates, in (𝐵 + 1)4 = (𝑘4 + 𝑘2 cosh4 (𝜎∕2)). Notably, the second term only starts to dominate once cosh2 (𝜎∕2) exceeds 𝑘, so that for the concentration parameters of practical interest the sample complexity retains the 𝑘4 behavior of Theorem 2. This gives rise to two distinct regimes of algorithmic effectiveness: • Localized Regime (Large 𝜎): When the distribution is highly concentrated around 𝒙, the relevance-fidelity factor (1 + 𝑒−𝜎 )𝑘 from Lemma 4 approaches 1, aligning the empirical fidelity objective almost perfectly with the true relevance error. In this setting, the exact MIP formulation from Section 4.2 is the method of choice, as the sample complexity of Theorem 2 does not depend on 𝜎 at all—its cost being the worst-case exponential runtime of the solver rather than the number of samples. In contrast, IHT suffers a double theoretical penalty: its sample complexity grows with cosh4 (𝜎∕2), and its approximation guarantee degrades with cosh2 (𝜎∕2), making the gradient-based approach eventually infeasible both statistically and computationally. • Intermediate Regime (Moderate 𝜎): For moderate 𝜎, this factor remains small enough to allow for a practical sample size 𝑚. In this “sweet spot,” the data matrix geometry is regular enough to guarantee linear convergence, and the relevance-fidelity factor remains close enough to 1 for high-quality explanations. In this regime, the IHT explainer emerges as a scalable, strictly polynomial-time alternative, avoiding the worst-case exponential runtime associated with the NP-hard exact MIP formulation. F. Koriche et al.: Preprint submitted to Elsevier

Page 23 of 35

Probabilistic Linear Explanations

Ultimately, this theoretical duality shows that formal explainability cannot rely on a one-size-fits-all approach: the choice between exact combinatorial search and approximate gradient-based methods should be guided by the locality properties of the underlying distribution.

6. Experiments This section compares our two explainers with LIME, MAPLE, and a convex relaxation of (MIP) across 16 OpenML benchmarks in both classification and regression. Section 6.1 details the evaluation protocol. Section 6.2 examines how effectively each explainer minimizes empirical fidelity error and at what cost to admissibility: MIP and IHT yield statistically indistinguishable results, while LIME’s occasional fidelity advantage stems from its systematic violation of the anchoring condition. Section 6.3 addresses the metrics most relevant to users, demonstrating that empirical error reliably predicts out-of-sample performance, and that the ranking reverses when considering relevance—the objective our framework optimizes. Parameter sweeps over sparsity budget, locality, and sample size confirm that this reversal is robust, and that no amount of data can compensate for an unanchored explanation.

6.1. Experimental Setup All experiments were implemented in Python, and executed on a computing node equipped with a 24-core Intel Xeon Gold 6252 processor (2.10 GHz), 384 GB RAM, and an NVIDIA Quadro RTX 8000 GPU (48 GB VRAM). The full source code, environment specifications, and scripts for reproducing these experiments are publicly available.5

Datasets. We evaluate our methods on 16 tabular datasets from the OpenML repository, as summarized in Table 1. These datasets span a diverse range of application domains, including biology (Cancer Drug Response Methylation, Musk, NCI 60 Thioguanine), healthcare (Heart Disease, Medical Charges, Postoperative), environment (Forest Fires, Seoul Bike Sharing, Nomao), sociology (Speed Dating, Student Performance), finance (Adult, California Housing, Credit Approval), and law (COMPAS, Communities and Crime). The datasets are divided into two groups: the first eight correspond to binary classification tasks, while the remaining eight are used for continuous regression. To create the Boolean hypercube representation required by our framework, we encode each attribute as a set of indicator columns, mapping {0, 1} to {−1, +1}. Attributes with at most 𝐾 = 4 distinct values are one-hot encoded regardless of type, so a numerical attribute with few values is left undiscretized. Categorical attributes with more than 𝐾 categories are reduced to the 𝐾 −1 most frequent, with all others grouped into a residual category. Numerical attributes with more than 𝐾 values are discretized into 𝐾 quantile-based bins using a standard 𝐾-bins strategy. We discard constant columns, as well as numerical attributes whose quantile edges are not all distinct, a situation that can occur for highly skewed or zero-inflated variables. This explains why the indicator counts in Table 1 may be smaller than the raw attribute count from the original data. Each discretized numerical attribute contributes exactly 𝐾 indicators, while a retained categorical attribute contributes at most 𝐾; the quantity 𝑑 − (𝐾 ⋅ Num) therefore measures the contribution of the categorical attributes alone, which ranges from two indicators per attribute (COMPAS) to three (Speed Dating). For regression tasks, target variables are rescaled to [−1, +1] to match the assumed codomain. After preprocessing, the average dimensionality is 𝑑 ≈ 220 (median 44), ranging from 𝑑 = 12 (Medical Charges) to 𝑑 = 1652 (Cancer Drug Response Methylation).

Models. For each dataset, we employ a Multi-Layer Perceptron (MLP) implemented using Scikit-Learn as the

black-box model 𝑓 . For classification, we use the default MLPClassifier settings: a single hidden layer of 100 neurons, ReLU activation, the Adam optimizer, and up to 200 training iterations; the model outputs its predicted label, mapped to {−1, +1}. For regression, we use MLPRegressor with three hidden layers (150, 100, and 50 neurons), ReLU activation, the Adam optimizer, and up to 500 training iterations; the real-valued outputs are clipped to [−1, +1] to ensure 𝑓 matches the assumed codomain. Model quality is evaluated using the normalized quadratic loss ̂ 2 , which, in the binary case, coincides with the zero-one error. The resulting 5-fold cross-validation 𝓁(𝑦, 𝑦) ̂ = 41 (𝑦 − 𝑦) estimates are given in the final column of Table 1.

Explanation Tasks. In our experiments, an explanation task is defined by a tuple (𝑓 , 𝒙, 𝑘, 𝑚, 𝜎), where 𝑓 denotes the

neural network trained on the benchmark and 𝒙 is a reference instance sampled uniformly at random from the test set. Sampling from 𝒙,𝜎 is exact and requires no burn-in: conditioned on 𝒙, the coordinates of 𝒛 are mutually independent, 5 Source code and replication scripts: https://github.com/FredericKoriche/ProbabilisticLinearExplanations.git

F. Koriche et al.: Preprint submitted to Elsevier

Page 24 of 35

Probabilistic Linear Explanations Table 1 Overview of the 16 OpenML benchmark datasets utilized in our experiments. Total Raw denotes the number of original attributes in each OpenML file; Cat and Num indicate the counts of categorical and numerical attributes retained after preprocessing, encoded via one-hot encoding and 𝐾 = 4 quantile-based binning, respectively. Dim represents the resulting Boolean dimension. The first eight datasets are classification tasks, and the last eight are regression tasks. The final column shows the neural network’s generalization error, assessed by 5-fold cross-validation using the normalized quadratic loss 𝓁(𝑦, 𝑦) ̂ = 14 (𝑦 − 𝑦) ̂ 2. Dataset

ID

Instances

Cat

Num

Total Raw

Dim (𝑑)

CV Loss

Classification Credit Approval COMPAS Postoperative Heart Disease Adult Nomao Speed Dating Musk

29 42192 40683 43672 179 1486 40536 1116

690 5278 88 1190 48842 34465 8378 6598

4 8 8 6 9 41 58 1

2 1 0 5 2 6 2 160

15 13 8 11 14 118 120 167

17 20 23 37 42 149 180 644

0.145 0.405 0.373 0.111 0.163 0.054 0.169 0.015

Regression Medical Charges California Housing Forest Fires Seoul Bike Sharing Student Performance NCI 60 Thioguanine Communities and Crime Cancer Drug Resp. Meth.

44146 44024 44962 46328 42352 46132 46286 46139

163065 20640 517 8760 395 60 1994 475

0 0 2 5 20 0 1 0

3 8 8 8 4 47 91 413

3 8 12 17 32 48 122 808

12 32 40 46 69 188 367 1652

0.003 0.020 0.005 0.008 0.021 0.027 0.030 0.035

each retaining its reference value 𝑥𝑖 with probability 1∕(1 + 𝑒−𝜎 ) and flipping with probability 𝑝𝜎 = 𝑒−𝜎 ∕(1 + 𝑒−𝜎 ). This interpretation provides an intuitive understanding of the locality parameter: 𝜎 = 0 yields the uniform distribution over the hypercube (𝑝𝜎 = 12 ), 𝜎 = 1 perturbs about a quarter of the coordinates (𝑝𝜎 ≈ 0.27), and 𝜎 = 1.75 perturbs about one coordinate in seven (𝑝𝜎 ≈ 0.15). Our experiments span small sparsity budgets (𝑘 ≤ 8), locality regimes ranging from uniform to sharply concentrated neighborhoods, and sample sizes up to 𝑚 = 105 ; the specific grid for each sweep is detailed in the corresponding subsection. For each configuration, all explainers receive the same sample so that performance differences reflect optimization rather than sampling variability. Finally, because the difficulty of an explanation task depends on the reference instance, we select 10 independent instances 𝒙 per benchmark and, for each configuration (𝑓 , 𝑘, 𝑚, 𝜎), report the mean and dispersion of each explainer’s performance across these instances.

Explainers. To solve the explanation tasks, we employ two explainers: MIP, which solves the mixed-integer

formulation (MIP), and IHT, a scalable Iterative Hard Thresholding algorithm. The MIP formulation is solved using the Gurobi solver with default parameters, imposing a strict 120-second time limit. If the time limit is reached, the solver returns its incumbent solution, which remains admissible but may not be optimal; we therefore report the frequency of such timeouts alongside the affected results. The IHT explainer is implemented in PyTorch, leveraging CUDAaccelerated tensor operations for efficient gradient computations and sparse projections. Consequently, wall-clock comparisons between the two methods should be interpreted with the CPU/GPU asymmetry in mind. For IHT, we use the constant step size 𝜂 = 18 cosh2 (𝜎∕2), which is exactly the inverse restricted smoothness constant 1∕𝛽 prescribed 19 by Assumption 1. The algorithm is run for 5000 iterations, initialized at the admissible explanation 𝒘0 = 𝑓 (𝒙) 𝒆1 , in accordance with our convergence analysis. We compare our MIP and IHT methods with three complementary baselines. The first, MIPrelax , serves as an ablation baseline by implementing the continuous convex relaxation of (MIP). This approach is modeled in CVXPY and solved using the Gurobi backend. MIPrelax substitutes the combinatorial sparsity constraint with the continuous surrogate ‖𝒘‖1 ≤ 𝑘, while directly enforcing the anchoring hyperplane 𝒘 ⋅ 𝒙 = 𝑓 (𝒙) and the variable bounds

F. Koriche et al.: Preprint submitted to Elsevier

Page 25 of 35

Probabilistic Linear Explanations Table 2 ̂ across classification and regression benchmarks (𝑘 = 5, 𝑚 = 5000, 𝜎 = 1.0), ordered by Empirical fidelity error (𝖥) dimensionality (𝑑). The last two columns report methods that do not satisfy the sparsity budget; MIPrelax is a relaxation, and its value is a valid lower bound on the optimum of (MIP) rather than an achievable score. Sparse (‖𝒘‖0 ≤ 𝑘) Dataset

Dense

IHT

MIP

LIME

MIPrelax

MAPLE

Classification Credit Approval COMPAS Postoperative Heart Disease Adult Nomao Speed Dating Musk

0.107 ± 0.010 0.105 ± 0.068 0.141 ± 0.048 0.162 ± 0.029 0.152 ± 0.044 0.226 ± 0.064 0.205 ± 0.070 0.124 ± 0.074

0.107 ± 0.010 0.104 ± 0.068 0.141 ± 0.048 0.162 ± 0.030 0.152 ± 0.044 0.224 ± 0.063 0.203 ± 0.071 0.124 ± 0.074

0.108 ± 0.009 0.094 ± 0.048 0.130 ± 0.027 0.162 ± 0.026 0.144 ± 0.029 0.190 ± 0.039 0.181 ± 0.044 0.118 ± 0.062

0.090 ± 0.005 0.081 ± 0.044 0.108 ± 0.024 0.121 ± 0.015 0.114 ± 0.021 0.124 ± 0.023 0.126 ± 0.032 0.077 ± 0.040

0.284 ± 0.060 0.199 ± 0.179 0.246 ± 0.102 0.225 ± 0.046 0.269 ± 0.066 0.188 ± 0.036 0.179 ± 0.041 0.095 ± 0.035

Regression Medical Charges California Housing Forest Fires Seoul Bike Sharing Student Performance NCI 60 Thioguanine Communities and Crime Cancer Drug Resp. Meth.

0.001 ± 0.000 0.020 ± 0.002 0.005 ± 0.000 0.015 ± 0.009 0.012 ± 0.003 0.022 ± 0.021 0.023 ± 0.006 0.019 ± 0.009

0.001 ± 0.000 0.020 ± 0.001 0.005 ± 0.000 0.015 ± 0.009 0.012 ± 0.003 0.022 ± 0.021 0.023 ± 0.006 0.018 ± 0.008

0.001 ± 0.000 0.019 ± 0.001 0.004 ± 0.000 0.010 ± 0.002 0.011 ± 0.001 0.012 ± 0.001 0.018 ± 0.003 0.009 ± 0.001

0.001 ± 0.000 0.012 ± 0.001 0.003 ± 0.000 0.007 ± 0.002 0.007 ± 0.000 0.007 ± 0.001 0.011 ± 0.001 0.003 ± 0.000

0.013 ± 0.008 0.027 ± 0.003 0.005 ± 0.001 0.014 ± 0.006 0.017 ± 0.004 0.009 ± 0.002 0.014 ± 0.002 0.008 ± 0.001

𝒘 ∈ [−1, +1]𝑑 . The other two baselines are state-of-the-art model-agnostic explainers: LIME (Ribeiro et al., 2016) and MAPLE (Plumb et al., 2018).6 Neither LIME nor MAPLE was developed for our hypothesis space, and as such, neither is guaranteed to produce an admissible explanation: the anchoring condition 𝒘⋅𝒙 = 𝑓 (𝒙) is not inherently enforced by either method, and while LIME permits a hard sparsity constraint ‖𝒘‖0 ≤ 𝑘, MAPLE does not. We intentionally retain both methods at their default algorithmic settings to assess them as originally published. As a result, they optimize a less constrained objective and may sometimes report lower fidelity errors than methods restricted to 𝒙,𝑘 —a benefit that is only meaningful if the explanation is admissible. For each configuration, we therefore report the rate at which each baseline violates the anchoring and sparsity constraints, and interpret its fidelity accordingly.

6.2. Optimization Quality and Admissibility This initial set of experiments focuses on the optimization aspect: given a fixed sample, how effectively does each ̂ and at what cost to admissibility and computation? We set the sparsity explainer minimize the empirical fidelity error 𝖥, budget to 𝑘 = 5, the sample size to 𝑚 = 5000, and the locality parameter to 𝜎 = 1.0, running all five explainers on ten reference instances for each of the 16 benchmarks—resulting in 160 explanation tasks in total. Table 2 presents the resulting fidelity errors, with datasets ordered by dimensionality 𝑑 as in Table 1. To avoid multiple tables, support sizes, anchoring violation rates, and runtimes are discussed in the main text and shown in Figure 3.

Admissibility comes first. Two explainers can be fairly compared on 𝖥̂ only if they search the same solution space, so we begin by reporting where each explainer actually lands. IHT, MIP, and LIME consistently yield exactly five nonzero coefficients across all 160 tasks. By construction, MIPrelax produces dense solutions: its 𝐿1 surrogate disperses weight across the entire hypercube, resulting in supports that range from 13 coefficients on Medical Charges to 1412.6 ± 54.8 on Cancer Drug Resp. Meth. MAPLE is also dense, but much less predictable—its support size varies from a single coefficient on Medical Charges to 598.0 ± 648.7 on Cancer Drug Resp. Meth., with standard deviations often matching the mean. This variability makes it difficult for users to anticipate the length of the returned explanation.

6 Because the official MAPLE repository has not been updated since its original release, we manually ported its source code to Python 3.14 to ensure compatibility with our experimental environment. The port, included in our replication code, contains only compatibility fixes.

F. Koriche et al.: Preprint submitted to Elsevier

Page 26 of 35

Probabilistic Linear Explanations

The anchoring constraint clearly distinguishes the methods. IHT, MIP, and MIPrelax enforce this condition structurally and satisfy it exactly on every task. LIME violates the constraint in all 160 cases: its average gap |𝒘 ⋅ 𝒙 − 𝑓 (𝒙)| ranges from 0.026 on Medical Charges to 0.783 on Nomao, and in 7 out of 80 classification tasks, it exceeds 1.0—meaning the surrogate at 𝒙 is further from 𝑓 (𝒙) than the distance between the two classes. MAPLE violates anchoring on 84% of the tasks; its only clean benchmark, Medical Charges, is also where it collapses to a one-coefficient constant model and where its fidelity error is an order of magnitude higher than all other methods. By contrast, one constraint is rarely active. Across all 160 tasks, the largest coefficient produced by IHT is 0.935, and by MIP is 0.956. For MIP, this indicates that the bound ‖𝒘‖∞ ≤ 1 is simply inactive at the optimum. For IHT, the measurement reveals something more: Section 5.3 deliberately omits the box from the feasible set, retaining it only on ̂ so the iterates are free to leave it—but they never do. This is the empirical counterpart of the reference optimum 𝒘, Lemma 8(ii), which controls the magnitude without imposing it. Only MIPrelax reaches the bound exactly, as expected for an 𝐿1 surrogate that distributes weight across many coordinates. LIME, though unconstrained in this regard, also stays within the box on every task (maximum 0.917), whereas MAPLE exceeds the bound on four occasions, reaching up to 1.270. At these sparsity levels, the box constraint does not distinguish between the explainers; instead, the admissibility question reduces to sparsity and anchoring.

IHT matches the exact solver. With these considerations in mind, Table 2 reveals that our two explainers are nearly indistinguishable in terms of optimization quality. Their fidelity errors match to three decimal places on 12 out of 16 benchmarks. On the remaining four, IHT lags by at most 6% in relative terms (Cancer Drug Resp. Meth., 0.019 versus 0.018)—a difference smaller than the natural variability across reference instances. Furthermore, for three of these four benchmarks, the MIP solver was halted by the time limit, meaning its reported value is merely the incumbent rather than a certified optimum. Thus, the apparent advantage is not even definitive. The MIPrelax column serves a distinct purpose. Any 𝒘 meeting ‖𝒘‖0 ≤ 𝑘 and ‖𝒘‖∞ ≤ 1 automatically satisfies ‖𝒘‖1 ≤ 𝑘, so the relaxation expands the feasible set and its optimum provides a valid lower bound for the optimum of (MIP). This makes the ablation informative: on benchmarks where MIP cannot certify optimality, the two columns bracket the unknown optimum—for example, on Musk, it lies within [0.077, 0.124]. Although this bracket is loose, as expected from an 𝐿1 surrogate at 𝑘 = 5, it demonstrates that neither explainer is leaving significant fidelity unachieved.

LIME’s advantage is an anchoring artifact. LIME reports a lower fidelity error than MIP on several benchmarks—

for example, 0.094 versus 0.104 on COMPAS and 0.144 versus 0.152 on Adult. This result warrants an explanation rather than a disclaimer, and the answer is straightforward. When the MIP is solved to optimality, it achieves the ̂ over 𝒙,𝑘 ; thus, any 𝒘 with a strictly lower score must fall outside that feasible set. Since LIME adheres minimum 𝖥 to the sparsity budget and, as discussed above, remains well within the box constraint, the only remaining explanation is anchoring. This advantage does not indicate that LIME is a better optimizer—it simply means LIME is optimizing a different, less constrained problem. Figure 2 quantifies this trade-off. For each task, we plot LIME’s relative fidelity advantage over MIP, 1 − ̂ LIME )∕𝖥(𝒘 ̂ MIP ), against its anchoring gap. These two quantities are strongly correlated across both task families 𝖥(𝒘 (𝑟 = 0.92), and the scatter cloud passes through the origin: when LIME produces an explanation close to the anchoring hyperplane, its advantage disappears or even turns slightly negative, as predicted by theory. Across the 160 tasks, the advantage is actually negative for 56 of them: despite exploring a larger solution space, LIME fails to outperform the anchored optimum in over a third of the cases. Conversely, the largest advantages correspond to anchoring gaps approaching or exceeding 1.0—that is, to explanations that misstate the very prediction they are supposed to explain. The color scale confirms that this effect is not a consequence of dimensionality: it appears at every scale, and the largest gaps occur on high-dimensional benchmarks where a sparse anchored fit is most challenging. All tasks are shown without filtering; the few points with both a small gap and a sizeable relative advantage originate from Medical Charges and Forest Fires, the two benchmarks where every explainer is already within 5 × 10−3 of a perfect fit, making the ratio of two near-zero errors difficult to interpret meaningfully. Together with the elimination argument above, these measurements quantify the magnitude of the discrepancy and establish that this deviation is what yields LIME’s apparent fidelity advantage.

Scalability. The two explainers we propose differ not in solution quality, but in computational cost. Figure 3 presents wall-clock time versus dimension on a doubly logarithmic scale. IHT shows nearly constant runtime: from 𝑑 = 12 to 𝑑 = 1652, it stays between 0.055 and 0.141 seconds, a result of a fixed number of iterations over tensors whose F. Koriche et al.: Preprint submitted to Elsevier

Page 27 of 35

Probabilistic Linear Explanations

(a) Classification Tasks

(b) Regression Tasks

̂ LIME )∕𝖥(𝒘 ̂ MIP ), is plotted against its anchoring gap |𝒘⋅𝒙−𝑓 (𝒙)|. Figure 2: LIME’s relative fidelity advantage over MIP, 1−𝖥(𝒘 Each point corresponds to a single explanation task, with all 80 tasks from each family shown, and color indicates the dimension 𝑑. The advantage grows with anchoring violation and vanishes as the gap closes, illustrating that LIME’s lower fidelity error arises from leaving the hypothesis space, not from more effective optimization within it.

size scales linearly with 𝑑. In contrast, MIP demonstrates the steepest increase. It is the fastest for small instances (0.088 ± 0.013 seconds on Credit Approval), overtakes LIME and MAPLE around 𝑑 ≈ 40, and then saturates at the time limit for larger dimensions. The certification rate highlights this transition: the solver proves optimality on all ten instances of every benchmark with 𝑑 ≤ 69, but only two of the 60 instances drawn from benchmarks with 𝑑 ≥ 149—both on Nomao, the smallest of them—and none at all beyond 𝑑 = 180. LIME and MAPLE occupy the middle ground, with MAPLE consistently more computationally demanding, reaching 57.7 seconds on Cancer Drug Response Methylation, while IHT completes the same task in just 0.137 seconds. This practical division of labor aligns closely with our theoretical expectations. Below the certification threshold, MIP is both exact and efficient, making it the natural choice. Beyond that threshold, it becomes an uncertified anytime heuristic, while IHT still provides explanations of statistically indistinguishable quality in just a tenth of a second, across all tested dimensions.

(a) Classification Tasks

(b) Regression Tasks

Figure 3: Wall-clock time per explanation task as a function of Boolean dimension 𝑑, displayed on a doubly logarithmic scale. Each point corresponds to a benchmark, with a least-squares fit shown for each explainer. MIPrelax is omitted for clarity. MIP points clustered near the horizontal band at 120 seconds represent instances stopped by the time limit. IHT stands out in that its computational cost remains nearly constant as 𝑑 increases.

F. Koriche et al.: Preprint submitted to Elsevier

Page 28 of 35

Probabilistic Linear Explanations Table 3 Out-of-sample fidelity error (𝖥) and relevance error (𝖱) for the three explainers that respect the sparsity budget (𝑘 = 5, 𝑚 = 5000, 𝜎 = 1.0), ordered by dimensionality. MIPrelax and MAPLE are omitted, as neither returns a 𝑘-sparse explanation. Lower is better in both blocks. Fidelity (𝖥) Dataset

Relevance (𝖱)

IHT

MIP

LIME

IHT

MIP

LIME

Classification Credit Approval COMPAS Postoperative Heart Disease Adult Nomao Speed Dating Musk

0.107 ± 0.009 0.106 ± 0.069 0.140 ± 0.046 0.163 ± 0.029 0.153 ± 0.044 0.228 ± 0.064 0.210 ± 0.072 0.127 ± 0.076

0.107 ± 0.009 0.106 ± 0.068 0.140 ± 0.046 0.162 ± 0.030 0.152 ± 0.043 0.227 ± 0.064 0.208 ± 0.074 0.127 ± 0.075

0.108 ± 0.008 0.095 ± 0.047 0.130 ± 0.026 0.161 ± 0.025 0.144 ± 0.030 0.190 ± 0.039 0.182 ± 0.045 0.119 ± 0.062

0.008 ± 0.017 0.057 ± 0.127 0.070 ± 0.099 0.067 ± 0.050 0.097 ± 0.097 0.228 ± 0.131 0.190 ± 0.159 0.103 ± 0.094

0.007 ± 0.015 0.070 ± 0.130 0.069 ± 0.098 0.070 ± 0.048 0.094 ± 0.095 0.225 ± 0.125 0.196 ± 0.160 0.102 ± 0.094

0.043 ± 0.047 0.137 ± 0.239 0.142 ± 0.207 0.124 ± 0.073 0.136 ± 0.159 0.405 ± 0.244 0.257 ± 0.222 0.133 ± 0.128

Regression Medical Charges California Housing Forest Fires Seoul Bike Sharing Student Performance NCI 60 Thioguanine Communities and Crime Cancer Drug Resp. Meth.

0.001 ± 0.000 0.020 ± 0.002 0.005 ± 0.000 0.015 ± 0.009 0.012 ± 0.003 0.023 ± 0.021 0.024 ± 0.006 0.020 ± 0.010

0.001 ± 0.000 0.020 ± 0.001 0.005 ± 0.000 0.015 ± 0.009 0.012 ± 0.003 0.023 ± 0.022 0.024 ± 0.006 0.020 ± 0.009

0.001 ± 0.000 0.019 ± 0.001 0.004 ± 0.000 0.010 ± 0.002 0.011 ± 0.001 0.013 ± 0.001 0.018 ± 0.003 0.009 ± 0.001

0.001 ± 0.001 0.017 ± 0.004 0.005 ± 0.001 0.022 ± 0.019 0.014 ± 0.006 0.038 ± 0.056 0.032 ± 0.017 0.039 ± 0.028

0.001 ± 0.001 0.017 ± 0.003 0.005 ± 0.001 0.022 ± 0.019 0.014 ± 0.006 0.039 ± 0.056 0.032 ± 0.017 0.038 ± 0.027

0.001 ± 0.001 0.022 ± 0.007 0.006 ± 0.001 0.026 ± 0.019 0.019 ± 0.008 0.052 ± 0.066 0.038 ± 0.017 0.047 ± 0.030

6.3. Generalization: Fidelity and Relevance The previous subsection compared explainers on the empirical objective they all minimize. We now turn to the two questions that matter most to users: does a low empirical error lead to a low error on the neighborhood distribution 𝒙,𝜎 , and does a faithful explanation provide insight into the associated prediction? Table 3 addresses both, under the same configuration as before (𝑘 = 5, 𝑚 = 5000, 𝜎 = 1.0), focusing on the three explainers that respect the sparsity budget; MIPrelax and MAPLE, which return dense explanations, are omitted from this subsection.

The empirical error is a faithful proxy. Across all 16 × 3 entries of the fidelity block, the out-of-sample error 𝖥

̂ in Table 2 by more than 0.005, with the largest difference observed on never deviates from its empirical counterpart 𝖥 Speed Dating. This consistency holds across three orders of magnitude in dimension, from 𝑑 = 12 to 𝑑 = 1652, and the gap shows no discernible trend with increasing 𝑑. Empirical risk minimization over 𝒙,𝑘 is thus already tight at 𝑚 = 5000, so the optimization quality measured in Section 6.2 directly translates to generalization quality; the 𝑚-sweep below shows how much earlier this regime is reached.

Fidelity and relevance do not rank the explainers the same way. LIME keeps its fidelity advantage out of sample, for the reason established in Section 6.2: it optimizes over a strictly larger space. The relevance block reverses the verdict. LIME is worse than MIP on 15 of the 16 benchmarks, and ties on the sixteenth, Medical Charges, where every explainer is already at an error of 0.001. The margin is substantial rather than marginal: LIME’s relevance error is 1.3 times that of MIP at the median, twice as large on Nomao (0.405 against 0.225) and on Postoperative (0.142 against 0.069), and six times as large on Credit Approval (0.043 against 0.007). This reversal is the empirical content of the relevance-to-fidelity lemma of Section 4.1. Fidelity averages the surrogate’s error over the whole neighborhood, so an explanation can be faithful on average while being wrong precisely where the user reads it, namely on the subcube its own support designates. Anchoring is what forbids that failure mode, and the fidelity LIME gains by abandoning it is paid for here.

The relevance bound is conservative, and it needs admissibility. At 𝑘 = 5 and 𝜎 = 1.0, the multiplicative factor (1 + 𝑒−𝜎 )𝑘 from Lemma 4 is 4.79. The measured ratio 𝖱∕𝖥 remains well below this value for both of our explainers: it F. Koriche et al.: Preprint submitted to Elsevier

Page 29 of 35

Probabilistic Linear Explanations

(a) COMPAS (classification, 𝑑 = 20)

(b) Adult (classification, 𝑑 = 42)

(c) California Housing (regression, 𝑑 = 32)

(d) Student Performance (regression, 𝑑 = 69)

Figure 4: Relevance error 𝖱 as a function of sparsity budget 𝑘, with 𝑚 = 5000 and 𝜎 = 1.0, averaged over ten reference instances per benchmark. Top row: classification; bottom row: regression. Note the differing vertical scales between rows.

never exceeds 1.95 for IHT and 1.90 for MIP, and stays below 1 on every classification benchmark. The bound is thus conservative, as expected from a worst-case change-of-measure argument, and a small fidelity error is, in practice, a reliable certificate of relevance at the sparsity budgets explanations can afford. LIME falls outside the scope of Lemma 4, and the measurements confirm this is not a mere technicality. Its ratio rises to 4.00 on NCI 60 Thioguanine, just under the bound, and to 5.22 on Cancer Drug Resp. Meth., exceeding it. On the highest-dimensional benchmark of the study, the guarantee is not just unavailable to LIME—it is invalid. The proof makes the reason clear: relevance measures the error against 𝑓 (𝒙) on the subcube singled out by the support, and anchoring ensures the surrogate matches 𝑓 (𝒙) there. Without anchoring, the conditioned quantity is no longer the one fidelity controls. An unanchored explanation may be faithful on average over the neighborhood, yet say nothing reliable on the subcube its own support selects—and no amount of fidelity can repair this disconnect. Admissibility is not just what makes an explanation legible to the reader; it is what makes the cheaper criterion a true proxy for the more demanding one.

Varying the sparsity budget. Figure 4 varies 𝑘 from 1 to 8 across four representative benchmarks, two from each task family, with 𝑚 and 𝜎 held constant. For both IHT and MIP, the relevance error decreases monotonically throughout this range on all benchmarks. This deserves explicit mention, as it is not guaranteed: while a larger budget allows for a better fit, it also reduces the size of the subcube on which relevance is measured by a factor of 1 + 𝑒−𝜎 per additional coordinate, causing the bound in Lemma 4 to degrade exponentially. Across the budgets a practical explanation can offer, the benefit of increased fit prevails, and the lemma’s exponential warning does not materialize. LIME behaves differently. Its relevance error exceeds both anchored explainers at every budget and on every benchmark, by a factor that stays around 1.3 on California Housing and Student Performance and around 1.5 on F. Koriche et al.: Preprint submitted to Elsevier

Page 30 of 35

Probabilistic Linear Explanations

(a) COMPAS (classification, 𝑑 = 20)

(b) Adult (classification, 𝑑 = 42)

(c) California Housing (regression, 𝑑 = 32)

(d) Student Performance (regression, 𝑑 = 69)

Figure 5: Relevance error 𝖱 as a function of the concentration parameter 𝜎, with 𝑘 = 5 and 𝑚 = 5000, averaged over ten reference instances per benchmark. Top row: classification; bottom row: regression. Note that vertical scales differ between panels.

Adult. On COMPAS, its curve is not even monotonic, and at 𝑘 = 8, its relevance error is roughly four times that of MIP. Thus, increasing the budget does not repair an unanchored explanation and may worsen it: the anchoring gap is not a fixed penalty that a richer hypothesis space eventually eliminates. For our two explainers, the sweep reinforces the findings of Section 6.2: their curves closely track each other at every budget, with IHT slightly lower at the smallest budgets on some benchmarks and both converging as 𝑘 increases. Exact coincidence is not expected, since MIP is optimal for empirical fidelity error while the sweep evaluates relevance; nonetheless, the remaining differences are well within the dispersion across reference instances.

Varying the locality. Figure 5 varies 𝜎 from 0 to 1.75, moving from the uniform distribution to a neighborhood that

perturbs roughly one coordinate in seven, with 𝑘 = 5 and 𝑚 = 5000. For three of the four benchmarks, the relevance error for every explainer decreases steadily with increasing 𝜎: by a factor of about three on Adult and 1.7 on the two regression tasks. A more concentrated neighborhood is easier to track with five coefficients, and the conversion factor (1 + 𝑒−𝜎 )𝑘 from Lemma 4 tightens accordingly, dropping from 32 at 𝜎 = 0 to 2.23 at 𝜎 = 1.75. Guarantee and measurement thus improve together as the explanation becomes more local. COMPAS is the exception; Table 1 suggests why: it has by far the largest cross-validation error in the study, so the model exhibits little local structure for any explainer to capture, and increasing locality reveals none. LIME remains above both anchored explainers at every value of 𝜎 and across all benchmarks, by a factor ranging from about 1.2 to 1.6. The uniform case 𝜎 = 0 is particularly telling: the penalty is not an artifact of a specific locality regime that LIME was not designed for—it is already present, and on Adult, it is widest when the neighborhood covers the entire hypercube. F. Koriche et al.: Preprint submitted to Elsevier

Page 31 of 35

Probabilistic Linear Explanations

(a) COMPAS (classification, 𝑑 = 20)

(b) Adult (classification, 𝑑 = 42)

(c) California Housing (regression, 𝑑 = 32)

(d) Student Performance (regression, 𝑑 = 69)

Figure 6: Relevance error 𝖱 as a function of sample size 𝑚, shown on a logarithmic grid, with 𝑘 = 5 and 𝜎 = 1.0, averaged over ten reference instances per benchmark. Top row: classification; bottom row: regression. Note that vertical scales differ between panels.

Our two explainers, meanwhile, remain indistinguishable throughout the entire range. This serves as a mild but useful validation of Section 5.3: the step size 𝜂 prescribed by Assumption 1 depends on 𝜎 via cosh2 (𝜎∕2), and IHT follows the exact solver at every point on the grid with no additional tuning beyond that formula.

Varying the sample size. Figure 6 sweeps 𝑚 from 102 to 105 with 𝑘 = 5 and 𝜎 = 1.0. For IHT and MIP, the results align with predictions from empirical risk minimization: relevance error drops sharply during the first decade, then levels off, plateauing at 𝑚 = 1000 for California Housing and Student Performance, and at 𝑚 = 5000 for Adult and COMPAS. Beyond this point, estimation error essentially disappears, leaving only the approximation error of a fivecoefficient linear model, which no amount of additional sampling can reduce. This context also clarifies the guarantee 11 in Theorem 3: for Adult, at 𝜀 = 10−2 and 𝛿 = 0.05, the required sample size is on √ the order of 10 —eight orders of 3 4 magnitude above the 10 sufficient in practice. The (𝐵 + 1) term, with 𝐵 = 𝑘 + 4 𝑘∕𝛼, accounts for most of this gap: the analysis pays for a bound holding uniformly over all 𝑘-sparse explanations and the worst-case restricted curvature, neither of which is realized in these instances. LIME’s curves behave quite differently: they remain flat from the outset, are irregular, and on COMPAS and Adult, are not even monotonic—the value at 𝑚 = 105 is no better than at 𝑚 = 500. This contrast underscores the subsection’s central point. For anchored explanations, variance is the limiting factor, which can be reduced with more samples; for LIME, it is the violated constraint, which no amount of data can fix. The comparison at both ends of the grid makes this concrete: on all four benchmarks, IHT and MIP with just 𝑚 = 100 samples already achieve lower relevance error than LIME with 𝑚 = 105 , illustrating a thousandfold increase in data that yields no benefit for LIME.

F. Koriche et al.: Preprint submitted to Elsevier

Page 32 of 35

Probabilistic Linear Explanations

7. Conclusion We have introduced a unified framework for probabilistic explainability that encompasses both binary classification and continuous regression within a single formalism. Our framework produces sparse, anchored linear explanations over the Boolean hypercube, which generalize traditional subset-based approaches by incorporating both the magnitude and direction of each feature’s contribution. Identifying the optimal explanation is NPPP -hard for neural networks, so we address the two sources of this complexity separately. The counting challenge is managed by introducing the fidelity error, an unconditional surrogate whose gap to the relevance error is bounded by (1 + 𝑒−𝜎 )𝑘 over a parameterized family of neighborhood distributions. The combinatorial aspect is tackled either exactly with a Mixed Integer Programming formulation or approximately with an Iterative Hard Thresholding method, whose projection onto the anchoring hyperplane is both exact and efficiently computable. Both methods offer sample complexity guarantees that scale polynomially with the sparsity budget and logarithmically with the dimensionality. Experimentally, the two approaches yield nearly identical solution quality, differing primarily in computational cost: MIP certifies optimality in moderate dimensions, while IHT efficiently scales to the largest benchmarks. Compared to LIME and MAPLE, our results support the core claim of the framework: explanations that depart from the hypothesis space may appear more faithful, but they are less relevant—and no amount of data can compensate for that loss. This work opens several avenues for future research, of which we highlight three.

Richer explanation classes. Linear explanations are effective at capturing the magnitude and direction of individual

feature contributions, but they cannot model interactions between features. A natural next step is to explore sparse, lowdegree polynomial explanations, where a small set of monomials—such as binomials reflecting pairwise correlations— replace the linear terms. In the space {−1, +1}𝑑 , products of coordinates also take values in {−1, +1}, so these polynomials can be interpreted as sparse, anchored linear functions over an expanded set of features. Both (MIP) and Algorithm 1 can be applied directly to this lifted problem. Lemma 4 generalizes as well, with the exponent 𝑘 replaced by 𝑞𝑘 for degree-𝑞 monomials, since fixing the support literals determines all relevant monomials. However, the analysis in Section 5.1 does not transfer straightforwardly: the lifted features introduce dependencies, resulting in a non-isotropic covariance structure and requiring new restricted eigenvalue bounds. Whether these analytical challenges can be overcome without substantially increasing sample complexity for small degrees remains open.

Optimizing relevance directly. Our approach achieves relevance by optimizing fidelity, prompting the question of whether this indirect route is necessary. Specifically, could one minimize the empirical relevance error—the average loss over samples satisfying 𝒘 ⋅ 𝒛 = 𝒘 ⋅ 𝒙—directly? Statistically, this seems feasible: Lemma 4 shows that the conditioning event has probability at least (1 + 𝑒−𝜎 )−𝑘 , so estimation only increases sample size requirements by a constant factor. The main challenge, however, is algorithmic. The conditioning event depends discontinuously on 𝒘, making the empirical objective piecewise constant across a complex arrangement of hyperplanes—devoid of gradients, and unsuitable for methods relying on convexity or thresholding. An exact formulation is equally problematic, requiring an indicator for each sample tied to an equality constraint, which in turn renders the objective a ratio of linear forms in these indicators. Underlying all of this is the NPPP -hardness established in Theorem 1, which persists even with empirical relaxation. While the surrogate fidelity-based approach appears more promising in practice, discovering a tractable formulation for the direct problem—even under strong restrictions on 𝑓 —would significantly advance our understanding.

Beyond binarized domains. Our framework currently operates on the Boolean hypercube, mapping categorical

and numerical attributes via the binarization process outlined in Section 6.1. While this approach is common in model-agnostic explainability, it creates a disconnect between the sparsity motivated by cognitive considerations and user reasoning: the sparsity budget 𝑘 counts indicator literals, whereas users typically think in terms of attributes. ∏ Bridging this gap would require extending the framework to a product space 𝑖 Ω𝑖 of finite attribute domains, allowing explanations to be expressed at the attribute level. The core analytical tools—relying only on the product structure of 𝒙,𝜎 and the additive decomposition of Hamming distance—remain applicable under this generalization. The key modification lies in the sparsity constraint, which becomes a group constraint over the one-hot encoded block of each attribute, necessitating a group analogue of the GSHP projection step. A further, more ambitious generalization would be to relax the assumption of coordinate independence, enabling the analysis to follow the data manifold rather than rely on metric balls around 𝒙. However, this would break the covariance isotropy underpinning Section 5.1, requiring a fundamentally new restricted eigenvalue analysis adapted to the resulting dependence structure. F. Koriche et al.: Preprint submitted to Elsevier

Page 33 of 35

Probabilistic Linear Explanations

Acknowledgments. This work has benefited from the support of the AI Chair EXPEKCTATION (ANR-19-CHIA0005-01) of the French National Research Agency. It was also partially supported by TAILOR, a project funded by EU Horizon 2020 research and innovation programme under GA No 952215.

References Agarwal, A., Negahban, S., Wainwright, M.J., 2010. Fast global convergence rates of gradient methods for high-dimensional statistical recovery, in: Advances in Neural Information Processing Systems 23: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS). Agarwal, S., Jabbari, S., Agarwal, C., Upadhyay, S., Wu, S., Lakkaraju, H., 2021. Towards the unification and robustness of perturbation and gradient based explanations, in: Proceedings of the 38th International Conference on Machine Learning (ICML), pp. 110–119. Alvarez-Melis, D., Jaakkola, T., 2018. On the robustness of interpretability methods, in: Proceedings of the 2018 ICML Workshop on Human Interpretability in Machine Learning (WHI). Arenas, M., Barceló, P., Orth, M.A.R., Subercaseaux, B., 2022. On computing probabilistic explanations for decision trees,, in: Advances in Neural Information Processing Systems 35: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS). Audemard, G., Bellart, S., Bounia, L., Koriche, F., Lagniez, J., Marquis, P., 2022. Trading complexity for sparsity in random forest explanations, in: Annual AAAI Conference on Artificial Intelligence, pp. 5461–5469. Barceló, P., Monet, M., Pérez, J., Subercaseaux, B., 2020. Model interpretability through the lens of computational complexity, in: Advances in Neural Information Processing Systems 33: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS). Bartlett, P.L., Jordan, M.I., McAuliffe, J.D., 2006. Convexity, classification, and risk bounds. Journal of the American Statistical Association 101, 138–156. Bassan, S., Huang, X., Katz, G., 2026. Unifying formal explanations: A complexity-theoretic perspective, in: International Conference on Learning Representations (ICLR). Beale, E.M.L., Kendall, M.G., Mann, D.W., 1967. The discarding of variables in multivariate analysis. Biometrika 54, 357–366. Bertsimas, D., King, A., Mazumder, R., 2016. Best subset selection via a modern optimization lens. The Annals of Statistics 44, 813–852. Bertsimas, D., Weismantel, R., 2005. Optimization Over Integers. Dynamic Ideas. Blanc, G., Lange, J., Tan, L., 2021. Provably efficient, succinct, and precise explanations, in: Advances in Neural Information Processing Systems 34: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS), pp. 6129–6141. Blumensath, T., Davies, M.E., 2008. Iterative thresholding for sparse approximations. Journal of Fourier Analysis and Applications 14, 629–654. Blumensath, T., Davies, M.E., 2009. Iterative hard thresholding for compressed sensing. Applied and Computational Harmonic Analysis 27, 265–274. Bounia, L., Koriche, F., 2023. Approximating probabilistic explanations via supermodular minimization, in: Proceedings of the 39th International Conference on Uncertainty in Artificial Intelligence (UAI), pp. 216–225. Burkart, N., Huber, M.F., 2021. A survey on the explainability of supervised machine learning. Journal of Artificial Intelligence Research 70, 245–317. Candes, E., Tao, T., 2005. Decoding by linear programming. IEEE Transactions on Information Theory 51, 4203–4215. Cooper, M., Marques-Silva, J., 2023. Tractability of explaining classifier decisions. Artificial Intelligence 316, 103841. Darwiche, A., Hirth, A., 2020. On the reasons behind decisions, in: European Conference on Artificial Intelligence (ECAI), pp. 712–720. Fligner, M., Verducci, J., 1993. Probability models and statistical analyses for ranking data. volume 80. Springer. Garreau, D., von Luxburg, U., 2020. Explaining the explainer: A first theoretical analysis of LIME, in: Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 1287–1296. Ghorbani, A., Abid, A., Zou, J., 2019. Interpretation of neural networks is fragile, in: Proceedings of the 33rd Conference on Artificial Intelligence (AAAI), pp. 3681–3688. Guidotti, R., 2024. Counterfactual explanations and how to find them: literature review and benchmarking. Data Mining and Knowledge Discovery 38, 2770–2824. Halford, G., Wilson, W., Phillips, S., 1998. Relational complexity: A measure of capacity limitations in processing. Behavioral and Brain Sciences 21, 803–831. Hocking, R.R., Leslie, R.N., 1967. Selection of the best subset in regression analysis. Technometrics 9, 531–540. Hui, L., Belkin, M., 2021. Evaluation of neural architectures trained with square loss vs cross-entropy in classification tasks, in: Proceedings of the 9th International Conference on Learning Representations (ICLR). Ignatiev, A., 2020. Towards trustable explainable AI, in: Proceedings of the 29th International Joint Conference on Artificial Intelligence (IJCAI), pp. 5154–5158. Ignatiev, A., Cooper, M., Siala, M., Hebrard, E., Marques-Silva, J., 2020. Towards formal fairness in machine learning, in: International Conference on Principles and Practiceof Constraint Programming (CP), p. 846–867. Ignatiev, A., Narodytska, N., Marques-Silva, J., 2019. Abduction-based explanations for machine learning models, in: Annual AAAI Conference on Artificial Intelligence, pp. 1511–1519. Izza, Y., Huang, X., Ignatiev, A., Narodytska, N., Cooper, M., Marques-Silva, J., 2023. On computing probabilistic abductive explanations. International Journal of Approximate Reasoning 159, 108939. Jain, P., Tewari, A., Kar, P., 2014. On iterative hard thresholding methods for high-dimensional m-estimation, in: Advances in Neural Information Processing Systems 27: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS), pp. 685–693. Jalali, A., Johnson, C., Ravikumar, P., 2011. On learning discrete graphical models using greedy methods, in: Advances in Neural Information Processing Systems 24: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS), pp. 1935–1943. Johnson-Laird, P., 2010. Mental models and human reasoning. Proceedings of the National Academy of Sciences 107, 18243–18250.

F. Koriche et al.: Preprint submitted to Elsevier

Page 34 of 35

Probabilistic Linear Explanations Kakade, S.M., Sridharan, K., Tewari, A., 2008. On the complexity of linear prediction: Risk bounds, margin bounds, and regularization, in: Advances in Neural Information Processing Systems 21: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS), pp. 793–800. Koriche, F., Lagniez, J.M., Mengel, S., Tran, C., 2024. Learning model agnostic explanations via constraint programming, in: Machine Learning and Knowledge Discovery in Databases. Research Track (ECML/PKDD), pp. 437–453. Kozachinskiy, A., 2023. Inapproximability of sufficient reasons for decision trees. CoRR abs/2304.02781. URL: https://doi.org/10.48550/ arXiv.2304.02781, arXiv:2304.02781. Kyrillidis, A., Becker, S., Cevher, V., Koch, C., 2013. Sparse projections onto the simplex, in: Proceedings of the 30th International Conference on Machine Learning (ICML), pp. 235–243. Lage, I., Chen, E., He, J., Narayanan, M., Kim, B., Gershman, S.J., Doshi-Velez, F., 2019. Human evaluation of models built for interpretability, in: Annual AAAI Conference on Human Computation and Crowdsourcing (HCOMP), pp. 59–67. Lakkaraju, H., Kamar, E., Caruana, R., Leskovec, J., 2019. Faithful and customizable explanations of black box models, in: Proceedings of the 2019 AAAI/ACM Conference on AI, Ethics, and Society (AIES), pp. 131–138. Li, J., Nagarajan, V., Plumb, G., Talwalkar, A., 2021. A learning theoretic perspective on local explainability, in: Proceedings of the 9th International Conference on Learning Representations (ICLR). Littman, M., Goldsmith, J., Mundhenk, M., 1998. The computational complexity of probabilistic planning. Journal of Artificial Intelligence Research 9, 1–36. Lundberg, S.M., Lee, S., 2017. A unified approach to interpreting model predictions, in: Advances in Neural Information Processing Systems 30: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS), pp. 4765–4774. Mallows, C.L., 1957. Non-null ranking models. i. Biometrika 44, 114–130. Marden, J., 1996. Analyzing and modeling rank data. CRC Press. Marques-Silva, J., Ignatiev, A., 2022. Delivering trustworthy AI through formal XAI. Proceedings of the 36th Annual AAAI Conference on Artificial Intelligence , 12342–12350. Miller, G.A., 1956. The magical number seven, plus or minus two: Some limits on our capacity for processing information. Psychological Review 63, 81–97. Molnar, C., 2025. Interpretable Machine Learning: A Guide for Making Black Box Models Explainable. 3rd ed., leanpub.com. Mukherjee, A., Basu, A., 2017. Lower bounds over Boolean inputs for deep neural networks with ReLU gates. CoRR abs/1711.03073. arXiv:1711.03073. Narayanan, M., Chen, E., He, J., Kim, B., Gershman, S., Doshi-Velez, F., 2018. How do humans understand explanations from machine learning systems? An evaluation of the human-interpretability of explanation. CoRR abs/1802.00682. arXiv:1802.00682. Natarajan, B.K., 1995. Sparse approximate solutions to linear systems. SIAM J. Comput. 24, 227–234. Negahban, S., Yu, B., Wainwright, M., Ravikumar, P., 2009. A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers, in: Advances in Neural Information Processing Systems: Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS). Ordyniak, S., Paesani, G., Szeider, S., 2023. The parameterized complexity of finding concise local explanations, in: International Joint Conference on Artificial Intelligence (IJCAI), pp. 3312–3320. Parberry, I., 1996. Circuit complexity and feedforward neural networks, in: Smolensky, P., Mozer, M., Rumelhart, D. (Eds.), Mathematical Perspectives on Neural Networks, pp. 85–111. Plumb, G., Molitor, D., Talwalkar, A., 2018. Model agnostic supervised local explanations, in: Advances in Neural Information Processing Systems 31. Proceedings of the Annual Conference on Neural Information Processing Systems (NeurIPS), pp. 2520–2529. Ribeiro, M.T., Singh, S., Guestrin, C., 2016. "why should I trust you?": Explaining the predictions of any classifier, in: ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 1135–1144. Ribeiro, M.T., Singh, S., Guestrin, C., 2018. Anchors: High-precision model-agnostic explanations, in: Annual AAAI Conference on Artificial Intelligence (AAAI), pp. 1527–1535. Rifkin, R., Klautau, A., 2004. In defense of one-vs-all classification. Journal of Machine Learning Research 5, 101–141. Shalev-Shwartz, S., Ben-David, S., 2014. Understanding Machine Learning: From Theory to Algorithms. Cambridge University Press. Shalev-Shwartz, S., Srebro, N., Zhang, T., 2010. Trading accuracy for sparsity in optimization problems with sparsity constraints. SIAM Journal of Optimization 20, 2807–2832. Slack, D., Hilgard, S., Jia, E., Singh, S., Lakkaraju, H., 2020. Fooling LIME and SHAP: Adversarial attacks on post hoc explanation methods, in: Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society (AIES), p. 180–186. Subercaseaux, B., Arenas, M., Meel, K.S., 2025. Probabilistic explanations for linear models, in: Annual AAAI Conference on Artificial Intelligence, pp. 20655–20662. Suykens, J.A., Vandewalle, J., 1999. Least squares support vector machine classifiers. Neural processing letters 9, 293–300. Wainwright, M., 2019. High-Dimensional Statistics: A Non-Asymptotic Viewpoint. Cambridge University Press. Wäldchen, S., MacDonald, J., Hauch, S., Kutyniok, G., 2021. The computational complexity of understanding binary classifier decisions. Journal of Artificial Intelligence Research 70, 351–387. Yuan, X., Li, P., Zhang, T., 2017. Gradient hard thresholding pursuit. Journal of Machince Learning Research 18, 166:1–166:43. Zhao, X., Huang, W., Huang, X., Robu, V., Flynn, D., 2021. BayLIME: Bayesian local interpretable model-agnostic explanations, in: Proceedings of the 37th International Conference on Uncertainty in Artificial Intelligence (UAI), pp. 887–896.

F. Koriche et al.: Preprint submitted to Elsevier

Page 35 of 35