1 Introduction↩︎

Most influential sets (MIS) are subsets of the training data whose removal induces the largest change in a target quantity of interest, such as a regression coefficient, treatment effect, or prediction. Formally, for a target functional \(\phi\) and a leave-set-out estimator, the size-\(k\) MIS, \({\mathbb{S}}_k^{\max}\), maximizes the leave-set-out discrepancy over all subsets \({\mathbb{S}}\) with \(\lvert {\mathbb{S}}\rvert = k.\) Set influence differs qualitatively from singleton influence because observations interact; some points reinforce each other, while others mask each others’ effects. MIS therefore help diagnose which subsamples drive inferences and how models amplify or suppress particular groups.

In applied work, MIS provide target-specific diagnostics for interpretability, accountability, fairness, robustness, and data curation [1]–[8].

Despite their relevance, the study of MIS has historically been constrained by computation. Exact identification requires maximizing over \(\binom{n}{k}\) subsets, which is computationally infeasible, even for moderate datasets. Classical influence diagnostics focus on individual data points cook1979Influential?, and practical extensions to sets are largely insensitive to higher-order interactions [9]. Consequently, the literature has often favored robust estimation — controlling rather than identifying worst-case sensitivity huber_robust_2009?.

Recent work on influence functions — infinitesimal approximations that are widely used as diagnostic tools [10], [11] — has renewed interest in MIS. This has yielded influence-function-based approximations [12], greedy heuristics and failure modes [13]–[15], influence bounds [16]–[18], and formal tests for excessive influence [19]. However, influence functions remain inaccurate for extreme influence and for sets of data points [20], [21], and efficient and accurate methods for MIS selection have remained unavailable.

In this paper, we develop an efficient algorithm for size-\(k\) MIS selection for influence targets whose leave-set-out effects admit a ratio representation. The key is that these effects can be written as linear-fractional functions of the removed set, transforming the combinatorial search into a one-dimensional parametric optimization problem. Using [22]’s method, each subproblem reduces to selecting the top-\(k\) scores, giving \(\mathcal{O}(n)\) cost per iteration and finite termination. For fixed residualized inputs, the algorithm returns a globally optimal set for the corresponding univariate ratio objective. In partial linear models, this gives exact oracle MIS selection when the nuisance functions are known. With estimated nuisance functions, Neyman orthogonality supports first-stage stability; a separation condition yields consistent selection of the oracle MIS.

1.1 Contributions↩︎

We make four contributions.

• Reduction. We identify a class of leave-set-out influence objectives with linear-fractional structure and show that size-\(k\) MIS selection for these reduces to a one-parameter sequence of top-\(k\) selection problems.

• Algorithm. We apply Dinkelbach’s method to obtain a finite-step algorithm with \(\mathcal{O}(n)\) cost per iteration. For fixed residualized inputs, the algorithm returns a globally optimal set for the univariate ratio objective.

• Theory. The algorithm is exact for any fixed residualized inputs. In partial linear models, the scaled objective uniformly approximates a first-order oracle orthogonal-score objective under generated-score and denominator stability; value consistency and exact set recovery follow under a separation condition.

• Evidence. We benchmark the method against enumeration and greedy baselines and apply it to randomized experiments and datasets from machine learning, biology, economics, and statistics, recovering MIS at scales where enumeration is infeasible.

1.2 Outline↩︎

The remainder of this paper is organized as follows. 2 introduces the partial linear setting and formalizes the MIS problem. 3 develops the linear-fractional reduction, presents the algorithm, and states the theoretical guarantees. 4 evaluates the method in simulations and applications. 5 discusses relevance, limitations, and implications, and 6 concludes.¹

2 Partial Linear Setting↩︎

We study influential sets using residualization in partial linear models [23], [24]. This section introduces the model and estimator, then formalizes influence and \(k\)-most influential sets.

2.1 Model and Estimation↩︎

Consider the PLM \[\label{eq:plm} y_i = x_i \beta_0 + g_0(Z_i) + u_i, \quad \mathbb{E}[u_i \!\mid\! x_i, Z_i] = 0,\tag{1}\] where \((y_i, x_i, Z_i)_{i=1}^n\) are i.i.d.draws from \(P,\) \(x_i\in\mathbb{R}\) is the treatment, \(\beta_0 \in \mathbb{R}\) is the parameter of interest, and \(g_0 : \mathbb{R}^d \to \mathbb{R}\) is an unknown nuisance function of covariates \(Z_i \in \mathbb{R}^d\). Assume \(\mathbb{E}[u_i^2 \!\mid\! Z_i] = \sigma^2(Z_i) < \infty.\)

Let \(h_0(Z_i) \mathrel{\vcenter{:}}= \mathbb{E}[ x_i \!\mid\! Z_i ],\) write \(x_i = h_0(Z_i) + v_i,\) and define \(m_0(Z_i) \mathrel{\vcenter{:}}= \mathbb{E}[y_i\!\mid\! Z_i] = g_0(Z_i) + \beta_0 h_0(Z_i).\)

We estimate \(\beta_0\) using residualized outcomes and treatments: \[\label{eq:residuals} \tilde{y}_i \mathrel{\vcenter{:}}= y_i - \hat{m}(Z_i), \qquad \tilde{x}_i \mathrel{\vcenter{:}}= x_i - \hat{h}(Z_i),\tag{2}\] where \((\hat{m}, \hat{h})\) are flexible cross-fitting first-stage estimators. The residualized (Robinson) second-stage estimator is the ordinary least squares (OLS) coefficient from regressing \(\tilde{y}_i\) on \(\tilde{x}_i,\) \[\label{eq:beta95hat} \hat{\beta} \mathrel{\vcenter{:}}= \arg\min_{\beta \in \mathbb{R}} \sum_{i=1}^n\left( \tilde{y}_i-\tilde{x}_i\beta \right)^2 = \frac{ \sum_{i=1}^n \tilde{x}_i \tilde{y}_i }{ \sum_{i=1}^n \tilde{x}_i^2 }.\tag{3}\] This estimator exists when \(\sum_{i=1}^n \tilde{x}_i^2 > 0;\) otherwise a ridge-stabilized version can be used.

2.2 Influence and Most Influential Sets↩︎

Let \([n]\mathrel{\vcenter{:}}= \{1, \dots, n\}\). For any index set \({\mathbb{S}}\subseteq [n]\), let \(\hat{\beta}_{-{\mathbb{S}}}\) denote the residualized estimator after dropping \({\mathbb{S}}\) from the second stage and recomputing on \([n] \setminus {\mathbb{S}}.\)

The target functional determines the direction and scale of change. We focus on signed affine targets, \(\phi(\beta) = a \beta + b\). Because \(b\) cancels and non-zero \(a\) only rescales, and possibly reverses, the ordering of sets, we set \(a = 1\) in the exposition. Targets with \(a<0\) correspond to directed decreases, while predicted-value targets at treatment level \(x_0\) correspond to \(a = x_0.\)

The influence of \({\mathbb{S}}\) on \(\hat{\beta}\) in 3 admits a ratio form.

This is an exact finite-sample deletion identity, and not an infinitesimal influence-function approximation. It separates two mechanisms; the numerator \[W({\mathbb{S}}) \mathrel{\vcenter{:}}= \sum_{s\in{\mathbb{S}}}\tilde{x}_s\tilde{r}_s\] is the score removed from the full-sample normal equation, while the denominator \[G({\mathbb{S}}) \mathrel{\vcenter{:}}= \sum_{i\notin{\mathbb{S}}}\tilde{x}_i^2\] is the residual curvature that remains after removal. Thus, a set can be influential either because it removes observations with large aligned scores, or because it removes high-curvature observations that amplify the effect of the remaining score imbalance.

Although \(W({\mathbb{S}})\) and \(G({\mathbb{S}})\) are additive separately, their ratio is not. For \(i \notin {\mathbb{S}}\), \[\frac{W({\mathbb{S}})+w_i}{G({\mathbb{S}})-c_i} - \frac{W({\mathbb{S}})}{G({\mathbb{S}})} = \frac{ w_iG({\mathbb{S}})+c_iW({\mathbb{S}}) }{ G({\mathbb{S}})\{G({\mathbb{S}})-c_i\} },\] where \(w_i = \tilde{x}_i \tilde{r}_i\) and \(c_i = \tilde{x}_i^2\). The first term reflects \(i\)’s direct score contribution; the second depends on the score already accumulated by \({\mathbb{S}}\). These interactions generate two phenomena. First, influence can be joint: observations with limited singleton influence may become highly influential when removed together. Second, influence can be masked: observations may appear unimportant in singleton rankings because their effect emerges only after other observations are removed. We use these terms informally in the main text; formal definitions of joint influence, masking, and greedy failure are given in 9.

2.3 Limitations of Existing Approaches↩︎

Standard methods for identifying MIS face obstacles:²

• Enumeration is computationally prohibitive. Exact search requires evaluating \(\binom{n}{k}\) sets; for example, evaluating \(\binom{100}{10}\) sets would take roughly 200 days at \(1\mu\)s per set, and \(k=11\) would increase this to roughly \(4.5\) years.

• First-order approximations, including influence functions [12], are efficient but rank observations using local or singleton information, and can therefore miss joint influence and masking effects.

• Greedy selection [15] builds a set sequentially by adding the observation with the largest current marginal gain. This implicitly assumes nested influential sets across \(k\). Here, however, marginal gains depend on the current set through both \(W({\mathbb{S}})\) and \(G({\mathbb{S}})\), and early choices can exclude the optimal size-\(k\) set.

3 Method↩︎

Finding MIS requires balancing additive score contributions in the numerator against amplification from the shrinking denominator in 1. Conditional on fixed residualized inputs, this trade-off yields a linear-fractional optimization problem over subsets.

3.1 Linear-fractional Reduction↩︎

Map the influence components in 1 to weights and costs, \[w_i \mathrel{\vcenter{:}}= \tilde{x}_i \tilde{r}_i, \qquad c_i \mathrel{\vcenter{:}}= \tilde{x}_i^2, \qquad T \mathrel{\vcenter{:}}= \sum_{i=1}^n c_i.\] For any set \({\mathbb{S}}\) of size \(k\), \[\hat{\beta} - \hat{\beta}_{-{\mathbb{S}}} = \frac{W({\mathbb{S}})}{G({\mathbb{S}})},\] where \(W({\mathbb{S}}) \mathrel{\vcenter{:}}= \sum_{i\in{\mathbb{S}}} w_i,\) and \(G({\mathbb{S}}) \mathrel{\vcenter{:}}= T - \sum_{i\in{\mathbb{S}}}c_i.\)

1 ensures that the estimator and fractional objective are well defined on the feasible class. If the condition fails for some sets, the same reduction applies to a ridge-stabilized denominator, replacing \(G({\mathbb{S}})\) by \(G({\mathbb{S}}) + \lambda\) for \(\lambda > 0\).

3.2 Algorithm and Complexity↩︎

To maximize the fixed-input ratio over size-\(k\) sets, we use [22]’s method. For an auxiliary parameter \(\eta\), define \[\label{eq:dinkelbach95problem} F_\eta({\mathbb{S}}) \mathrel{\vcenter{:}}= W({\mathbb{S}}) - \eta\, G({\mathbb{S}})\tag{4}\] For fixed \(\eta\), \[F_\eta({\mathbb{S}}) = \sum_{i\in{\mathbb{S}}} (w_i + \eta c_i) - \eta T,\] so the term \(-\eta T\) is constant in \({\mathbb{S}}\). Maximizing \(F_\eta({\mathbb{S}})\) over \(|{\mathbb{S}}| = k\) is therefore equivalent to selecting the \(k\) largest scores \[s_i(\eta) \mathrel{\vcenter{:}}= w_i + \eta c_i.\] 2 alternates between this top-\(k\) selection and the ratio update \[\eta^{(t+1)} = \frac{ W(S^{(t+1)}) }{ G(S^{(t+1)}) }.\] With exact arithmetic, finite termination occurs when the selected set no longer changes, equivalently when the Dinkelbach residual is zero. The tolerance \(\tau\) is used as a numerical stopping rule.

3.2.0.1 Warm-start

At the optimum, the Dinkelbach parameter equals the maximal fixed-input influence, \[\eta^\star = \max_{|{\mathbb{S}}|=k} \frac{W({\mathbb{S}})}{G({\mathbb{S}})}.\] This identity can be used to reduce iterations when computing MIS paths across subset sizes. When tracing \(k=1, \dots, K\), 3 initializes the run for size \(k\) at the solution value from size \(k - 1\).

3.2.0.2 Computational Complexity

Each iteration of 2 has three steps. First, the scores \(s_i(\eta)=w_i+\eta c_i\) are computed in \(\mathcal{O}(n)\) time. Second, the top-\(k\) scores are selected, which can be done in expected \(\mathcal{O}(n)\) time using linear-time selection, or in \(\mathcal{O}(n \log k)\) time using a size-\(k\) heap. Once the set is selected, \(W({\mathbb{S}})\) and \(G({\mathbb{S}})\) are computed in \(\mathcal{O}(k)\) time. Thus, with linear-time selection, each iteration costs \(\mathcal{O}(n)\) time and \(\mathcal{O}(n)\) memory.

Let \(I_k\) denote the number of Dinkelbach iterations for subset size \(k\). The total cost for one size-\(k\) problem is therefore \(\mathcal{O}(I_k n)\). For an MIS path over \(k = 1, \dots, K\), the cost is \(\mathcal{O}\!\left( n \sum_{k=1}^K I_k \right),\) with warm starts from 3 typically reducing \(I_k\). Because the feasible class is finite, the algorithm terminates finitely; empirically, only a small number of iterations is needed.

3.3 Optimization Exactness of 2↩︎

The algorithm solves the fixed-input fractional optimization problem exactly. It converts the ratio problem into a sequence of linear maximization problems over a finite feasible class, each solved exactly by top-\(k\) selection.

Under absolute continuity and non-degeneracy of the realized weights and costs, ties between distinct feasible ratios occur with probability zero. In that case, \(\hat{\mathbb{S}}_k^{\max}\) is almost surely unique.

3.4 Statistical Validity for Partial Linear Models↩︎

2 solves the fixed-input fractional optimization problem exactly. In a partial linear model, however, the inputs \[w_i = \tilde{x}_i \tilde{r}_i, \qquad c_i = \tilde{x}_i^2\] are generated by first-stage residualization. The statistical question is whether the empirical problem is close to the oracle problem based on the unknown residuals \[v_i = x_i - h_0(Z_i), \qquad u_i = y_i - x_i \beta_0 - g_0(Z_i).\]

For fixed \(k\), the relevant first-order oracle objective is additive. Let \[\phi_i \mathrel{\vcenter{:}}= \frac{v_i u_i}{\mu_v} \quad\text{with} \quad \mu_v \mathrel{\vcenter{:}}= \mathbb{E}[v_i^2] > 0,\] and define \[Q_n^{\mathrm{or}}({\mathbb{S}}) \mathrel{\vcenter{:}}= \sum_{i\in{\mathbb{S}}} \phi_i, \qquad |{\mathbb{S}}|=k.\] Let \[{\mathbb{S}}_{k}^{\mathrm{or}} \in \arg\max_{|{\mathbb{S}}|=k} Q_n^{\mathrm{or}}({\mathbb{S}})\] denote an oracle first-order \(k\)-MIS. The empirical scaled objective is \[\widehat Q_n({\mathbb{S}}) \mathrel{\vcenter{:}}= n\{\hat{\beta}-\hat{\beta}_{-{\mathbb{S}}}\} = n \frac{ \sum_{i\in{\mathbb{S}}}\tilde{x}_i\tilde{r}_i }{ \sum_{j\notin{\mathbb{S}}}\tilde{x}_j^2 }.\]

For fixed \(k\), define \[E_{{\mathbb{S}}} \mathrel{\vcenter{:}}= \sum_{i\in{\mathbb{S}}}(\tilde{x}_i\tilde{r}_i-v_i u_i), \qquad A_{{\mathbb{S}}} \mathrel{\vcenter{:}}= \sum_{i\in{\mathbb{S}}} v_i u_i,\] and \[\begin{align} \delta_{n,k} \mathrel{\vcenter{:}}=& \sup_{{\mathbb{S}}\subset \left[ n \right], |{\mathbb{S}}|=k}|E_{{\mathbb{S}}}|,\\ B_{n,k} \mathrel{\vcenter{:}}=& \sup_{{\mathbb{S}}\subset \left[ n \right], |{\mathbb{S}}|=k}|A_{{\mathbb{S}}}|,\\ \rho_{n,k} \mathrel{\vcenter{:}}=& \sup_{{\mathbb{S}}\subset \left[ n \right], |{\mathbb{S}}|=k} \left|\frac{1}{n}\sum_{j\notin{\mathbb{S}}}\tilde{x}_j^2-\mu_v\right|.\\ \end{align}\]

Primitive sufficient conditions include cross-fitting, bounded moments, non-degenerate treatment residual variance, and first-stage nuisance rates strong enough that the generated score error over fixed-size subsets is \(o_p(1)\); details are given in the Appendix.

4 Experiments↩︎

We evaluate the method in simulations and applications. The simulations verify exact optimization against enumeration in small problems, and they measure runtime and iteration counts at scales where enumeration is impossible. The applications then use MIS traces as diagnostics in randomized experiments, linearized prediction tasks, and benchmark regression problems.

Across settings, we compute MIS in both signed directions, identifying sets that maximally increase and maximally decrease the target. We summarize results with influence traces, plotting the optimal influence value as a function of \(k\). We also track non-nestedness events, \(\hat{\mathbb{S}}_k \not\subset \hat{\mathbb{S}}_{k+1}\), and set overlap using the Jaccard index \(J({\mathbb{A}},{\mathbb{B}}) \mathrel{\vcenter{:}}= |{\mathbb{A}}\cap{\mathbb{B}}|/|{\mathbb{A}}\cup{\mathbb{B}}|.\)

4.1 Simulations↩︎

We use simulations for two main purposes: (i) to validate theoretical predictions and (ii) to characterize runtime and iteration counts at scale.

4.1.0.1 Optimization exactness

To validate 1, we compare 2 with enumeration and greedy baselines in problems where exhaustive search is feasible (\(n \leqslant 50\), \(k \leqslant 3\)). We use data-generating processes (DGPs) designed to stress the selection problem, including heteroskedasticity, heavy tails, mixtures, nonlinearities, endogeneity, autocorrelated errors, and masking. Across more than 10,000 Monte Carlo replications per design, 2 always matches the enumerated optimum, while greedy selection sometimes returns suboptimal sets.

4.1.1 Residualized Performance↩︎

We next assess the residualized setting motivated by 2. The simulation design contains two seeded influential groups of three observations each, built-in masking [15], and nonparametric confounding. We vary \(n \in \{1000, 5000, 10000\}\), estimate the nuisance functions with gradient boosting and \(5\)-fold cross-fitting, and compare the recovered empirical MIS with the known oracle-style score benchmark. Details are given in Appendix 10.

The residualized influence values track the oracle values closely, even at moderate sample sizes. Set recovery is more demanding; strongly influential sets are recovered reliably, whereas weakly separated or moderately influential sets converge more slowly. This matches 2 — value approximation is stable under weaker conditions than exact recovery of the maximizing set.

4.1.2 Runtime and Feasibility↩︎

We benchmark an implementation of 2 on synthetic univariate regressions across a grid of sample sizes \(n\) and set sizes \(k\), recording wall-clock time (across 100 runs) and the number of Dinkelbach iterations. As shown in 4, the median wall-clock time for \(n = 10^6, k = 10^5\) remains below \(200\) ms.³ Across the entire grid, 2 converged in a median of three iterations, with a maximum of six.

To stress-test convergence, we further scale the implementation and use adversarial warm starts. At \(n=10^8\), runtime is roughly five seconds for \(k=10^6\) and 60 seconds for \(k=10^7\). At \(n=10^9\), where memory becomes limiting, runtime is about ten seconds for \(k=10^6\), 80 seconds for \(k=10^7\), and fifteen minutes for \(k=10^8\). For the \(n=10^8,k=10^7\) case, replacing the warm start with \(\eta^{(0)}\in\{10^3,10^6,10^9,10^{12},10^{18}\}\) delays convergence by only one iteration, from five to six.

4.2 Applications↩︎

We use the method as a set-level sensitivity diagnostic in three classes of applications. In randomized experiments, the target is a fixed residualized univariate ratio and the algorithm identifies exact empirical MIS. In linearized prediction of semantic similarity, MIS quantify which training subsets most affect a local scalar prediction. In observational benchmark regressions, the PLM theory gives a first-order oracle interpretation when residualized-score stability is credible. In all applications we compute signed MIS in both directions.

4.2.1 Randomized Controlled Trials↩︎

We revisit the seven microcredit RCTs studied by [27] and [12]. These trials are useful stress tests because average treatment effect (ATE) estimates can be sensitive to small influential subsets — in two trials removing a single observation flips the sign of the estimated ATE, and sets of size 15 suffice across all trials. Previous work has used approximations, as enumeration is infeasible at the smallest size (\(n \geqslant 961\)). Our algorithm traces exact \(k\)-MIS over the full range of \(k\) at negligible computational cost.

5 shows the Mongolia trial of [26]. The increasing and decreasing traces differ substantially. In the increasing direction, non-nestedness is frequent but mostly local: the greedy trace remains close in value even as set overlap declines. In the decreasing direction, one non-nestedness event causes a large and persistent shift in the optimal set. Because greedy selection cannot revise earlier choices, it remains trapped on a suboptimal path and the influence gap accumulates.

6 reports the remaining six trials. Non-nestedness is common across trials, confirming that optimal MIS paths are generally not nested in \(k\). Greedy selection often approximates influence magnitudes well, but it need not recover the maximizing sets. The Philippines trial features the clearest asymmetry — the increasing trace is initially stable and later accelerates, while the decreasing trace follows the opposite pattern. These results suggest that value approximation and set recovery should be evaluated separately, consistent with 2.

4.2.2 Linearized Text Embeddings↩︎

Next, we use MIS for a linearized prediction task. We predict semantic similarity for \(n=5{,}749\) English sentence pairs from the STS benchmark [28]. Each pair is mapped to frozen sentence embeddings using a pre-trained transformer model [29], and the scalar regressor is the cosine similarity between the two embeddings. For three query pairs with low, medium, and high similarity, we study the local prediction target \(\phi(\beta) = x_0 \beta\).

1 reports the full-sample fitted score, the fitted scores after removing the \(5{,}000\)-MIS in each direction, and the smallest \(k\) for which the fitted score exceeds \(1.0\) for three specific sentences (with minimal, average, and maximal similarity). Because the target scales with \(x_0\), high-similarity query pairs are mechanically more sensitive to removal of influential training subsets. The traces reveal substantial sensitivity for the medium- and high-similarity pairs, including cases where deleting a subset pushes the linearized score outside the nominal \([0,1]\) range. Further details are provided in Appendix 10.

4.2.3 Additional Applications↩︎

We also apply the method to benchmark datasets from economics, biology, statistics, and machine learning; details are in Appendix 10. Across these datasets, greedy selection usually tracks the optimal influence value closely, but exact and greedy sets often differ. Most non-nestedness events are short-lived, in contrast to the persistent failure seen in the Mongolia trial. Warm-starting along the MIS path in 3 keeps iteration counts small, with at most four iterations in these applications. Runtime is dominated by fitting the original model rather than by MIS selection.

5 Discussion↩︎

This paper makes most influential set analysis computationally explicit for a broad class of fixed-input univariate objectives. The key reduction is algebraic — once leave-set-out influence is written as a score-over-curvature ratio, the search over \(\binom{n}{k}\) subsets reduces to repeated top-\(k\) selection along a one-dimensional parameter path. For the fixed residualized inputs covered by 1, the result is the exact global optimum, not a greedy or local approximation.

The statistical interpretation depends on how those inputs are generated. In randomized experiments and other fixed-input settings, the empirical MIS is the target. In partial linear models, 2 shows that the scaled residualized objective has a first-order oracle interpretation when residualized scores and denominators are stable. This distinction is important: orthogonality supports value stability, but exact recovery of the selected set also requires separation between the best and second-best oracle sets.

5.1 Practical Use↩︎

We recommend reporting signed MIS traces over a range of \(k\). The increasing and decreasing traces answer different robustness questions and can implicate different subsets. Large early jumps indicate sensitivity to few observations; gradual traces suggest dispersed influence; sharp composition changes signal non-nestedness and potential failure of greedy approximations. The implicated sets should then be inspected directly, especially when they contain duplicated records, unusual covariate patterns, high-leverage observations, or observations tied to a common source.

5.1.0.1 Robustness Auditing

MIS traces provide target-specific worst-case deletion diagnostics. Compared with singleton influence measures, they reveal whether sensitivity is concentrated in one record or distributed across interacting observations. Compared with ad-hoc deletion checks, each set is the most adverse set of its size for the chosen target. Compared with greedy approximations, the exact trace is not constrained to follow a nested deletion path. Our experiments show that greedy methods may approximate the influence value while missing the maximizing set, and in some cases a single non-nestedness event produces persistent error.

5.1.0.2 Data Curation

MIS can guide preprocessing and data cleaning without relying on blunt transformations such as trimming or winsorization. The selected sets identify which observations actually drive the downstream estimate or prediction. This can prioritize manual review, reveal duplicated or contaminated records, and document whether proposed cleaning rules address the relevant sources of sensitivity.

5.1.0.3 Prediction Diagnostics

For linearized prediction targets, MIS provide dataset-level explanations that complement feature-level attributions. A large singleton effect indicates that one example matters on its own. A large set effect indicates that the prediction is anchored by a subset whose collective role may not be visible from individual records. This is useful for identifying duplicated, highly similar, or mutually masking examples that jointly support a prediction.

5.2 Limitations↩︎

The computational reduction requires a linear-fractional leave-set-out objective. This covers residualized univariate least-squares targets and related affine functionals, but non-linear estimands, vector-valued targets, and general \(M\)-estimators may require approximation, local linearization, or different optimization subproblems.

The PLM result is also deliberately first-order. For fixed \(k\), the scaled residualized objective is close to an oracle orthogonal-score objective under residualized-score stability. This does not imply that arbitrary first-stage learners produce identical MIS, nor that exact set recovery is guaranteed without a gap. When first-stage learners are unstable for high-leverage or extreme-score observations, the recovered set remains the exact empirical MIS for the residualized regression but may not estimate the oracle first-order MIS.

Finally, exact MIS are diagnostic rather than prescriptive. A highly influential set may reveal contamination, model misspecification, subgroup heterogeneity, or genuine scientific signal. The method identifies the subset responsible for sensitivity; substantive interpretation still requires domain knowledge.

5.3 Extensions↩︎

Several extensions are natural. Ridge-stabilized denominators preserve the same basic reduction and can improve behavior when residual curvature is small. For broader classes of smooth estimators, approximate score-over-curvature representations may yield local or asymptotic analogues of the present algorithm. Structured versions of MIS are also promising; replacing the top-\(k\) step with a constrained selection subproblem would allow influential clusters, time-contiguous blocks, group-balanced removals, or stratified deletion rules.

A second direction is inference on the traces themselves. Exact MIS computation makes it feasible to separate ordinary sampling variability from unusually influential subsets, and can provide candidate sets for formal tests of excessive influence [19]. This shifts the next question from how to find influential sets to which influential sets are substantively or statistically anomalous.

5.4 Broader Implications↩︎

Efficient MIS computation turns set-level sensitivity analysis from a bespoke forensic exercise into a routine diagnostic. The method complements robust estimation, which limits the effect of contamination; meanwhile MIS identify which observations are responsible for sensitivity and how that sensitivity accumulates with set size. In high-stakes empirical and machine-learning pipelines, such inspectable counterfactual subsets can improve transparency, support data provenance, and clarify the relationship between data and conclusions.

6 Conclusion↩︎

We introduced an exact algorithm for identifying most influential sets when leave-set-out effects have a linear-fractional representation. The method reduces the combinatorial search over \(\binom{n}{k}\) subsets to a finite sequence of top-\(k\) selection problems, with \(\mathcal{O}(n)\) cost per iteration and few iterations in practice.

The guarantees separate computation from statistical interpretation. For fixed residualized inputs, the algorithm returns a globally optimal empirical MIS. For partial linear models, the scaled residualized influence objective is first-order equivalent to an oracle orthogonal-score objective under residualized-score stability; value consistency follows directly, and exact set recovery follows under a separation condition.

Empirically, the method scales to datasets far beyond the reach of enumeration. Across randomized experiments, semantic similarity prediction, and benchmark regressions, exact MIS traces reveal frequent non-nestedness and show where greedy approximations succeed or fail. These traces provide actionable diagnostics for robustness auditing, data curation, and prediction analysis. Future work can extend the framework to structured deletion constraints, vector-valued or nonlinear targets, and formal inference for anomalous influence.

Impact Statement↩︎

This paper develops methods for identifying most influential sets, whose removal induces the greatest change in a target quantity of interest. The primary goal is to improve robustness auditing and clarify how empirical conclusions depend on specific subsets of data.

Potential positive impacts include improved accountability, data provenance, and reliability. By identifying the subsets that drive estimates or predictions, the method can help practitioners detect contamination, duplicated records, unstable estimates, subgroup dependence, or other sources of sensitivity. It can also guide targeted data cleaning and make robustness analyses more reproducible.

A potential risk is adversarial use. An actor could use influence information to identify records whose inclusion or removal would steer a model or estimate in a desired direction. We view transparency as the appropriate mitigation. Set-level sensitivity already exists whether or not it is measured; efficient diagnostics make these dependencies visible, auditable, and easier to report.

The computational footprint is modest relative to enumeration. In our benchmarks, the method scales to millions of observations on a single thread, and the runtime is typically dominated by fitting the original model rather than by MIS selection.

7 Proof of 1↩︎

7.0.0.1 Problem Formulation

Let \(n \in \mathbb{N}\) and \(k \in \{1, \ldots, n\}\). Given weights \(w_i, c_i \in \mathbb{R}\) for \(i = 1, \ldots, n\) and \(T \in \mathbb{R}\) such that \[T - \sum_{i \in {\mathbb{S}}} c_i > 0 \quad \text{for all } {\mathbb{S}}\subset [n] \text{ with } |{\mathbb{S}}| = k,\] we consider the fractional programming problem: \[\label{eq:main95problem} \max_{{\mathbb{S}}\subset [n], \, |{\mathbb{S}}| = k} \frac{W({\mathbb{S}})}{G({\mathbb{S}})},\tag{5}\] where \(W({\mathbb{S}}) \mathrel{\vcenter{:}}= \sum_{i \in {\mathbb{S}}} w_i\) and \(G({\mathbb{S}}) \mathrel{\vcenter{:}}= T - \sum_{i \in {\mathbb{S}}} c_i\).

Dinkelbach’s Method↩︎

For \(\eta \in \mathbb{R}\), define the parametric objective function: \[\label{eq:parametric} F_\eta({\mathbb{S}}) \mathrel{\vcenter{:}}= W({\mathbb{S}}) - \eta \, G({\mathbb{S}}) = \sum_{i \in {\mathbb{S}}}(w_i + \eta c_i) - \eta T.\tag{6}\] For each \(\eta \in \mathbb{R}\), define \[z(\eta) \mathrel{\vcenter{:}}= \max_{{\mathbb{S}}\subset[n],\,|{\mathbb{S}}|=k} \{W({\mathbb{S}})-\eta G({\mathbb{S}})\}.\] Because \(G({\mathbb{S}}) > 0\) on the feasible class, \(\eta^\star = \max_{|{\mathbb{S}}|=k}\frac{W({\mathbb{S}})}{G({\mathbb{S}})}\) if and only if \(z(\eta^\star)=0\). Moreover, \(z(\eta)>0\) for \(\eta<\eta^\star\) and \(z(\eta)<0\) for \(\eta>\eta^\star\).

The key observation [22], [30] is that solving the fractional program equation 5 is equivalent to finding the parameter \(\eta^*\) such that \(z(\eta^*) = 0\).

7.0.0.2 Proof

Main Argument: By Dinkelbach’s fundamental result, the algorithm terminates when \(z(\eta^*) = 0\), which certifies that \(\eta^* = \Delta_k^{\max}\) and the corresponding selection is globally optimal. Lemmas 1–3 guarantee that this occurs in finite iterations. Thus, the algorithm outputs an exact global optimum. 0◻

8 Proof of 2↩︎

Problem Formulation↩︎

Define \[D_n({\mathbb{S}}) \mathrel{\vcenter{:}}= \sum_{j\notin{\mathbb{S}}}\tilde{x}_j^2, \qquad d_n({\mathbb{S}}) \mathrel{\vcenter{:}}= \frac{D_n({\mathbb{S}})}{n} - \mu_v,\] \[A_{{\mathbb{S}}} \mathrel{\vcenter{:}}= \sum_{i\in{\mathbb{S}}} v_i u_i, \qquad E_{{\mathbb{S}}} \mathrel{\vcenter{:}}= \sum_{i\in{\mathbb{S}}}\bigl(\tilde{x}_i \tilde{r}_i - v_i u_i\bigr),\] so that \(\sum_{i\in{\mathbb{S}}}\tilde{x}_i \tilde{r}_i = A_{{\mathbb{S}}}+E_{{\mathbb{S}}}\), and \[\widehat Q_n({\mathbb{S}}) = n\bigl(\hat{\beta} - \hat{\beta}_{-{\mathbb{S}}}\bigr) = \frac{A_{{\mathbb{S}}}+E_{{\mathbb{S}}}}{\mu_v + d_n({\mathbb{S}})},\qquad Q_n^{\mathrm{or}}({\mathbb{S}}) = \frac{A_{{\mathbb{S}}}}{\mu_v}.\] Write \[\begin{align} \delta_{n,k} \mathrel{\vcenter{:}}= \sup_{|{\mathbb{S}}|=k}|E_{{\mathbb{S}}}|,\; B_{n,k} \mathrel{\vcenter{:}}= \sup_{|{\mathbb{S}}|=k}|A_{{\mathbb{S}}}|,\; \rho_{n,k} \mathrel{\vcenter{:}}= \sup_{|{\mathbb{S}}|=k}|d_n({\mathbb{S}})|, \end{align}\] all controlled by Assumption 2; in particular \(\rho_{n,k}=o_p(1)\).

Proof↩︎

8.0.0.1 Step 1: uniform approximation.

A direct rearrangement gives the identity \[\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}}) \;=\; \frac{E_{{\mathbb{S}}}}{\mu_v+d_n({\mathbb{S}})} \;-\; \frac{A_{{\mathbb{S}}}\,d_n({\mathbb{S}})}{\mu_v\bigl(\mu_v+d_n({\mathbb{S}})\bigr)},\] hence \[\bigl|\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}})\bigr| \;\le\; \frac{|E_{{\mathbb{S}}}|}{|\mu_v+d_n({\mathbb{S}})|} \;+\; \frac{|A_{{\mathbb{S}}}|\,|d_n({\mathbb{S}})|}{\mu_v\,|\mu_v+d_n({\mathbb{S}})|}.\] Since \(\rho_{n,k}=o_p(1)\), the event \(\bigl\{\inf_{|{\mathbb{S}}|=k}\,(\mu_v+d_n({\mathbb{S}}))\ge \mu_v/2\bigr\}\) has probability tending to one. On this event, taking suprema yields \[\sup_{|{\mathbb{S}}|=k} \bigl|\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}})\bigr| \;\le\; \frac{2\,\delta_{n,k}}{\mu_v} \;+\; \frac{2\,B_{n,k}\,\rho_{n,k}}{\mu_v^{2}} \;=\; o_p(1).\]

8.0.0.2 Step 2: value consistency.

Let \[\hat{\mathbb{S}}_k\in\arg\max_{|{\mathbb{S}}|=k}\widehat Q_n({\mathbb{S}}), \qquad {\mathbb{S}}_k^{\mathrm{or}}\in\arg\max_{|{\mathbb{S}}|=k}Q_n^{\mathrm{or}}({\mathbb{S}}).\] By optimality of \(\hat{\mathbb{S}}_k\), \(\widehat Q_n(\hat{\mathbb{S}}_k)\ge \widehat Q_n({\mathbb{S}}_k^{\mathrm{or}})\), and combining this with the uniform approximation gives \[Q_n^{\mathrm{or}}(\hat{\mathbb{S}}_k) \;\ge\; Q_n^{\mathrm{or}}({\mathbb{S}}_k^{\mathrm{or}}) \;-\; 2\sup_{|{\mathbb{S}}|=k}\bigl|\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}})\bigr|.\] Hence the oracle value of the empirical maximizer is within \(o_p(1)\) of the oracle optimum.

8.0.0.3 Step 3: selection consistency.

Assume the oracle maximizer is unique, with gap \[\Gamma_{n,k} \;\mathrel{\vcenter{:}}=\; Q_n^{\mathrm{or}}({\mathbb{S}}_k^{\mathrm{or}}) \;-\; \max_{\substack{|{\mathbb{S}}|=k\\ {\mathbb{S}}\ne{\mathbb{S}}_k^{\mathrm{or}}}} Q_n^{\mathrm{or}}({\mathbb{S}}).\] If \(\hat{\mathbb{S}}_k\ne{\mathbb{S}}_k^{\mathrm{or}}\), then \(Q_n^{\mathrm{or}}(\hat{\mathbb{S}}_k)\le Q_n^{\mathrm{or}}({\mathbb{S}}_k^{\mathrm{or}})-\Gamma_{n,k}\), and Step 2 forces \[2\sup_{|{\mathbb{S}}|=k}\bigl|\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}})\bigr| \;\ge\; \Gamma_{n,k}.\] Therefore, provided \[\Gamma_{n,k}^{-1} \sup_{|{\mathbb{S}}|=k}\bigl|\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}})\bigr| \;=\; o_p(1), \label{eq:ratio95condition}\tag{7}\] we have \(\Pr\!\bigl(\hat{\mathbb{S}}_k={\mathbb{S}}_k^{\mathrm{or}}\bigr)\to 1\).

8.1 When does condition equation 7 hold?↩︎

Since \(\rho_{n,k} = o_p(1)\) gives \(\inf_{|{\mathbb{S}}|=k}\bigl\{\mu_v+d_n({\mathbb{S}})\bigr\}=\mu_v+o_p(1)\), the bound of Step 1 sharpens to the rate \[\sup_{|{\mathbb{S}}|=k}\bigl|\widehat Q_n({\mathbb{S}}) - Q_n^{\mathrm{or}}({\mathbb{S}})\bigr| \;=\; O_p\!\left( \frac{\delta_{n,k}}{\mu_v} \;+\; \frac{B_{n,k}\,\rho_{n,k}}{\mu_v^{2}} \right).\] Both summands being nonnegative, condition equation 7 is equivalent to \[\Gamma_{n,k}^{-1}\bigl(\delta_{n,k} + B_{n,k}\,\rho_{n,k}\bigr) \;=\; o_p(1),\] which in turn is equivalent to the two margin-dominance conditions holding jointly: \[\boxed{ \delta_{n,k} \;=\; o_p(\Gamma_{n,k}), \qquad B_{n,k}\,\rho_{n,k} \;=\; o_p(\Gamma_{n,k}). }\]

9 Joint Influence, Masking, and Greedy Failure↩︎

The following definitions are stated for a generic set objective \[H_k({\mathbb{S}}), \qquad {\mathbb{S}}\subset[n], \quad |{\mathbb{S}}|=k,\] where, in the main applications, \[H({\mathbb{S}}) = \frac{W({\mathbb{S}})}{G({\mathbb{S}})} = \frac{ \sum_{i\in{\mathbb{S}}} w_i }{ T - \sum_{i\in{\mathbb{S}}}c_i }.\] For signed analyses in the opposite direction, the same definitions apply after replacing \(H\) by the corresponding direction-specific objective.

For the linear-fractional objective in the main text, marginal gains depend on the current set. Writing \(w_i=\tilde{x}_i\tilde{r}_i\), \(c_i=\tilde{x}_i^2\), \(W({\mathbb{S}})=\sum_{j\in{\mathbb{S}}}w_j\), and \(G({\mathbb{S}})=T-\sum_{j\in{\mathbb{S}}}c_j\), we have \[H({\mathbb{S}}\cup\{i\})-H({\mathbb{S}}) = \frac{ w_iG({\mathbb{S}})+c_iW({\mathbb{S}}) }{ G({\mathbb{S}})\{G({\mathbb{S}})-c_i\} }.\] Thus the gain from adding \(i\) depends not only on its own score \(w_i\) and curvature \(c_i\), but also on the accumulated score \(W({\mathbb{S}})\) and remaining curvature \(G({\mathbb{S}})\) of the current set. This set dependence allows joint influence, masking, non-nested optimal paths, and greedy failure.

10 Additional Results↩︎

Influence Convergence Simulation↩︎

We consider the partial linear model from 2. To evaluate the finite sample behavior of 2, we generate data with the following structure. Covariates are drawn as \(Z_{ij} \stackrel{\text{iid}}{\sim} \text{Uniform}(0,1)\) for \(j = 1, \ldots, p_z\). Treatment assignment incorporates confounding through \[X_i = h(\mathbf{Z}_i) + v_i,\] where \(h(Z) = 0.5 z_1 - 0.2 z_2\) and \(v_i\) follows a mixture distribution to create high leverage sets: \[v_i \sim \begin{cases} \mathcal{N}(8, 0.25^2) & \text{for } i = n, n-1, n-2, \\ \mathcal{N}(6, 0.25^2) & \text{for } i = n-3, n-4, n-5,\\ \mathcal{N}(0, 1) & \text{for } 1, \ldots, n-6. \end{cases}\]

The nonparametric component in \(Y\) exhibits sparsity, depending on only four covariates: \[g(\mathbf{z}) = 2\sin(2\pi z_1) + 1.5 z_2^2 + z_3 z_4.\]

10.0.0.1 Heterogeneous Treatment Effects.

To ensure high influence for some observations, we generate outcomes with treatment effect heterogeneity, violating constant \(\beta\): \[Y_i = g(\mathbf{Z}_i) + \beta_i X_i + \varepsilon_i,\] where \[(\beta_i, \varepsilon_i) \sim \begin{cases} (-0.5, \mathcal{N}(0, 1)) & \text{for } i = 1, \ldots, n-6 \\ (0.1, \mathcal{N}(0, 0.1^2)) & \text{for } i = n-5, n-4, n-3 \\ (0.4, \mathcal{N}(0, 0.1^2)) & \text{for } i = n-2, n-1, n. \end{cases}\] This design extends the illustration of masking and joint influence in [15]. Each setup is resampled, for a total of one thousand Monte Carlo draws. 7 displays the mean and \(95\%\) confidence bounds of the simulation results.

10.0.0.2 Simulation Results

While influence converges quickly and uniformly to the oracle, the associated \({\mathbb{S}}^\star\) converges much slower. While even for medium sized sample size residualization works well with an average Jaccard similarity of almost one for all sample sizes, accuracy is substantially worse if added observations are not influential. Thus knowing the correct \(k\) is beneficial for set accuracy in finite samples. If inference on \({\mathbb{S}}^\star\) is not desired but merely robustness, even for moderate \(n\) residualization recovers the true \(\Delta_k^{\max}\) well.

10.1 Linearized Text Embeddings↩︎

The sentence-pairs and similarity scores are:

• ‘Romney picks Ryan as vice presidential running mate: source’

• ‘Romney to tap Ryan as vice presidential running mate’

• ‘An older man is standing outside in front of a truck.’

• ‘A woman dressed in green is roller skating outside at an event.’

• ‘A new study, conducted in Europe, found the medicine worked just as well as an earlier disputed study, sponsored by ImClone Systems, said it did.’

• ‘Doctors concluded Erbitux, the cancer drug that enmeshed ImClone Systems in an insider trading scandal, worked just as well as an earlier company-sponsored study said it did.’

Traces of the predicted similarity scores as \(k\) increases are shown in 8.

10.2 Additional Applications↩︎

The datasets considered are listed in 2; results are provided as Supplementary Material.

10.3 Comparison with Influence Function Approximations↩︎

9 compares exact \(\Delta(\mathbb{S})\) with their AMIP approximation [12].

References↩︎

---

1. The Appendix contains proofs and additional details; code is provided at https://github.com/nk027/findingMIS.↩︎

2. See the discussion in [13]–[15] for further discussion.↩︎

3. The implementation uses a size-\(k\) heap for ordering, and is implemented in R and Rcpp [25], on a Ryzen AI Max+ Pro 395, using a single thread. We also benchmarked an implementation of the greedy algorithm by [15] on a subset of the grid (a runtime of \(200\) ms is achieved for \(n = 10^4\) and \(k = 50\)), as well as the approximations implemented by [12]; across a range of settings for \((n, k)\), our method is typically faster to compute, often by an order of magnitude.↩︎