# A Survey of Data Attribution: Methods, Applications, and Evaluation in the Era of Generative AI

Junwei Deng, Yuzheng Hu, Pingbang Hu, Ting-Wei Li, Shixuan Liu, Jiachen

T. Wang, Dan Ley, Qirun Dai, Benhao Huang, Jin Huang, et al.

► To cite this version:

Junwei Deng, Yuzheng Hu, Pingbang Hu, Ting-Wei Li, Shixuan Liu, et al. A Survey of Data Attribution: Methods, Applications, and Evaluation in the Era of Generative AI. 2025. hal-05230469

HAL Id: hal-05230469

https://hal.science/hal-05230469v1

Preprint submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A SURVEY OF DATA ATTRIBUTION: METHODS, APPLICATIONS,
AND EVALUATION IN THE ERA OF GENERATIVE AI

WORKING PAPER

Junwei Deng¹* Yuzheng Hu¹*† Pingbang Hu¹* Ting-Wei Li¹* Shixuan Liu¹* Jiachen T. Wang²

Dan Ley³ Qirun Dai⁴ Benhao Huang⁵ Jin Huang⁶ Cathy Jiao⁵ Hoang Anh Just⁷ Yijun Pan⁶

Jingyan Shen⁸ Yiwen Tu⁹ Weiyi Wang⁶ Xinhe Wang⁵ Shichang Zhang³ Shiyuan Zhang¹

Ruoxi Jia⁷ Himabindu Lakkaraju³ Hao Peng¹ Weijing Tang⁵ Chenyan Xiong⁵ Jieyu Zhao¹⁰

Hanghang Tong¹ Han Zhao¹ Jiaqi W. Ma¹†

¹University of Illinois Urbana-Champaign ²Princeton University ³Harvard University
⁴University of Chicago ⁵Carnegie Mellon University ⁶University of Michigan ⁷Virginia Tech
⁸New York University ⁹University of California San Diego ¹⁰University of Southern California

August 29, 2025

ABSTRACT

Training data is the fuel of modern artificial intelligence (AI), fundamentally shaping the capabilities, limitations, and biases of AI systems. The emergence of large-scale generative models has elevated the importance of understanding how data influences their behaviors, bringing the field of *data attribution* to the forefront. This survey provides a comprehensive overview of data attribution, covering its methods, applications, and evaluation protocols, with a particular emphasis on the challenges and opportunities arising in the era of generative AI. We start by introducing a conceptual framework for attribution centered on three core questions: *what to attribute* (model behaviors), *attribute to what* (training entities), and *how to attribute* (influence measures). Within this framework, we systematically review major attribution approaches, including those based on influence functions, weighted marginal contributions, training dynamics, and simulators. We then examine key applications of data attribution, such as data selection, fact tracing, adversarial attacks and defenses, and the emerging data economy. Finally, we critically assess common evaluation criteria, including the quality of counterfactual predictions, utility in downstream tasks, and computational efficiency. We conclude with a forward-looking perspective on the future of data attribution, highlighting key open challenges and promising directions for future research.

**Keywords** Data Attribution · Generative AI

# 1 Introduction

Training data is the cornerstone of modern artificial intelligence (AI), significantly shaping model behaviors through its scale, quality, and composition. These data characteristics profoundly influence critical aspects of AI models, including performance (Kaplan et al., 2020), robustness (Gowal et al., 2021), privacy (Carlini et al., 2021), and safety (Qi et al., 2024). Consequently, precisely understanding how variations in training data affect model outcomes has become increasingly crucial. This has led to the emergence of *data attribution*, a field focused on systematically quantifying the contribution of individual training data samples or data subsets a model’s predictions, behaviors, or internal states. Since the seminal work by Koh and Liang (2017), which introduced influence functions (Hampel, 1974; Hampel et al., 2011) as a practical tool for data attribution, there has been a surge of new developments in this field, offering more scalable and precise techniques (Ghorbani and Zou, 2019; Pruthi et al., 2020; Hara et al., 2019; Park et al., 2023; Wang and Jia, 2023a; Wang et al., 2024b).

<footnote>*Co-first authors.
†Corresponding to {yh46, jiaqima}@illinois.edu.</footnote>
The advent of large-scale generative AI models, such as diffusion models (Ho et al., 2020b) and large language models (LLMs) (Brown et al., 2020), has placed an even greater emphasis on the importance of training data, in turn posing new demands on data attribution. These demands arise from several interconnected challenges and opportunities unique to the generative paradigm. First, the immense scale of these models, often comprising billions of parameters trained on vast datasets, renders many traditional attribution methods computationally infeasible and necessitates the development of highly efficient and scalable techniques (Schioppa et al., 2022; Kwon et al., 2024; Choe et al., 2024). Second, the complexity of generative tasks, from autoregressive language modeling to diffusion-based image generation, expands the scope of attribution itself. This requires new definitions for both the model behaviors to attribute (e.g., generation quality, factual recall (Akyürek et al., 2022), or safety (Cohen and Giryes, 2024)) and the training entities to which we attribute them (e.g., individual tokens (Yeh et al., 2022), data domains (Xie et al., 2023), or abstract concepts (Cohen-Wang et al., 2024) beyond single examples). Third, generative AI creates novel and critical applications for data attribution. These include verifying the factual grounding of generated text (Akyürek et al., 2022), ensuring compliance with copyright and privacy regulations by tracing generated content (Deng et al., 2024b), understanding emergent phenomena like in-context learning (Cohen-Wang et al., 2024), and diagnosing complex failure modes such as hallucination (Lin et al., 2024b). Finally, evaluating attribution methods in this new landscape is itself a major challenge. Traditional evaluation strategies that depend on repeated model retraining are hindered by prohibitive computational costs, while the inherent stochasticity of model training further complicates the reliability of assessments (Nguyen et al., 2023).

This survey provides a comprehensive overview of the field of data attribution, examining its **methods**, **applications**, and **evaluation**, with a special focus on its evolution in the context of generative AI. For methods (Section 2), we first introduce a conceptual framework for understanding data attribution, organized around three fundamental questions: “What to attribute?” (*model behaviors*), “Attribute to what?” (*training entities*), and “How to attribute?” (*influence measures*). This framework provides a foundation for understanding data attribution methods. Building on this framework, we systematically review state-of-the-art attribution methods, which we categorize according to the principles of their influence measures (i.e., “How to attribute?”). For applications (Section 3), we explore the diverse and emerging application scenarios of data attribution in generative AI, including *data selection*, *fact tracing*, *adversarial attacks and defenses*, and *data economy*. For evaluation (Section 4), we critically assess the major types of evaluation criteria from the literature, including *counterfactual prediction of model behavior*, *utility in downstream tasks*, and *computational costs*, and highlight their strengths and limitations. Finally, we conclude the survey with a discussion of the survey scope (Section 5), and a forward-looking perspective on the future of data attribution in the era of generative AI (Section 6). Figure 1 provides a visual roadmap of these components and their interconnections.

<figure><img src="image_2.png" /><figcaption>Figure 1: Overview of the key components in this survey, comprising methods, applications, and evaluation. The methods section presents a high-level conceptual framework consisting of training entities, influence measures, and model behaviors. The applications section explores various sub-applications that can leverage data attribution. The evaluation section categorizes existing evaluation methods into distinct types. Corresponding section numbers are indicated in parentheses.</figcaption></figure>
## 2 Methods

This section provides a systematic study of training data attribution methods. We begin by introducing the fundamental concepts of data attribution. We then review mainstream attribution methods, discussing their definitions, underlying principles, and challenges in scaling to large generative models. Finally, we compare these approaches in terms of efficiency, attribution quality, and applicability, highlighting their trade-offs and practical considerations.

**Definition.** *Data attribution* measures the *influence* of each *training entity* on the *model behavior*.

Section 2.1.3

Section 2.1.2

Section 2.1.1

### 2.1 Foundations of Data Attribution

<figure><img src="image_3.png" /><figcaption>Figure 2: An overview of the three fundamental components in data attribution. Models are trained on *training entities*, which leads to various types of *model behaviors*. Data attribution quantifies the *influence measure* of each training entity on these behaviors.</figcaption></figure>

We decompose data attribution into three fundamental components: **model behavior**, which represents the outcome to be attributed; **training entity**, the source responsible for these behaviors; and **influence measure**, which quantifies the contribution of training entities to the observed model behavior. See Figure 2 for an overview. We systematically review each component, providing the necessary foundation for understanding different training data attribution methods.

#### 2.1.1 Model Behavior: What to Attribute?

The definition of model behavior centers around the question: **what to attribute?** This involves identifying the *measurable* effects of a model, typically tied to a well-defined *target function*. These effects can be *static*, such as a trained model’s prediction or loss evaluated on a specific test sample or an entire test set, or *dynamic*, capturing how these predictions or losses evolve throughout the training process.

Traditionally, data attribution has been primarily applied to standard supervised learning tasks, such as classification (Koh and Liang, 2017; Yeh et al., 2018; Pruthi et al., 2020). However, the emergence of advanced generative models, including diffusion models (Ho et al., 2020b; Nichol and Dhariwal, 2021; Song and Ermon, 2019; Song et al., 2021b,c) and large language models (LLMs) (Brown et al., 2020; Ouyang et al., 2022; Achiam et al., 2023; Touvron et al., 2023a,b; Dubey et al., 2024), has significantly expanded the scope of model behavior. These models, particularly due to their strong multi-task capabilities (Radford et al., 2019), introduce richer and more diverse definitions of model behaviors, giving rise to novel attribution challenges (e.g., fact tracing (Akyürek et al., 2022; Chang et al., 2025), also see Section 3 for more examples). To provide a structured overview, Table 1 summarizes commonly used target functions across different generative models, along with their representative references.

#### 2.1.2 Training Entity: Attribute to What?

In parallel, the definition of training entity addresses the question: **attribute to what?** At its core, this focuses on identifying the *source* that drives the model's behavior, providing insights into its underlying mechanisms. The vast majority of the literature attributes model behavior to *individual training samples*, which is the most natural and
<table><caption>Table 1: Commonly used target functions for different types of generative models.</caption><thead><tr><th>Model Type</th><th>Target Functions</th><th>Representative References</th></tr></thead><tbody><tr><td>Decoder-Only LLMs</td><td>Negative log likelihood<br>Multi-class Margin</td><td>Choe et al. (2024)<br>Park et al. (2023)</td></tr><tr><td>Diffusion Models</td><td>DDPM loss<br>Norm-based noise prediction<br>Noise prediction</td><td>Georgiev et al. (2023)<br>Zheng et al. (2023)<br>Lin et al. (2025)</td></tr><tr><td>Other Models</td><td>CLIP loss</td><td>Park et al. (2023)</td></tr></tbody></table>

widely applicable formulation. Accordingly, our discussion in Sections 2.2 to 2.6 will primarily consider attribution to individual samples unless otherwise noted.

While individual sample attribution is the standard, some studies explore alternative formulations. In LLMs, attribution can be performed at a finer granularity, such as *token-level* attribution (Lin et al., 2024b), where the influence of individual tokens is assessed. Another approach is *concept-level* attribution (Cohen-Wang et al., 2024), which aims to identify broader semantic structures such as key phrases, document sections, or contextual passages that are influential. Conversely, attribution can be extended to *groups of samples* (Koh et al., 2019), often under the assumption of *additivity* (Hu et al., 2024a). Additionally, in settings involving multiple data distributions, some methods attribute model behavior to entire datasets or training domains (Wang et al., 2023a; Liu et al., 2025). We summarize different training entity formulations in the literature, along with their representative references, in Table 2.

<table><caption>Table 2: Granularity of training entities.</caption><thead><tr><th>Granularity</th><th>Description</th><th>Representative References</th></tr></thead><tbody><tr><td>Token</td><td>Attribution to individual tokens</td><td>Yeh et al. (2022); Lin et al. (2024b)</td></tr><tr><td>Concept</td><td>Attribution to broader semantic structures</td><td>Cohen-Wang et al. (2024); Chuang et al. (2025)</td></tr><tr><td>Sample</td><td>Attribution to individual samples</td><td>Grosse et al. (2023); Park et al. (2023), Wang et al. (2024g,h, 2025b)</td></tr><tr><td>Group</td><td>Attribution to groups of training samples</td><td>Koh et al. (2019); Basu et al. (2020), Hu et al. (2024a); Ley et al. (2024)</td></tr><tr><td>Domain</td><td>Attribution to training domains</td><td>Xie et al. (2023); Wang et al. (2023a), Kang et al. (2024); Liu et al. (2025)</td></tr></tbody></table>

### 2.1.3 Influence Measure: How to Attribute?

An influence measure serves as the *channel* between training entities and model behavior by providing a systematic method to quantify the impact of each training entity on the observed behavior. It answers the question of **how to attribute**, establishing a computational framework for linking sources to outcomes.

Influence measures typically assign *scores* to training entities, most commonly individual training samples. These scores reflect the importance of each training entity in shaping the model's behavior and are typically interpreted in two ways: 1) **Magnitude**: A larger absolute value indicates a stronger influence; 2) **Sign**: A positive score suggests a training entity reinforces the model's behavior, whereas a negative score indicates an opposing effect.

Within the extensive literature on data attribution, we identify the following primary categories of influence measures:

* **Influence function** (Section 2.2): Use first-order approximations of model parameters to efficiently approximate the effect of removing a single training sample (i.e., leave-one-out (LOO)) without requiring full retraining. The seminal work in this line is Koh and Liang (2017). Additionally TRAK (Park et al., 2023) is introduced as a variant of influence function.

* **Weighted marginal contribution** (Section 2.3): Quantify sample influence based on its expected marginal contribution when added to subsets of the training data, where each subset is assigned a specific weight. Methods such as Data Shapley (Ghorbani and Zou, 2019; Jia et al., 2019a) and Data Banzhaf (Wang and Jia, 2023a) fall into this category.
* **Training dynamic methods** (Section 2.4): Leverage intermediate model checkpoints to analyze how training samples dynamically affect model behaviors along the training trajectory. Two prominent methods in this category are TracIn (Pruthi et al., 2020) and SGD-influence (Hara et al., 2019).

* **Simulator** (Section 2.5): Use surrogate models (e.g., Datamodels (Ilyas et al., 2022)) to estimate the counterfactual effect of model training. We will also discuss the connections between this category and weighted marginal contribution methods.

* **Others** (Section 2.6): Methods that do not fit into the above categories, such as those based on unlearning (Wang et al., 2024h) or similarity (Yang et al., 2025).

Influence measures form the core of data attribution, providing a quantitative framework for assessing the impact of training entities and are typically agnostic to the choice of target functions. In Sections 2.2 to 2.6, we analyze specific methods in detail. For each category, we first review early approaches that established the core methodologies. We then examine the challenges of scaling these methods to large-scale generative models. Finally, we discuss recent advancements aimed at addressing these challenges and highlight open questions that remain in the field.

## 2.2 Category: Influence Function

The first category of methods focuses on the influence function (IF), a tool originating in robust statistics (Hampel, 1974; Hampel et al., 2011) and is closely related to cross validation (Debruyne et al., 2008) and jackknife (Giordano et al., 2019b). IF was first introduced to deep learning in the seminal work by Koh and Liang (2017). It provides a principled way to understand how changes in individual training data points affect model parameters and predictions.

**Notation.** Let $z = (x, y)$ denote a data point, where $x$ is the input and $y$ is the corresponding label. The loss incurred on a single data point is denoted by $\ell(z; \theta)$, where $\theta \in \mathbb{R}^p$ represents the model parameters. Given a training set $D = \{z_i\}_{i=1}^n$, the empirical risk is defined as $\mathcal{L}(\theta) = \frac{1}{n} \sum_{i=1}^n \ell(z_i; \theta)$, and the optimal model parameters are obtained by $\hat{\theta} = \arg \min_{\theta} \mathcal{L}(\theta)$. We denote $\theta_{-z}$ as the optimal model parameters obtained by $\arg \min_{\theta} \frac{1}{n} \sum_{i=1}^n \ell(z_i; \theta) - \frac{1}{n} \ell(z; \theta)$. We denote the Hessian of the empirical risk evaluated at $\theta$ by $H_{\theta} = \nabla_{\theta}^2 \mathcal{L}(\theta)$.

**Derivation of IF.** The key idea in Koh and Liang (2017) is to approximate the change in model parameters when a single training point $z$ is removed, i.e., $\hat{\theta}_{-z} - \hat{\theta}$, without the need for expensive retraining. To formalize this idea, consider a perturbed empirical risk that upweights $z$ by a small amount $\epsilon$:

$$ \mathcal{L}_{\epsilon,z}(\theta) = \frac{1}{n} \sum_{i=1}^{n} \ell(z_i; \theta) + \epsilon \ell(z; \theta), \quad (1) $$

and let $\hat{\theta}_{\epsilon,z}$ denote the minimizer of this perturbed objective. Here, $\epsilon = -\frac{1}{n}$ corresponds to leave-one-out (LOO) (i.e., $\hat{\theta}_{-z}$) while $\epsilon = 0$ corresponds to the original model parameters $\hat{\theta}$. When $\ell$ is strictly convex, we can apply the implicit function theorem (Krantz and Parks, 2002) and compute the local perturbation of $\hat{\theta}_{\epsilon,z}$ around $\epsilon = 0$ as

$$ \left. \frac{d\hat{\theta}_{\epsilon,z}}{d\epsilon} \right|_{\epsilon=0} = -H_{\hat{\theta}}^{-1} \nabla_{\theta} \ell(z; \hat{\theta}). $$

Building on this formula, IF leverages a first-order Taylor expansion to approximate the change in parameters as

$$ \hat{\theta}_{-z} - \hat{\theta} \approx \frac{1}{n} H_{\hat{\theta}}^{-1} \nabla_{\theta} \ell(z; \hat{\theta}). \quad (2) $$

In practice, we are often interested not only in how the parameters change but also in the impact of these changes on a specific quantity. As discussed in Section 2.1.1, this is typically called the target function, denoted as $f(\theta)$. For example, a common target function is the loss at a test point, $f(\theta) = \ell(z_{\text{test}}; \theta)$. By applying the chain rule, we propagate the parameter change's effect as follows

$$ f(\hat{\theta}_{-z}) - f(\hat{\theta}) \approx \nabla_{\theta} f(\hat{\theta})^{\top} (\hat{\theta}_{-z} - \hat{\theta}) \approx \frac{1}{n} \nabla_{\theta} f(\hat{\theta})^{\top} H_{\hat{\theta}}^{-1} \nabla_{\theta} \ell(z; \hat{\theta}). \quad (3) $$

This expression quantifies the influence of the training point $z$ on the target function $f$, encapsulating the central idea of IF as a first-order approximation to LOO retraining.
**Challenges.** Koh and Liang (2017) demonstrates that, although the derivation of IF relies on the loss function being strongly convex and the model is trained to convergence, IF can be adapted to non-convex models and still perform reasonably well. However, their experiments are restricted to small-scale neural networks. When scaling IF to larger and complex models, which are prevalent in generative AI, several challenges emerge. For example, Basu et al. (2021) shows that IF is sensitive to model size and various training hyperparameters. Meanwhile, Hammoudeh and Lowd (2022) argues that IF induces a low-loss penalty, causing confidently predicted training samples to be ranked as uninfluential. Bae et al. (2022) show that, in practice, IF in deep neural networks approximates a different measure, i.e., the proximal Bregman response function, rather than exact LOO. Finally, even putting aside these methodological concerns, the computational cost and memory footprint of IF, which involves the (inverse) Hessian and gradients for all training samples, become prohibitively large as both model and dataset sizes increase.

These challenges have spurred two research directions<fnref n="3" /> particularly relevant to generative AI: one that tackles non-convexity and non-convergence issues in large-scale models, and another that improves the scalability of IF’s. In Section 2.2.1 and Section 2.2.2, we examine these research avenues respectively. Moreover, in Section 2.2.4, we examine TRAK (Park et al., 2023), a variant of IF that has recently attracted significant attention. We conclude this section by summarizing key takeaways and outlining future work in Section 2.2.5.

### 2.2.1 Addressing Non-convexity and Non-convergence

Non-strongly convexity and non-convergence are two key bottlenecks that limit the application of IF to large-scale models. In this section, we review existing strategies aimed at mitigating these challenges and highlight several open questions for future research. We note that trajectory-specific LOO (Hara et al., 2019; Bae et al., 2024; Wang et al., 2025b), a line of work within the category of training dynamic methods, also aims to address these challenges. We discuss it in more detail in Section 2.4.2.

**Singularity of Hessian.** One of the primary challenges when applying influence functions to non-strongly convex models is that the Hessian matrix may be singular or not positive-definite. In such cases, the inverse of the Hessian, which is required for calculating the IF as in Equation (2), becomes unstable or undefined. To address this issue, a common strategy, introduced by Koh and Liang (2017), is to add a damping term $\lambda I$ to the Hessian to make it positive-definite. This damping term can be interpreted as equivalent to adding an $L_2$-regularization to Equation (1), thereby constraining the model parameters to a small neighborhood around their optimal values. This constraint not only enhances numerical stability but also reinforces the validity of the first-order approximation, thereby improving IF’s attribution quality (Schioppa et al., 2022). However, tuning the regularization coefficient is challenging without a clear ground truth; indeed, Zhang and Zhang (2022) demonstrates that IF’s accuracy is sensitive to this choice, with too small a value significantly degrading performance.

Another widely adopted approach is to replace the Hessian with an approximation with benign properties. One such choice is the Fisher Information Matrix (FIM) (Fisher, 1922); specifically, for models trained with the negative log-likelihood objective, the FIM is equivalent to Hessian in expectation (a.k.a. Bartlett’s second identity (Bartlett, 1953)), where the expectation is taken with respect to the model’s predictive distribution. Unlike the Hessian in non-convex settings, the FIM is guaranteed to be positive semi-definite. Recent work (e.g., Teso et al., 2021; Kwon et al., 2024; Zhou et al., 2024) has replaced the vanilla Hessian with the empirical FIM in IF, achieving promising results. An alternative to the FIM is the Generalized Gauss-Newton (GGN) matrix (Bae et al., 2022; Grosse et al., 2023; Mlodozieńiec et al., 2025), which is equivalent to the FIM in many practical cases (Martens, 2020). Both the FIM and GGN have been shown to offer a more robust and stable approximation of the Hessian, thereby enhancing the quality of attributions derived from IF. Additionally, the simpler formulation of the FIM as a sum of outer products of gradients leads to improved computational efficiency. We will discuss this further in Section 2.2.2.

**Randomness from model training.** Recent studies have shown that the validity of IF can be compromised by the inherent randomness in model training (K and Søgaard, 2021; Nguyen et al., 2023). This randomness stems from factors such as random initialization, batch sampling, stochastic optimization methods (e.g., dropout (Srivastava et al., 2014)), and even the limitations of floating-point precision. Such randomness can inject significant noise into the influence scores, potentially overwhelming the true signal and rendering the attributions less informative.

One common strategy to address this challenge is to use an ensemble approach (Dietterich, 2000), which smooths out noise and improves the signal-to-noise ratio by averaging influence scores across multiple models trained on the same dataset. Early work in this area involved independently trained models (K and Søgaard, 2021; Park et al., 2023;

<footnote n="3">³Beyond the research discussed here, IF’s have been extended to analyze the impact of a subset or group of samples (Koh et al., 2019; Hu et al., 2024a), incorporate higher-order information (Giordano et al., 2019a; Basu et al., 2020), and refine influence ranking through normalization techniques (Hammoudeh and Lowd, 2022; Barshan et al., 2020).</footnote>
Bae et al., 2024), while more recent studies have demonstrated that ensembling different versions of a single model, obtained through variations in dropout, fine-tuning procedures, or checkpoint selection, can further improve attribution quality (Deng et al., 2024a; Wu et al., 2024b). On the other hand, while ensembling is a convenient, plug-and-play technique, it also incurs additional computational cost.

**Non-convergence gap.** In addition to the issues discussed above, another significant concern is that neural networks may not be trained to full convergence, either intentionally (since early stopping is widely recognized to enhance generalization (Prechelt, 2002; Goodfellow, 2016)) or due to computational constraints. In this context, Bae et al. (2022) systematically analyzes this issue along with the previous ones and argue that the practical implementation of IF in deep neural networks does not approximate LOO but instead tracks a different measure, which they term the proximal Bregman response function (PBRF). This quantity can be viewed as a prediction-constrained version of LOO. Notably, while their work does not provide a solution to bridge the non-convergence gap, it offers insight into what IF is tracking and contends that this measure can still be valuable for data attribution. More recently, Wu et al. (2024b) proposed Debias and Denoise Attribution (DDA), which addresses the challenge of non-convergence in fine-tuning LLMs. Specifically, DDA improves influence function estimates through a debiasing strategy that removes bias introduced by the pretrained model and a denoising strategy that smooths variations in influence scores across multiple checkpoints.

### 2.2.2 Improving Scalability

Aside from non-convexity and non-convergence, the computational challenges associated with IF present significant bottlenecks in large-scale machine learning applications. In this section, we outline the primary hurdles in scaling IF and review techniques to address them.

Recall that computing IF involves the following inverse-Hessian-Vector Product (iHVP):

$$ \underbrace{H_{\hat{\theta}}^{-1} v}_{O(p^3)} \quad \text{where} \quad H_{\hat{\theta}} = \underbrace{\frac{1}{n} \sum_{i=1}^{n} \nabla_{\theta}^2 \ell(z_i, \hat{\theta})}_{O(np^2)} \quad \text{and} \quad v = \underbrace{\nabla_{\theta} \ell(z_{\text{test}}, \hat{\theta})}_{O(p)} $$

The standard approach of computing a single iHVP (i.e., for a single $v$) requires calculating the empirical Hessian and then inverting it. This results in a time complexity of $O(np^2 + p^3)$ and a space complexity of $O(p^2)$, which are generally prohibitive even for small neural networks. Moreover, once the iHVP is obtained, the default strategy for identifying influential samples involves scanning the entire training set and computing the per-sample gradient, then taking the dot-product with the iHVP as in Equation (3) to obtain the influence for each sample. This increases the time complexity by $O(np)$.

To reduce computational costs, existing research focuses primarily on three directions: 1) restricting the set of model parameters considered (reducing $p$); 2) subsampling the training dataset (reducing $n$); and 3) efficiently approximating the iHVP instead of computing it directly using the standard approach.<fnref n="4" /> Among these, the third direction has attracted substantially more research attention. We note that the techniques from these directions are orthogonal and can often be combined.

**Restricting model parameters.** A common approach to scaling up influence functions is to restrict the computation to a subset of model parameters. For instance, several works compute influences only on the last layer (Koh and Liang, 2017; Yeh et al., 2018; Guo et al., 2021). However, Basu et al. (2021) shows that influence functions computed on a single layer may not fully capture the overall influence of training examples.

In the context of LLMs, Grosse et al. (2023) compute influences only for MLP layers while excluding others such as the attention layers. Another promising strategy is to leverage parameter-efficient fine-tuning paradigms common in LLMs, such as Low-Rank Adaptation (LoRA) (Hu et al., 2022a). Performing attribution solely on the LoRA layers substantially reduces the number of parameters under consideration, an approach that has gained increasing popularity in recent studies (Kwon et al., 2024; Zhou et al., 2024).

**Filtering training samples.** Another strategy for scaling up influence functions is to pre-filter the training dataset to identify a subset of samples that are likely to be influential, thereby avoiding the expensive computation of influence scores for every training example. Guo et al. (2021) reduce the search space to the top-$k$ nearest neighbors of the test point in the embedding space. Grosse et al. (2023) instead leverage TF-IDF (Ramos et al., 2003), a token-matching

<footnote n="4">⁴Methods in this category also involve (more advanced) strategies of reducing $p$.</footnote>
technique, to select sequences most related to a given query, which is conceptually similar to the kNN approach of Guo et al. (2021). However, Grosse et al. (2023) also note that token-matching methods may fail to capture the higher-level abstractions inherent in language models. Other lines of work involve attributing at the level of groups of training points (Ley et al., 2024), and can wrap existing IF methods, reporting speedups based on group granularity.

**Efficiently approximating iHVP.** The two directions above do not tackle the fundamental challenge of computing the iHVP, which is the primary bottleneck in scaling IF's. Here, we discuss three approaches to efficiently approximate the iHVP: 1) using compute-friendly approximations of the Hessian; 2) employing iterative algorithms; and 3) applying gradient projection methods.

To derive the time and space complexities of different methods, we consider the model to be a multilayer perceptron (MLP) with $L$ layers, each having a weight matrix $W_l \in \mathbb{R}^{a \times a}$ of identical shape across layers. The number of parameters per layer is $d = a^2$. At layer $l$, we denote $s_l \in \mathbb{R}^a$ as the pre-activation, $a_l \in \mathbb{R}^a$ as the activation, and $\theta_l = \text{vec}(W_l)$ as the flattened weights. We further define the shorthand $g_{s_l} = \nabla_{s_l} \ell(\hat{\theta}; z)$ and $g_{\theta_l} = \nabla_{\theta_l} \ell(\hat{\theta}; z)$ to indicate loss gradients w.r.t. $s_l$ and $\theta_l$. We continue to use $n$ to denote the total number of training points.

*Compute-friendly approximations of the Hessian.* A common strategy to reduce the computational cost of iHVP is to approximate the Hessian with a more tractable matrix. As discussed in “Singularity of Hessian,” one popular alternative is the Fisher Information Matrix (FIM), which approximates the second-order Hessian with first-order gradient information only. Another widely adopted technique is to use a layer-wise block-diagonal approximation of the FIM, ignoring interactions between parameters across different layers. Concretely, the full Hessian is approximated by $\text{diag}\{G_1, \dots, G_L\}$, where

$$G_l = \mathbb{E}_z \underbrace{[\nabla_{\theta_l} \ell(\hat{\theta}; z) \quad \nabla_{\theta_l} \ell(\hat{\theta}; z)^\top]}_{O(d)} = \underbrace{\mathbb{E}_z [g_{\theta_l} g_{\theta_l}^\top]}_{O(n d^2)}.$$

We refer to this approach as *FIM + block-diagonal approximation*, or *BLOCK DIAGONAL* for short. This approach alone is capable of cutting asymptotic iHVP runtime and storage notably. Specifically, computing and inverting all $G_\ell$ results in an overall time complexity of $O(nLd^2 + Ld^3)$, as opposed to $O(nL^2d^2 + L^3d^3)$ when all $p = Ldp = Ld$ parameters are considered. Meanwhile, the Hessian and its inverse would typically require $O(p^2) = O(L^2d^2)$ in space, whereas computing these layer-wise requires a space complexity of $O(Ld^2)$ instead.<fnref n="5" />

We further discuss a few advanced extensions of this approach in several recent works.

* K-FAC (Martens and Grosse, 2015), or Kronecker-Factored Approximate Curvature, is an efficient approximation to the aforementioned FIM + block-diagonal. The key idea behind K-FAC is that the parameter gradients within each layer of an MLP can be further factored into $g_{\theta_l} = a_{l-1} \otimes g_{s_l}$, where $\otimes$ denotes the Kronecker product:

$$G_l = \mathbb{E}_z [(a_{l-1} \otimes g_{s_l}) (a_{l-1} \otimes g_{s_l})^\top] = \mathbb{E}_z [(a_{l-1} a_{l-1}^\top) \otimes (g_{s_l} g_{s_l}^\top)] \\ \approx \underbrace{\mathbb{E}_z [a_{l-1} a_{l-1}^\top]}_{O(na^2)} \otimes \underbrace{\mathbb{E}_z [g_{s_l} g_{s_l}^\top]}_{O(na^2)} = A_{l-1} \otimes S_l.$$

This final step assumes approximate independence between the activations $a_{l-1}$ and the backpropagated gradients $g_{s_l}$. While this is not strictly true, Martens and Grosse (2015) justify this simplification using heuristic arguments based on graphical models, supported empirically—and it is theoretically exact in simpler linear-Gaussian regimes (Bernacchia et al., 2018). Since the inverse of the block-diagonal matrix is also block-diagonal, the $G^{-1}v$ approximation can be performed layer-wise, and the separation of $A_{l-1}$ and $S_l$ enables the following useful application of identities, where $V_l \in \mathbb{R}^{a \times a}$ and $v_l = \text{vec}(V_l)$:

$$G_l^{-1} v_l = (A_{l-1} \otimes S_l)^{-1} v_l = \underbrace{(A_{l-1}^{-1} \otimes S_l^{-1})}_{O(a^3)} v_l = \text{vec}(\underbrace{(S_l^{-1} V_l A_{l-1}^{-1})}_{O(a^3)}).$$

Estimating $A_{l-1}$ and $S_l$ thus requires $O(na^2)$ runtime, while inversion requires $O(a^3)$. The additional matrix multiplications cost $O(a^3)$. Space requirements scale with $O(a^2)$, the size of each matrix.

Putting this together for $L$ layers, K-FAC enjoys a time complexity of $O(nLd+Ld^{3/2})$, owing to the Kronecker factorization of each layer's parameters into input and output dimensions, as well as a space complexity of $O(Ld) = O(p)$, which is the size of the model.

<footnote n="5">⁵We assume that in practice $L \ll d$, and $O(Ld)$ storage requirements to compute gradients thus vanish. Similarly, computing $G_l$ technically requires $O(nLd + nd^2) = O(nd^2)$ time when performed layer-wise.</footnote>
* EK-FAC (Grosse et al., 2023), the eigenvalue-corrected K-FAC approximation originally proposed by George et al. (2018), eigendecomposes $A_{l-1}$ and $S_l$ and efficiently computes the iHVP by inverting the resulting diagonal form. Asymptotic time and space complexities remain the same as vanilla K-FAC, though the resulting approximation is provably better. EK-FAC also provides a way to handle damped iHVPs in the form of $(G_l + \lambda_l I)^{-1}$, which is equivalent to adding $\lambda_l$ to each of the eigenvalues.

* DATAINF (Kwon et al., 2024) proposes to swap the order of operations in the expression for $(G_l + \lambda_l I)^{-1}$, from the standard inverse of the FIM expectation, to an expectation over per-sample FIM inverses, via

$$ \left( \frac{1}{n} \sum_{i=1}^{n} g_{i, \theta_l} g_{i, \theta_l}^{\top} + \lambda_l I \right)^{-1} \approx \frac{1}{n} \sum_{i=1}^{n} \left( g_{i, \theta_l} g_{i, \theta_l}^{\top} + \lambda_l I \right)^{-1}, $$

where we extend our shorthand to $g_{i, \theta_l} = \nabla_{\theta_l} \ell(\theta, z_i)$. While a theoretical justification for this rearrangement is yet to be established, Kwon et al. (2024) derive error bounds under mild assumptions (bounded gradients and positive damping), showing that the approximation error decreases with stronger damping and is controlled by the layer's parameter dimension. The main practical benefit is that it enables direct use of the Sherman-Morrison formula, yielding a much more efficient closed form expression for $(G_l + \lambda_l I)^{-1}$ of

$$ \frac{1}{n} \sum_{i=1}^{n} (g_{i, \theta_l} g_{i, \theta_l}^{\top} + \lambda_l I)^{-1} = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{\lambda_l} \left( I - \frac{g_{i, \theta_l} g_{i, \theta_l}^{\top}}{\lambda_l + g_{i, \theta_l}^{\top} g_{i, \theta_l}} \right). $$

The full influence estimate $-v_l^\top (G_l + \lambda_l I)^{-1} g_{k, \theta_l}$ for the $l^{th}$-layer reduces the expression to$^6$

$$ \frac{1}{n\lambda_l} \sum_{i=1}^{n} \underbrace{\left( \frac{(v_l^\top g_{i, \theta_l})(g_{i, \theta_l}^\top g_{k, \theta_l})}{\lambda_l + g_{i, \theta_l}^\top g_{i, \theta_l}} \right)}_{O(d)} - v_l^\top g_{k, \theta_l}, $$

which now consists only of inner-products between layer-wise gradients in $\mathbb{R}^d$ space, reducing runtime complexity to $O(n d)$ per layer. For an estimate of the influence of a single training point on a single test point, DataInf has an overall time complexity of $O(nLd)$ and space complexity of $O(Ld)$, i.e., for each gradient term in the expression.

*Iterative algorithms.* Another major source of computational cost is the inversion of the Hessian (or its approximation). Even with block-diagonal approximations, direct inversion can be expensive. Iterative algorithms have been developed to compute the iHVP more efficiently.

* LISSA: To approximate $H_{\hat{\theta}}^{-1} v$, Koh and Liang (2017) employs LiSSA, an algorithm for stochastic optimization proposed in Agarwal et al. (2017)). Given a set of training batches $z_{s_1}, \dots, z_{s_t}$ sampled uniformly at random, the method iterates

$$ \tilde{H}_0^{-1} v = v; \quad \tilde{H}_j^{-1} v = v + (I - \nabla_{\theta}^2 \ell(z_{s_j}, \hat{\theta})) \tilde{H}_{j-1}^{-1} v, $$

for $j = 1, \dots, t$. This procedure is repeated $r$ times, and the results are averaged to yield the final iHVP approximation. The efficiency of LiSSA partly arises from the fact that the Hessian-vector product $\nabla_{\theta}^2 \ell(z, \hat{\theta}) v$ can be computed in $O(p)$ time (Pearlmutter, 1994), for $O(rt)$ samples. Practically, Koh and Liang (2017) note that $rt \approx n$ in their settings gives accurate results, resulting in an overall online time complexity of $O(np)$ and space complexity of $O(p)$. Subsequent works have reduced the effective number of repetitions $r$ by parallelizing the computation (Guo et al., 2021) or using advanced sampling techniques (e.g., centroid-based sampling via $K$-means) to improve the robustness of the Hessian approximation (Koh et al., 2024).

* SCHULZ: Zhou et al. (2024) adopt Schulz's iteration (Herzberger and Petković, 1990) to compute the inverse of the Hessian approximation $G$. The (per-layer) iterative update is given by

$$ \tilde{G}_0^{-1} \approx G^{-1}; \quad \tilde{G}_j^{-1} = \underbrace{\tilde{G}_{j-1}^{-1} (2I - G \tilde{G}_{j-1}^{-1})}_{O(d^3)}. $$

Although it does not improve the asymptotic complexity of matrix inversion, this method comes with numerical stability guarantees and is relatively insensitive to the initial guess.<fnref n="7" /> Computing $G$ requires $O(n d^2)$ time per

<footnote n="6">⁶N.B. the original DataInf notation has been heavily adapted to match the format presented in this section.</footnote>

<footnote n="7">⁷Initializing with the identity or a random gaussian matrix resulted in convergence within $t < 20$ steps according to Zhou et al. (2024).</footnote>
layer, and further applying t iterative updates as shown above at $O(d^3)$ time across L layers results in a time complexity of $O(nLd^2 + tLd^3)$ and a space complexity of $O(d^2)$.

**Gradient projection.** A third approach to reduce iHVP cost is to project the gradients onto a lower-dimensional subspace, thereby reducing the size of the matrices involved in the computation. We denote $\tilde{p}$ and $\tilde{d}$ as the projected number of parameters, a fraction of the original number of parameters $p$ and $d$. We do likewise for projected activations $\tilde{s}_l$ and $\tilde{a}_{l-1}$. Two common techniques are:

* RANDOM (Schioppa et al., 2022; Park et al., 2023): In this approach, a projection matrix $P \in \mathbb{R}^{\tilde{p} \times p}$ is drawn from a suitable distribution (e.g., Gaussian), where $\tilde{p} \ll p$ represents the reduced dimensionality. Replacing gradients $g_\theta$ with projected gradients $Pg_\theta$, the Hessian is approximated by:

$$ \tilde{G} = \frac{1}{n} \sum_{i=1}^{n} \underbrace{(Pg_{\theta})(Pg_{\theta})^{\top}}_{O(p\tilde{p}+\tilde{p}^2)}, $$

which rearranges to $\tilde{G} = PGP^{\top}$, and can be interpreted as a restriction of the Hessian to a subspace spanned by the rows of $P$. Likewise, train and test gradients are projected onto this space in $O(p\tilde{p})$ time and space to compute iHVPs and influence scores.

The cost of computing and inverting $\tilde{G}$ is $O(np\tilde{p}+\tilde{p}^3)$ which simplifies to $O(np\tilde{p})$ assuming $\tilde{p} \ll p, n$. Should this be performed layerwise with the block diagonal approximation, we have $O(nLd + nLd\tilde{d})$ runtime, as gradients must still be computed for the entire model using $O(Ld)$ time. Space scales with $O(Ld + d\tilde{d})$ owing to an initial full gradient computation followed by projection occurring at the layer level.

* ARNOLDI: While the Johnson-Lindenstrauss lemma (William and Lindenstrauss, 1984) ensures that the inner products between gradients are approximately preserved under random i.i.d. projection, if the projection subspace is not invariant under $G$ (i.e., $G$ applied to vectors in this subspace generates significant components orthogonal to the subspace), the projected Hessian $PGP^{\top}$ implicitly neglects these orthogonal components, which leads directly to approximation errors. This motivates Arnoldi iteration. Schioppa et al. (2022) propose to use Arnoldi iteration (Arnoldi, 1951) to approximate the top-$\tilde{p}$ (in magnitude) eigenvalues of $H_\theta$, and use the corresponding eigenvectors as the rows of the projection matrix $P \in \mathbb{R}^{\tilde{p} \times p}$. Since this subspace is spanned by the eigenvectors of $H_\theta$, it is by definition $H_\theta$-invariant, i.e., vectors will remain in the subspace when transformed by $H_\theta$. Another highly desirable property of top-$\tilde{p}$ eigenvector projection is that $PH_\theta P^{\top}$ simplifies to the diagonal matrix of eigenvalues, $\Lambda$, which is trivially invertible in $O(\tilde{p})$ time. The iHVP then becomes

$$ (PH_{\tilde{\theta}}P^{\top})^{-1}(Pv_l) = \underbrace{\Lambda^{-1}}_{O(\tilde{p})} \underbrace{Pv_l}_{O(p\tilde{p})}. $$

While projecting onto the top-$\tilde{p}$ eigenvectors of the Hessian matrix is highly desirable$^8$, this step comes with extra computational cost, as we will now outline. The Arnoldi iteration starts by iteratively constructing the following $\tilde{p}^{\text{th}}$ order Krylov subspace$^9$ of the Hessian:

$$ \text{Span}\{v, H_{\tilde{\theta}}v, \underbrace{H_{\tilde{\theta}}^2 v, \dots, H_{\tilde{\theta}}^{\tilde{p}} v}_{O(Bp)}\}, $$

where $v \in \mathbb{R}^p$ is randomly drawn. Intuitively, the sequence above naturally converges to the largest eigenvector of $H_\theta$, and for sufficiently large $\tilde{p}$, we may also expect an orthonormal basis of the subspace to span good approximations of the top-$\tilde{p}$ eigenvectors of $H_\theta$. While iterating on the above, Arnoldi orthogonalizes each new term w.r.t. each of the previous terms, and subsequently normalizes, to form said orthonormal basis. It therefore uses $\tilde{p}$ iterations, each of which requires an HVP and its orthogonalization with respect to every existing basis vector, costing $O(Bp\tilde{p} + p\tilde{p}^2)$ runtime. Note that the HVP is not computed over the full training data, nor individual training samples as in LiSSA, but rather over batches of size $B$. Concurrently, a matrix $A$ is computed containing information about inner products between each new HVP and the basis vectors up until that iteration. Schioppa et al. (2022) then eigendecomposes $A$ and projects the eigenvectors onto the orthonormal basis to approximate the top-$\tilde{p}$ eigenvectors of $H_\theta$. Eigendecomposing $A$ entails $O(\tilde{p}^3)$ runtime,

<footnote>$^8$Theoretically, the top-$\tilde{p}$ eigenpairs also yield the best rank-$\tilde{p}$ approximation of the full (symmetric) Hessian: specifically, $P^{\top}\tilde{H}P$ is closest in both the Frobenius and spectral norms to $H_\theta$, when $P \in \mathbb{R}^{\tilde{p} \times p}$ consists of the top-$\tilde{p}$ eigenvectors, and $\tilde{H} \in \mathbb{R}^{\tilde{p} \times \tilde{p}}$ is the diagonal matrix of corresponding eigenvalues (Eckart and Young, 1936).</footnote>

<footnote>$^9$Technically speaking, Schioppa et al. (2022) constructs a $> \tilde{p}^{\text{th}}$ order subspace first, before selecting the top-$\tilde{p}$ eigenpairs.</footnote>
while projecting gradients with P requires $O(p\tilde{p})$, meaning the total runtime remains at $O(Bp\tilde{p} + p\tilde{p}^2)$, or equivalently $O(BLd\tilde{d} + Ld\tilde{d}^2)$ in the block-diagonal setting. Space scales with $O(p\tilde{p} + p\tilde{p}^2)$, or $O(Ld\tilde{d} + Ld\tilde{d}^2)$.

* **LOGRA** (Choe et al., 2024): Low-Rank Gradient Projection improves space and time complexity of gradient projection by separating layer-wise projections into two smaller projections at the input and output of the layer (similarly exploited in K-FAC). In doing so, the full projected gradient may be rewritten as the following, where $P_i$ and $P_o$ represent the projections applied to the input activations and output pre-activations, respectively:

$$Pg_{\theta_l} = (P_i \otimes P_o)(a_{l-1} \otimes s_l) = P_i a_{l-1} \otimes P_o s_l = \tilde{a}_{l-1} \otimes \tilde{s}_l,$$

or the Kronecker product of projected activations/gradients in the layer, where $\tilde{a}_{l-1}$, $\tilde{s}_l \in \mathbb{R}^{\tilde{a}}$ represent projected activations. As for the Hessian estimate $\tilde{G}_l = PG_lP^\top$, applying the same identities used to derive K-FAC yields:

$$ \tilde{G}_l = \mathbb{E}_z[(Pg_{\theta_l})(Pg_{\theta_l})^\top] = \mathbb{E}_z[(\tilde{a}_{l-1} \otimes \tilde{s}_l)(\tilde{a}_{l-1} \otimes \tilde{s}_l)^\top] \\ \approx \mathbb{E}_z[\tilde{a}_{l-1}\tilde{a}_{l-1}^\top] \otimes \mathbb{E}_z[\tilde{s}_l\tilde{s}_l^\top] = \tilde{A}_{l-1} \otimes \tilde{S}_l, $$

with forward and backward uncentered covariance matrices $\tilde{A}_{l-1}$, $S_l \in \mathbb{R}^{\tilde{a} \times \tilde{a}}$. Overall, the derivations in this section are directly analogous to those in K-FAC, with the final iHVP form of $\tilde{G}_l^{-1}\tilde{v}_l = \text{vec}(\tilde{S}_l^{-1}\tilde{V}_l\tilde{A}_{l-1}^{-1})$, where $\tilde{v}_l = P\tilde{v}_l$ is the projected test gradient w.r.t. the layer's parameters. We denote the matrix form (taking input and output shapes into account) of $\tilde{v}_l$ as $\tilde{V}_l$. Projecting onto each of these spaces separately requires $O(\sqrt{d\tilde{d}})$ time, in contrast to standard projection of $g_{\theta_l} \in \mathbb{R}^d$ to $\tilde{g}_{\theta_l} \in \mathbb{R}^{\tilde{d}}$ which requires $O(d\tilde{d})$. The full runtime complexity of LoGra is $O(nLd/B + nL\sqrt{d\tilde{d}} + L\sqrt{\tilde{d}}^3)$, owing to per-sample gradient computation (batch size $B$), the computation of $\tilde{A}_{l-1}$ and $\tilde{S}_l$ followed by inversion of $\tilde{G}_l$ for each layer. Assuming $\sqrt{\tilde{d}} \ll n$, the first term dominates, giving $O(nL\sqrt{d\tilde{d}})$ runtime for iHVP pre-computation. As for space complexity, we still require $O(Ld)$ to store the model, which dominates the smaller space requirements of $O(\sqrt{d\tilde{d}})$ for the projection matrices and $O(\tilde{d})$ for matrix inversions. If projected training gradients are cached for reuse, an additional disk storage cost of $O(nL\sqrt{\tilde{d}})$ is incurred. In summary, LoGra uses a combination of gradient projection, block-diagonal approximation, and K-FAC approximation.

### 2.2.3 Comparison of IHVP Computational Complexities

To clarify comparisons across methods, Table 3 presents complexities for computing the influence of a single training point on a single test point, $v$. We distinguish clearly between **pre-compute** complexity, to represent any Hessian related computational cost (such as IHVPs, $H_{\hat{\theta}}^{-1}v$, or inverse Hessians, $H_{\hat{\theta}}^{-1}$), and **per-vector (online)** complexity, the cost of attributing to each training point, given access to the pre-computed data relating to the specific test point. We assume that the methods utilize the FIM + block-diagonal approximation of the Hessian wherever possible (i.e., besides the standard LISSA method). Note that in many cases, parallelization can exchange runtime complexity for space complexity $O(nLd)$ for increased space complexity $O(nLd)$, beneficial in large-scale settings, i.e. $L$ can be transferred from runtime to space complexity if parallelized, or online runtime complexity $O(nLd)$ can be interchanged with increased space complexity.

* **EK-FAC**, **SCHULZ**, and **BLOCK DIAGONAL** exploit structured approximations, significantly reducing complexities from exact methods $O(nL^2d^2 + L^3d^3)$ to $O(nLd + Ld^{3/2})$, $O(nLd^2 + tLd^3)$, and $O(nLd^2 + Ld^3)$, respectively. Assuming terms with $n$ dominate in practical settings, the $O(nL^2d^2)$ pre-compute complexity for exact influence is reduced to $O(nLd^2)$ with **BLOCK DIAGONAL** and **SCHULZ**, and further reduced to $O(nLd)$ with EK-FAC.

* **DATAINF** shifts the IHVP computation primarily online. However, it can be structured to transfer complexity offline, reducing the online complexity to $O(Ld)$ while incurring $O(nLd)$ offline costs by factoring out gradient terms, denoted with DATAINF (W/PRE) in Table 3. This is necessary if attributing to multiple training points, allowing online time complexity to scale with $O(nLd)$ instead of $O(n^2Ld)$, as we further describe below.

* **LISSA** typically moves computation online with complexity $O(rtLd)$. However, pre-computation is possible by caching intermediate vectors ($s_{\text{test}} = H_{\hat{\theta}}^{-1}v_{\text{test}}$), allowing the online runtime to be reduced to $O(Ld)$ per training point, or $O(nLd)$ across multiple training points. To compute importance for multiple test points, LISSA's pre-compute complexity grows with $m$, the number of test points, to $O(mrtLd)$, as does online runtime, which scales linearly with $m$.
* Projection methods (RANDOM, ARNOLDI, LOGRA) further minimize pre-compute complexity, with slight increases in per-vector complexity, making them indispensable in applications where full gradient computations are prohibitive (e.g., LLMs). Of these, LoGra offers the most extensive speedups, combining principles from EK-FAC (factorization of layerwise gradients into forward and backward passes) with gradient projection to achieve $O(nL\sqrt{d\tilde{d}})$ pre-compute complexity.

**Attributing to multiple training points.** With $p = Ld$ parameters, attributing to multiple training points increases online compute from $O(Ld)$ to $O(nLd)$ for all methods (this may be transferred to space complexity when utilizing parallelization). Online cost may further increase to $O(mnp)$ if scores are required across $m$ test points.

**Pre-compute costs for multiple test points.** While computing the influence of multiple training points increase online runtime, computing these training data influences on multiple test points has varying effects on the pre-compute complexities across the methods. In general, methods that compute IHVPs directly typically scale poorly to multiple test points, as they save little re-usable information such as the result of $H_{\tilde{\theta}}^{-1}$.

* EXACT influence adds $O(mp^2)$ to compute $H_{\tilde{\theta}}^{-1}v$ for all $m$ test points, since $H_{\tilde{\theta}}^{-1}$ is explicitly computed. This new term remains dominated by the original $O(p^3)$ terms, and complexities are unchanged. In layer-wise settings, this increase is $O(mLd^2)$, which vanishes for BLOCK DIAGONAL and SCHULZ, where the pre-compute time complexities include larger $O(nLd^2)$ terms.

* EK-FAC requires $O(Ld^{3/2})$ pre-compute runtime for each new IHVP requiring. Thus, the precompute complexity with $m$ IHVP operations is increased to $O(nLd + mLd^{3/2})$, which remains unchanged if $m\sqrt{d} \ll n$.

* DATAINF allows for factorization of $v_{\ell}^{\top}$ and $g_{k,\theta_{\ell}}$ terms. As such, the remaining expression comprising $n$ training gradients can first be computed in $O(nLd)$ time. Pre-computing all IHVPs comes with an additional $O(mLd)$ pre-compute cost, which again vanishes.

* LISSA solves each particular IHVP directly, in $O(rtLd)$ time, which scales poorly to $O(mrtLd)$ for multiple test points, since, as noted above, there is little intermediate, re-usable information across test points.

* Projection methods (RANDOM, ARNOLDI, LOGRA) explicitly compute the (projected) Hessian inverse and thus incur negligible additional pre-compute costs for multiple IHVPs (as $n$ dominates $m$).

**Recommendations.** Traditional influence function research has primarily focused on addressing the computational bottleneck of Hessian inverse computation, which represents the main challenge for training runs with small-scale datasets and many training epochs. However, modern large-scale deep learning, especially for large language models,

<table><caption>Table 3: Comparison of time and space complexity of different methods to compute the influence of a single training point on a single test point. The model is an MLP with $L$ layers of identical shape, each containing $d$ parameters (or $\tilde{d}$ projected parameters).</caption><thead><tr><th>Category</th><th>Method</th><th>Pre-compute Time</th><th>Pre-compute Space</th><th>Time</th><th>Space</th></tr></thead><tbody><tr><td>Exact</td><td>NAIVE</td><td>O(nL<sup>2</sup>d<sup>2</sup> + L<sup>3</sup>d<sup>3</sup>)</td><td>O(L<sup>2</sup>d<sup>2</sup>)</td><td>O(Ld)</td><td>O(Ld)</td></tr><tr><td rowspan="4">Hessian Approx.</td><td>BLOCK DIAGONAL</td><td>O(nLd<sup>2</sup> + Ld<sup>3</sup>)</td><td>O(d<sup>2</sup>)</td><td>O(Ld)</td><td>O(Ld)</td></tr><tr><td>EK-FAC</td><td>O(nLd + Ld<sup>3/2</sup>)</td><td>O(Ld)</td><td>O(Ld)</td><td>O(Ld)</td></tr><tr><td>DATAINF</td><td>-</td><td>-</td><td>O(nLd)</td><td>O(Ld)</td></tr><tr><td>DATAINF (W/PRE)</td><td>O(nLd)</td><td>O(Ld)</td><td>O(Ld)</td><td>O(Ld)</td></tr><tr><td rowspan="2">Iterative</td><td>LISSA/FASTIF</td><td>-</td><td>-</td><td>O(rtLd)</td><td>O(Ld)</td></tr><tr><td>SCHULZ</td><td>O(nLd<sup>2</sup> + tLd<sup>3</sup>)</td><td>O(d<sup>2</sup>)</td><td>O(Ld)</td><td>O(Ld)</td></tr><tr><td rowspan="7">Projection</td><td>RANDOM</td><td>O(nLd&#x0303;)</td><td>O(Ld&#x0303;)</td><td>O(Ld&#x0303;)</td><td>O(Ld&#x0303;)</td></tr><tr><td>RANDOM (W/PRE)</td><td>O(nLd&#x0303;)</td><td>O(nLd&#x0303; + Ld&#x0303;)</td><td>O(Ld&#x0303;)</td><td>O(Ld&#x0303;)</td></tr><tr><td>ARNOLDI</td><td>O(BLd&#x0303; + Ld&#x0303;<sup>2</sup>)</td><td>O(Ld&#x0303; + Ld&#x0303;<sup>2</sup>)</td><td>O(Ld&#x0303;)</td><td>O(Ld&#x0303;)</td></tr><tr><td>ARNOLDI (W/PRE)</td><td>O(BLd&#x0303; + Ld&#x0303;<sup>2</sup> + nLd&#x0303;)</td><td>O(nLd&#x0303; + Ld&#x0303; + Ld&#x0303;<sup>2</sup>)</td><td>O(Ld&#x0303;)</td><td>O(Ld&#x0303;)</td></tr><tr><td>LOGRA</td><td>O(nL&#x0303;√d&#x0303;)</td><td>O(Ld)</td><td>O(Ld + √d&#x0303;)</td><td>O(Ld)</td></tr><tr><td>LOGRA (W/PRE)</td><td>O(nL&#x0303;√d&#x0303;)</td><td>O(nL&#x0303;√d&#x0303; + Ld)</td><td>O(L&#x0303;√d&#x0303;)</td><td>O(Ld)</td></tr></tbody></table>
presents additional challenges. Simply improving Hessian inverse efficiency is insufficient; the critical bottleneck has extended to computing dot-products between per-sample gradients and the iHVP, which requires per-sample gradient information. The gradient decomposition techniques from K-FAC literature can significantly improve the computational efficiency of this process. For contemporary ML applications requiring real-time influence estimation across numerous inference queries, we recommend a precomputation-storage paradigm: computing and storing iHVPs for each training sample to disk, then performing efficient dot-products with test gradients at inference time. However, this approach demands substantial storage capacity proportional to the training set size and model dimensionality. To mitigate storage requirements while preserving influence estimation quality, projection methods can be used to compress iHVPs into lower-dimensional representations. This storage-projection framework is so far the most practical approach for scalable influence computation in production environments, trading one-time preprocessing costs and storage overhead for dramatically reduced inference-time latency. Moving forward, there are promising opportunities to develop more effective projection techniques that maximize influence estimation accuracy under constrained storage budgets.

### 2.2.4 TRAK: A Variant of Influence Function

Recently, TRAK (Park et al., 2023) has emerged as one of the most popular data attribution approaches, finding successful application in a range of generative models due to its effectiveness and scalability. While TRAK is proposed as an efficient estimator of datamodels (Ilyas et al., 2022) (discussed in Section 2.5), it is closer to IF on a methodology level. In particular, we begin by reviewing the concept of the one-step Newton update—a fundamental idea underlying TRAK that has long been recognized for its close connection to IF (Pregibon, 1981; Wojnowicz et al., 2016; Koh et al., 2019). Next, we discuss the key design choices behind TRAK. Finally, we examine its applications across various generative models.

**One-step Newton update.** At the core of TRAK is the one-step Newton update, detailed below. Let $\mathcal{L}_{-z}$ denote the loss function computed without the point $z$, and let $H_{-z}$ be its corresponding Hessian. Starting from the full-data estimate $\hat{\theta}$, the Newton-Raphson method (Ypma, 1995) updates the parameters as

$$ \theta_{NS} := \hat{\theta} - H_{-z}^{-1} \nabla_{\theta} \mathcal{L}_{-z}(\hat{\theta}), \quad (4) $$

where, by the first-order optimality condition $\nabla_{\theta}\mathcal{L}(\hat{\theta}) = 0$, we have

$$ \nabla_{\theta}\mathcal{L}_{-z}(\hat{\theta}) = -\frac{1}{n}\nabla_{\theta}\ell(z, \hat{\theta}). $$

Comparing Equation (4) with Equation (2), we observe that the IF and one-step Newton update serve as approximations of one another, assuming that $H_{-z,\hat{\theta}}^{-1} \approx H_{\hat{\theta}}^{-1}$.

**Design choices of TRAK.** We first consider the original formulation of TRAK based on the one-step Newton update in Equation (4), and then describe the key design choices that enable efficient and scalable attribution for large-scale models.

For clarity, we derive TRAK in the context of a multi-class classification problem, which also forms the basis for its application to language models. Here, the loss function $\ell$ is the cross-entropy loss, and $\mathcal{L}$ represents the sum (rather than the average) of the loss over the training set $D$. Let $p_i$ denote the predicted probability of the correct class for sample $z_i$, and consider the target function

$$ f(z_{\text{test}}; \theta) = \log \frac{p}{1-p}, $$

referred to as the multi-class margin function (Engstrom et al., 2024) for a test point $z_{\text{test}}$. This function is particularly useful because the cross-entropy loss can be succinctly expressed as $\ell(y_i, f(z_i; \theta)) = \log(1 + e^{-f(z_i; \theta)})$. We also denote $g_i = \nabla_{\theta}f(z_i; \hat{\theta})$ and $g_{\text{test}} = \nabla_{\theta}f(z_{\text{test}}; \hat{\theta})$. It is straightforward to compute that

$$ \begin{aligned} \nabla_{\theta}\ell(z_i, \hat{\theta}) &= \frac{\partial\ell(y_i, f(z_i; \hat{\theta}))}{\partial f(z_i; \hat{\theta})} g_i = (p_i - 1)g_i, \\ H_{-z_i, \hat{\theta}} &= \sum_{j \neq i} \left( \frac{\partial\ell(y_j, f(z_j; \hat{\theta}))}{\partial f(z_j; \hat{\theta})} \nabla_{\theta}^2 f(z_j; \hat{\theta}) + \frac{\partial^2\ell(y_j, f(z_j; \hat{\theta}))}{\partial f(z_j; \hat{\theta})^2} g_j g_j^{\top} \right) \\ &= \sum_{j \neq i} \left( \frac{\partial\ell(y_j, f(z_j; \hat{\theta}))}{\partial f(z_j; \hat{\theta})} \nabla_{\theta}^2 f(z_j; \hat{\theta}) + p_j(1 - p_j)g_j g_j^{\top} \right). \end{aligned} $$
With these quantities, the ideal TRAK score of sample $z_i$ for a test sample $z_{test}$ is

$$ \mathcal{I}_{\text{TRAK}}^{\text{ideal}}(z_{\text{test}})_i = g_{\text{test}}^{\top} \left( \sum_{j \neq i} \frac{\partial \ell(y_j, f(z_j; \hat{\theta}))}{\partial f(z_j; \hat{\theta})} \nabla_{\theta}^2 f(z_j; \hat{\theta}) + p_j(1 - p_j)g_j g_j^{\top} \right)^{-1} g_i(p_i - 1). \quad (5) $$

We now discuss the key design choices adopted by Park et al. (2023).

*Linearize the target function.* The first step is to linearize the target function around the optimal model parameter $\hat{\theta}$. Formally, TRAK uses the approximation

$$ f(z_i; \theta) \approx f(z_i; \hat{\theta}) + \nabla_{\theta} f(z_i; \hat{\theta})^{\top} (\theta - \hat{\theta}). $$

This linearization effectively eliminates the potentially non-PSD component of the Hessian in Equation (5), since $\nabla_{\theta}^2 f(z_j; \hat{\theta}) = 0$. This step is equivalent to the FIM approximation discussed in Sections 2.2.1 and 2.2.2. Using the Sherman-Morrison formula (see Hu et al. (2024a) for details), the TRAK score simplifies to

$$ \frac{g_{\text{test}}^{\top} (G^{\top} R G)^{-1} g_i (1 - p_i)}{1 - h_i}, $$

where $G = [g_1, \dots, g_n]^\top \in \mathbb{R}^{n \times p}$ is the stacked gradients, $R$ is a diagonal matrix with entries $p_i(1 - p_i)$, and $h_i = 1 - g_i^{\top} (G^{\top} R G)^{-1} g_i p_i (1 - p_i)$ is the leverage score for sample $z_i$.

*Omit scaling terms.* Park et al. (2023) find empirically that both the leverage score in the denominator and the diagonal matrix $R$ have little effect on the final estimates, and thus omit them from the estimator. Ignoring the leverage score essentially reverts the formula to that of IF. The TRAK score thus further simplifies to

$$ g_{\text{test}}^{\top} (G^{\top} G)^{-1} g_i (p_i - 1). $$

*Random projection.* To reduce memory and computational costs, TRAK leverages random projection. Rather than storing and computing with the full high-dimensional gradient $g_i \in \mathbb{R}^p$, TRAK computes projected gradients $\phi_i = P^{\top} g_i$ and $\phi_{\text{test}} = P^{\top} g_{\text{test}}$, where $P \sim \mathcal{N}(0, 1)^{p \times k}$. The score for sample $z_i$ is then given by

$$ \phi_{\text{test}}^{\top} (\Phi^{\top} \Phi)^{-1} \phi_i (p_i - 1), $$

with $\Phi = [\phi_1, \dots, \phi_n]^\top \in \mathbb{R}^{n \times k}$ representing the stacked projected gradients. In practice, $k$ is chosen to be on the order of tens of thousands—substantially smaller than the original parameter count—thereby greatly enhancing computational efficiency. In fact, random projection allows TRAK to directly compute the inverse matrix without relying on additional acceleration techniques discussed in Section 2.2.2. Nonetheless, while this method boosts efficiency without significant performance loss empirically, the justification based on the Johnson-Lindenstrauss lemma, as noted by Schioppa et al. (2022), warrants further examination (see “gradient projection” in Section 2.2.2 for discussion).

*Ensemble.* Finally, TRAK employs an ensemble approach to enhance robustness. By averaging scores from different reference models, either retrained on different subsets or using distinct checkpoints, the method mitigates the inherent randomness of individual model training. Formally, given $M$ reference models $\{\hat{\theta}_m\}_{m=1}^M$ and $M$ corresponding projection matrices $\{P_m \sim \mathcal{N}(0, 1)^{p \times k}\}_{m=1}^M$, the actual TRAK score for a test sample $z_{test}$ is defined as

$$ \mathcal{I}_{\text{TRAK}}(z_{\text{test}}) = \left( \frac{1}{M} \sum_{m \in [M]} \left( \phi_{m, \text{test}}^{\top} (\Phi_m^{\top} \Phi_m)^{-1} \Phi_m^{\top} \right) \right) \cdot \left( \frac{1}{M} \sum_{m \in [M]} Q_m \right), \quad (6) $$

where

$$ \phi_{m,i} = P_m^{\top} \nabla_{\theta} f(z_i; \hat{\theta}_m), \quad \phi_{m, \text{test}} = P_m^{\top} \nabla_{\theta} f(z_{\text{test}}; \hat{\theta}_m), \quad \Phi_m = [\phi_{m,1}, \dots, \phi_{m,n}], $$

and $Q_m = \text{diag}(p_{i,m} - 1)$ with $p_{i,m}$ being the predicted probability of the correct class for sample $z_i$ according to $\hat{\theta}_m$. This $Q_m$ matrix is typically referred to as the *output-to-loss gradient* (Park et al., 2023), defined by

$$ Q = \text{diag} \left( \left\{ \frac{\partial \ell(y_i, f(z_i; \hat{\theta}))}{\partial f(z_i; \hat{\theta})} \right\} \right). $$

Note that the original implementation of TRAK constructs the $Q$ matrix with diagonal entries $1 - p_i$, differing in sign from the actual output-to-loss gradient. This is because the target function here, the multi-class margin function, is maximized (i.e., higher values are better), as opposed to common target functions in IF, such as the test loss, where lower values are preferable. Consequently, Park et al. (2023) consider the negative LOO difference $f(z_{test}; \hat{\theta}) - f(z_{test}; \hat{\theta}_{-z_i})$ so that the interpretation of TRAK scores align with the convention in data attribution: a positive score indicates a helpful sample and a negative score indicates a harmful one. For target functions that we seek to minimize (see Table 5 for examples), no sign modification is necessary.
**Applications of TRAK in generative models.** TRAK can be extended to generative tasks by designing a proper target function for different loss functions of various generative tasks.

*Language models.* The core challenge and research direction to apply TRAK on language models is how to scale to large (e.g., billion-level parameters) models and large datasets while preserving high effectiveness.

DsDm (Engstrom et al., 2024) utilizes the multi-class margin for token predictions, following the approach in TRAK, to attribute examples. To scale to large settings, the paper adopts an approach to use a small proxy model to compute the data scores and then apply them to downstream tasks on larger models. DsDm (Engstrom et al., 2024) utilizes TRAK as an estimator for linear datamodels on data selection tasks. Specifically, DsDm calculates individual data scores using a small proxy model, i.e., a pre-trained 125M parameter GPT-2 style language model on 6 billion tokens, and selects the top-K positively influential samples for training larger models with a maximum budget of 1.3B parameters.

Another study, TrackStar, combines several techniques to refine the computation of the gradient-based influence score and scales the computation to large setups of 8B pertaining LLMs. Specifically, TrackStar defines the influence score as

$$ \mathcal{I}(z, z_q) = \bar{G}_{\theta}(z) \cdot \bar{G}_{\theta}(z_q), \qquad (7) $$

where $\bar{G}_{\theta}(z)$ is the projected, Hessian-corrected, and unit-normalized gradient for example $z$, given model parameters $\theta$:

$$ \bar{G}_{\theta}(z) = \frac{G_{\theta}(z)}{\|G_{\theta}(z)\|_{2}} \qquad G_{\theta}(z) = R^{-\frac{1}{2}} P_{d} \frac{\nabla_{\theta} \ell(z, \theta)}{\sqrt{V}}. $$

where $V = \mathbb{E}_z[(\nabla_{\theta} \text{Loss}(z, \theta))^2]$, $R^{-1/2} \in \mathbb{R}^{d \times d}$ is the Hessian matrix using Gaussian-Newton approximation under random projection, and $\ell$ is defined as negative log likelihood rather than margin function.

*Diffusion models.* Several works focus on adapting TRAK to diffusion models (Ho et al., 2020a). The most active research direction focuses on the design of the target function $f_A$ due to the non-trivial adaptation of TRAK designed for classification tasks to diffusion models.

Georgiev et al. (2023) introduces JourneyTrak to attribute the conditional distribution $p_{\theta}(\mathbf{x} | \mathbf{x}_t)$ of image $\mathbf{x}$ generated by diffusion models conditioned on the intermediate latent states at each step $t$. It takes the training loss $\mathcal{L}_{\text{simple}}(\mathbf{x}, \theta) = \mathbb{E}_{\epsilon,t} [\|\epsilon - \epsilon_{\theta}(\sqrt{\bar{\alpha}_t}\mathbf{x} + \sqrt{1 - \bar{\alpha}_t}\epsilon, t)\|_2^2]$ proposed in DDPM (Ho et al., 2020a) as the target functions.

Additionally, Zheng et al. (2023) introduces D-Trak and empirically demonstrates that replacing $f_A$ with alternative norm-based functions when computing $\phi_m$ in Equation equation 6, while keeping $\mathcal{L}_{\text{simple}}$ as the loss function $\mathcal{L}$ and using $\mathcal{L}_{\text{simple}}$ for computing $Q_k$, can significantly improve performance compared to TRAK, in terms of counterfactual evaluation and Linear Datamodeling Scores (LDS). Recently, Lin et al. (2025) points out that the shift of loss function $\mathcal{L}_{\text{simple}}$ fails to capture the distributional changes of generated images, as it conducts an indirect distributional comparison. Therefore, Lin et al. (2025) proposes Diffusion Attribution Score (DAS) by further revising $f_A$ to be the noise predictor output $\epsilon_{\theta}(\mathbf{x}, t)$ of the generated sample and directly attributes the KL-divergence between the predicted distributions $D_{\text{KL}}(p_{\theta}(\mathbf{x})\|p_{\theta_{-1}}(\mathbf{x}))$, which mitigates the indirect comparison caused by using $\mathcal{L}_{\text{simple}}$ as $f_A$. Lin et al. (2025) also derives the parameter change for diffusion model trained on $\mathcal{L}_{\text{simple}}$, which coincidences with Mlodoeniec et al. (2025) and further improves the LDS. We summarize target functions used in these methods in Table 5.

*Multimodal tasks.* TRAK can also be applied to attribute multimodal tasks. For example, enabling an understanding of which training data contribute to the image-text pair association. For CLIP models, Park et al. (2023) first propose to reformulate CLIP loss as a classification problem,

$$ \ell(x_i, y_i; \theta) = -\log p_1(x_i, y_i; \theta) - \log p_2(x_i, y_i; \theta), $$

where $p_1(x_i, y_i; \theta)$ and $p_2(x_i, y_i; \theta)$ are defined as:

$$ p_1(x_i, y_i; \theta) = \frac{\exp(S_{ii})}{\sum_{1 \le j \le n} \exp(S_{ij})}, \qquad p_2(x_i, y_i; \theta) = \frac{\exp(S_{ii})}{\sum_{1 \le j \le n} \exp(S_{ji})}. $$

Here, $(x_i, y_i)$ denote an image-text pair and $S_{ij}$ denotes the pairwise similarity score between the $i$-th image and $j$-th text embeddings, typically computed using cosine similarity. Then they propose the target function defined as

$$ \begin{aligned} f_A(x_i, y_i; \theta) &= \log \frac{p_1(x_i, y_i; \theta)}{1 - p_1(x_i, y_i; \theta)} + \log \frac{p_2(x_i, y_i; \theta)}{1 - p_2(x_i, y_i; \theta)} \\ &= -\log \sum_{1 \le j \le n} \exp(S_{ij} - S_{ii}) - \log \sum_{1 \le j \le n} \exp(S_{ji} - S_{ii}). \end{aligned} $$

TRAK demonstrates strong empirical performance while significantly reducing the computation and memory cost compared to linear datamodels and original influence functions.
<table><caption>Table 4: Summary of methods, generative models, model behaviors and hyperparameters in the largest experiments. For image generation task, we report *n* as number of candidate images while for token prediction task, *n* is taken as the number of tokens in the candidate datasets.</caption><thead><tr><th>Method</th><th>Generative Models</th><th>Model Behaviors</th><th>p</th><th>n</th><th>m</th><th>k</th></tr></thead><tbody><tr><td>TRAK</td><td>Language Models</td><td>Token Prediction</td><td>3 × 10<sup>8</sup></td><td>1.5 × 10<sup>6</sup></td><td>10</td><td>4000</td></tr><tr><td rowspan="2">JourneyTrak</td><td>CLIP</td><td>Multimodal Contrastive Learning</td><td>ResNet-50</td><td>MS COCO</td><td>100</td><td>20000</td></tr><tr><td>Diffusion Models</td><td>Image Generation</td><td>3.57 × 10<sup>7</sup></td><td>CIFAR-10</td><td>100</td><td>16384</td></tr><tr><td>D-Trak</td><td>Diffusion Models</td><td>Image Generation</td><td>3.57 × 10<sup>7</sup></td><td>5 × 10<sup>3</sup></td><td>1</td><td>32768</td></tr><tr><td>DAS</td><td>Diffusion Models</td><td>Image Generation</td><td>2.55 × 10<sup>6</sup></td><td>1.25 × 10<sup>4</sup></td><td>1</td><td>32768</td></tr><tr><td>DsDM</td><td>Language Models</td><td>Token Prediction</td><td>11.9 × 10<sup>7</sup></td><td>5 × 10<sup>3</sup></td><td>1</td><td>32768</td></tr><tr><td></td><td></td><td></td><td>1.25 × 10<sup>8</sup></td><td>6 × 10<sup>9</sup></td><td>4</td><td>16384</td></tr></tbody></table>

<table><caption>Table 5: Summary of methods and their target functions. For language model tasks, *T* refers to the context length. For D-Trak method, we specify the target function used to compute the gradient term in TRAK. $\bar{p}$ defines the mean probability that the model correctly predicts the masked token or the next token in example *x*. We also list the *Q* matrix in the implementation for TRAK-based methods. † marks terms using negative output-to-loss gradient. We summarize the implementation of *Q* matrix in Table.</caption><thead><tr><th>Method</th><th>Target Function</th><th>Formulation</th><th>Q Matrix</th></tr></thead><tbody><tr><td>TRAK</td><td>Multi-class margin</td><td>log (p(x)/1-p(x))</td><td>diag({1 - p(x)})†</td></tr><tr><td>TRAK - Masked LM</td><td>Multi-class margin</td><td>∑<sub>j∈masked tokens</sub> log (p(x<sub>j</sub>|x<sub>-j</sub>:θ)/1-p(x<sub>j</sub>|x<sub>-j</sub>:θ))</td><td>diag({1 - p̄(x)})†</td></tr><tr><td>TRAK - CLIP</td><td>Multi-class margin</td><td>-log ∑<sub>1≤j≤n</sub> exp(S<sub>ij</sub> - S<sub>ii</sub>) - log ∑<sub>1≤j≤n</sub> exp(S<sub>ji</sub> - S<sub>ii</sub>)</td><td>diag({1 - p<sub>1</sub>(x,y)+p<sub>2</sub>(x,y)})†</td></tr><tr><td>JourneyTrak</td><td>DDPM loss function</td><td>L<sub>simple</sub> = 1/k ∑<sub>i=1</sub><sup>k</sup> ||ε<sub>i</sub> - ε<sub>θ</sub> (√α<sub>t</sub> x̂<sub>t</sub><sup>(t)</sup> + √1 - α<sub>t</sub> ε<sub>i</sub>, t) ||<sub>2</sub><sup>2</sup></td><td>I</td></tr><tr><td>D-Trak</td><td>Norm-based loss function</td><td>L<sub>square</sub>(x, θ) = E<sub>t,ε</sub> [||ε<sub>θ</sub>(x<sub>t</sub>, t) ||<sub>2</sub><sup>2</sup>]</td><td>I</td></tr><tr><td></td><td></td><td>L<sub>Avg</sub>(x, θ) = E<sub>t,ε</sub> [Avg(ε<sub>θ</sub>(x<sub>t</sub>, t))]</td><td>I</td></tr><tr><td></td><td></td><td>L<sub>p-norm</sub>(x, θ) = E<sub>t,ε</sub> [||ε<sub>θ</sub>(x<sub>t</sub>, t) ||<sub>p</sub>]</td><td>I</td></tr><tr><td>DAS</td><td>Noise prediction</td><td>ε<sub>θ</sub>(x<sub>t</sub>, t)</td><td>-</td></tr><tr><td>DsDM</td><td>Multi-class margin</td><td>∑<sub>j=2</sub><sup>T</sup> log (p(x<sub>j</sub>|x<sub>-j</sub>:θ))</td><td>diag({1 - p̄(x)})†</td></tr><tr><td>TrackStar</td><td>Negative log-likelihood</td><td>-∑<sub>j=2</sub><sup>T</sup> log p(x<sub>j</sub> | x<sub>-j</sub>:θ)</td><td>I</td></tr></tbody></table>

### 2.2.5 Key Takeaways and Future Directions

Recent advances in scaling influence functions for large models have converged on a common recipe: applying gradient projection techniques (Park et al., 2023; Chang et al., 2025) to reduce dimensionality, combined with matrix factorization (Grosse et al., 2023) or other Hessian approximation techniques (Kwon et al., 2024; Młodzieżiec et al., 2025) to make computations tractable. While these modifications have enabled influence estimation in settings that were previously infeasible, rigorous ablations are lacking; for instance, it remains unclear whether FIM or GGN yields better empirical performance. More fundamentally, after these efficiency-oriented modifications, it is unclear whether the resulting metric still approximates LOO or instead tracks a different quantity. Addressing these questions represents an important direction for future work.

## 2.3 Category: Weighted Marginal Contribution Methods

Another line of research quantifies influence through the *weighted marginal contribution* of training samples. In this framework, a sample's marginal impact is evaluated over subsets of the training dataset, with each subset assigned a specific weight. Unlike simple LOO which assesses a sample's importance only in relation to the full training dataset, this weighted aggregation captures more nuanced interactions between data points, leading to a more comprehensive and robust estimation of sample influences. Notable methods in this category include Data Shapley (Ghorbani and Zou, 2019; Jia et al., 2019b), Beta Shapley (Kwon and Zou, 2022), and Data Banzhaf (Wang and Jia, 2023a), which we collectively refer to as *probabilistic values* (Dubey and Weber, 1977) (the formal probabilistic interpretation will be discussed in Equation (8)). In Section 2.3.1, we first review these key concepts and discuss their respective strengths and limitations. Next, we explore efficient methods for computing probabilistic values in Section 2.3.2. Finally, we discuss the limitations of this research direction, particularly in the context of generative AI, and outline potential future work in Section 2.3.3.
We note that methods in this category are closely related to game theory, and many works are directly motivated by the desire to establish a principled framework for compensating data holders (Ghorbani and Zou, 2019). As a result, the term *data valuation* appears more frequently than data attribution, although they share the same technical meaning. Throughout this section, we use data valuation and data attribution interchangeably.

**Notation.** Given a training set $D = \{z_i\}_{i=1}^n$, let $S \subset D$ be a subset of the training data, and let $\theta_S$ denote the model parameters obtained by training on $S$. A utility function $U: 2^D \to \mathbb{R}$ quantifies the performance of the model trained on $S$, typically via an external validation set$^{10}$ $V = \{v_i\}_{i=1}^m$. When the context is clear, we write $U(S)$ as the evaluation of $S$ and ignore $V$. We denote the marginal contribution of a data point $z_i$ to a subset $S$ as $\Delta(z_i, U, S) = U(S \cup \{z_i\}) - U(S)$. For instance, the LOO error of $z_i$ can be expressed as $\Delta(z_i, U, D \setminus \{z_i\})$. Finally, we denote a *value function* (Jia et al., 2019b) by $\phi^U$, which assigns a value $\phi_i^U$ to each $z_i$, and omit the superscript when the context is clear.

### 2.3.1 Notions of Probabilistic Values

**Data Shapley.** (Ghorbani and Zou, 2019; Jia et al., 2019b) is one of the earliest principled frameworks for data valuation (Sim et al., 2022) in machine learning. It is based on the Shapley value (Shapley, 1953), a concept originating from cooperative game theory. Formally, the Shapley value for $z_i$ is computed as the weighted average of its marginal contributions over all possible subsets $S \subseteq D \setminus \{z_i\}$

$$ \phi_i = \frac{1}{n} \sum_{k=1}^{n} \binom{n-1}{k-1}^{-1} \sum_{S \subseteq D \setminus \{z_i\}, |S|=k-1} \Delta(z_i, U, S). $$

Ghorbani and Zou (2019) argue that Data Shapley offers an *equitable* division of value, whereas LOO does not. More specifically, Data Shapley is recognized as the *unique* value function that simultaneously satisfies the following four axioms, which are typically desirable for data valuation:

1. **Null Player:** If adding $z_i$ has no effect on any subset, it should receive a value of zero. Formally, if $U(S \cup \{z_i\}) = U(S), \forall S \subseteq D \setminus \{z_i\}$, then $\phi_i = 0$.

2. **Symmetry:** If $z_i$ and $z_j$ contribute equally to all subsets, they should receive the same value. Formally, if $U(S \cup \{z_i\}) = U(S \cup \{z_j\}), \forall S \subseteq D \setminus \{z_i, z_j\}$, then $\phi_i = \phi_j$.

3. **Additivity:** The values are additive under the utility functions. Formally, for two utility functions $U_1, U_2$ and $\alpha_1, \alpha_2 \in \mathbb{R}$, $\phi_i^{\alpha_1 U_1 + \alpha_2 U_2} = \phi_i^{\alpha_1 U_1} + \phi_i^{\alpha_2 U_2}$.

4. **Efficiency:** The value of the entire dataset is completely distributed among all samples. Formally, $\sum_{i=1}^n \phi_i = U(D)$.

Since its inception, Data Shapley has been applied to a variety of machine learning tasks (Shim et al., 2021; Liu et al., 2022; Rozemberczki et al., 2022) and has motivated several generalizations, including distributional Shapley (Ghorbani et al., 2020; Kwon et al., 2021; Li and Yu, 2024a) and class-wise Shapley (Schoch et al., 2022). However, despite these promising advances, Data Shapley faces several issues in both its theoretical formulation and practical implementation. These limitations have motivated the development of alternative data valuation schemes, which we discuss below.

**Beta Shapley.** Kwon and Zou (2022) argue that the efficiency axiom is not necessary for machine learning for two main reasons. First, data valuation based on utility functions defined on some test set is typically not directly tied to monetization. Second, the ordering of samples is often sufficiently informative for many applications (e.g., data selection or mislabeled data detection), so that knowing the exact absolute values adds little extra value. Based on these arguments, they propose removing the efficiency axiom, leading to a broader class of value functions known as *semivalues* (Dubey and Weber, 1977; Dubey et al., 1981). Semivalues satisfy the other three axioms, even though they may not always adhere to the efficiency axiom. Formally,

**Definition 2.1 (Semivalue).** A value function $\phi$ is considered a semivalue if and only if there exists a set of weights $\{\alpha_k^{(n)}, k = 1, \dots, n\}$ satisfying $\sum_{k=1}^n \binom{n-1}{k-1} \alpha_k^{(n)} = 1$, such that the value function $\phi$ can be expressed as follows:

$$ \phi_i := \sum_{k=1}^{n} \alpha_k^{(n)} \sum_{S \subseteq N \setminus \{z_i\}, |S|=k-1} \Delta(z_i, U, S). \quad (8) $$

<footnote>$^{10}$The choice (e.g., quality or representativeness) of the validation set can significantly impact the resulting data values. Motivated by this, several works have addressed the challenge of data valuation in the absence of reliable validation datasets (Kwon and Zou, 2023; Lin et al., 2024c; Jahagirdar et al., 2024).</footnote>
The condition $\sum_{k=1}^{n} \binom{n-1}{k-1} \alpha_k^{(n)} = 1$ naturally induces a probability distribution over all subsets $S \subseteq N \setminus \{z_i\}$, where subsets of the same size receive equal weight. Under this distribution, $\phi_i$ represents the *expected contribution* of $z_i$ when joining a randomly sampled subset $S$. This probabilistic perspective also motivates the term *probabilistic value*. Kwon and Zou (2022) further propose a simple characterization of the coefficients: for any probability measure $\xi$ on $[0, 1]$, a set of weights $\{\alpha_k^{(n)}\}_{k \in [n]}$ defined by

$$ \alpha_k^{(n)} = \int_0^1 t^{k-1}(1-t)^{n-k} d\xi(t) \qquad (9) $$

satisfies the condition $\sum_{k=1}^{n} \binom{n-1}{k-1} \alpha_k^{(n)} = 1$, making the corresponding $\phi$ a semivalue. Mathematically, Equation (9) is equivalent to sampling each data point independently with a probability $t$, where $t$ is drawn from the probability distribution $\xi$. This formulation also arises in the Average Marginal Effect (AME) (Lin et al., 2022), which is motivated by measuring multiple treatment effects using randomized experiments in the literature of causal inference (Imbens and Rubin, 2015; Egami and Imai, 2019).

Both Data Shapley and LOO can be viewed as specific instances of a semivalue, in particular as instantiations of Equation (9). Concretely, Data Shapley is obtained by setting $\xi = \text{Unif}([0, 1])$, which leads to

$$ \alpha_k^{(n)} = \frac{1}{n} \binom{n-1}{k-1}^{-1}. $$

In contrast, LOO is obtained by setting $\xi = \delta_1$ (the Dirac measure concentrated at $t = 1$, whose cumulative distribution function is $1(t \ge 1)$), yielding

$$ \alpha_k^{(n)} = 1[k = n]. $$

Beta Shapley, on the other hand, employs the beta distribution $\xi = \text{Beta}(\alpha, \beta)$, which results in

$$ \alpha_k^{(n)} = \frac{\text{Beta}(\alpha + n - k, \beta + k - 1)}{\text{Beta}(\alpha, \beta)}. $$

We note that Beta Shapley is a more general and flexible notion that encapsulates both Data Shapley and LOO as special cases: Data Shapley is recovered by setting $(\alpha, \beta) = (1, 1)$, while LOO is obtained by setting $\alpha = 1$ and taking $\beta \to \infty$. More generally, $\alpha$ and $\beta$ can be interpreted as the weights placed on subsets of small and large cardinality, respectively. The main thesis of Kwon and Zou (2022) is that the marginal contributions calculated at sufficiently large subsets is often negligible, making it difficult to distinguish samples that significantly enhance model performance. Therefore, they recommend choosing $\alpha > \beta$, which effectively downweights noisy contributions from larger subsets and provides a more robust estimation.

**Data Banzhaf.** Wang and Jia (2023a) observe that in machine learning, the computation of utility scores $U(S)$ can be noisy due to the inherent randomness in weight initialization and minibatch sampling. This further makes the calculated values $\phi_i$'s unreliable; as a matter of fact, different runs of the same data valuation algorithm may produce inconsistent value rankings. To quantify the brittleness, they introduce the concept of *safety margin*, defined as the maximum perturbation in utility scores that can be tolerated while preserving the relative value ordering of every pair of data points. They then propose the Data Banzhaf framework based on the Banzhaf value (Banzhaf III, 1964), a classic notion from cooperative game theory and an instantiation of a semivalue. Formally, it corresponds to setting $\xi = \delta_{0.5}$ in Equation (9), resulting in

$$ \alpha_k^{(n)} = \frac{1}{2^n} $$

for all $k$. In other words, Data Banzhaf assigns *uniform* weight to all subsets, regardless of their sizes. Wang and Jia (2023a) further demonstrate that Data Banzhaf achieves the largest safety margin among all semivalues, enabling it to *robustly* distinguish data quality in the presence of noise.

**Weighted Banzhaf values.** Li and Yu (2023) demonstrate that if one assumes a specific *structure* for the stochasticity (rather than relying on the safety margin, which is a *worst-case* robustness notion), Data Banzhaf may not be the most robust semivalue. Concretely, they show that when the stochasticity is parameterized by Kronecker noise, the uniquely most robust semivalue which minimizes the worst-case entropy belongs to the family of weighted Banzhaf values (*w*-weighted Banzhaf value), a generalized form of the Banzhaf value. Formally, this corresponds to setting $\xi = \delta_w$ in Equation (9), yielding

$$ \alpha_k^{(n)} = w^{k-1}(1-w)^{n-k}. $$

Notably, 0.5-weighted Banzhaf value exactly recovers the standard Banzhaf value used in Data Banzhaf.
## 2.3.2 Efficient Computations of Probabilistic Values

Having reviewed the definitions of various probabilistic values in Section 2.3.1, we now turn our attention to their computation. Taking Shapley values as an example, computing them *exactly* requires evaluating $U$ on $2^n$ different subsets of the training data, each corresponding to training a different model. This exponential complexity renders an exact calculation infeasible even for moderately sized datasets. In fact, Deng and Papadimitriou (1994) show that computing these values is NP-hard. In what follows, we discuss two main directions to efficiently compute probabilistic values: reducing the number of evaluations required and lowering the computational cost per evaluation.<fnref n="11" />

<table><caption>Table 6: Summary of applicability and convergence rates for various methods of estimating probabilistic values. Theoretically, OFA (Li and Yu, 2024b) achieves the best-known convergence rate for all probabilistic values on average. Empirically, Li and Yu (2024b) report that OFA generally outperforms other methods, although the complement estimator and MSR estimator perform comparably to OFA for the Shapley and Banzhaf values, respectively. For practitioners, since the target probabilistic value is typically known in advance, we recommend using the complement and MSR estimators for their simpler implementation.</caption><thead><tr><th>Method</th><th>Applicability</th><th>Convergence rate</th><th>Remark</th></tr></thead><tbody><tr><td>Permutation Sampling (Castro et al., 2009)</td><td>Shapley</td><td>O(n<sup>2</sup> log n)</td><td>—</td></tr><tr><td>Group Testing (Jia et al., 2019b)</td><td>Shapley</td><td>O(n(log n)<sup>2</sup>)</td><td>Empirically underperforms</td></tr><tr><td>Weighted Sampling Lift (Kwon and Zou, 2022)</td><td>All</td><td>O(n<sup>2</sup> log n)</td><td>—</td></tr><tr><td>Sparse Regression (Lin et al., 2022)</td><td>Partial</td><td>O(k log n)</td><td>Does not apply to Shapley value</td></tr><tr><td>Complement Estimator (Zhang et al., 2023a)</td><td>Shapley</td><td>Unknown</td><td>—</td></tr><tr><td>MSR (Wang and Jia, 2023a)</td><td>Banzhaf</td><td>O(n log n)</td><td>—</td></tr><tr><td>OFA (Li and Yu, 2024b)</td><td>All</td><td>O(n log n)</td><td>—</td></tr></tbody></table>

**Reducing the number of evaluations.** Since the definitions of probabilistic values take the form of a weighted sum, a natural strategy is to employ Monte Carlo sampling to approximate them, which inherently ensures *unbiasedness*. The central question for these approximation methods is: *how many evaluations are required to ensure that the estimated probabilistic values are close to the ground truth?* Here, closeness is typically characterized via the $(\epsilon, \delta)$-approximation:

$$ \Pr(\|\hat{\phi} - \phi\|_2 \le \epsilon) \ge 1 - \delta, $$

where $\hat{\phi}$ denotes the estimated score vector (estimator in short). We refer to the number of evaluations required (denoted by $T$) as a function of $n$ as the convergence rate for each method, which will be the main quantity of interest. In what follows, we review several representative approximation methods, discuss their applicability and convergence rates, and summarize them in Table 6. Before that, we make an additional comment.

**Remark 2.2.** Approximating probabilistic values, particularly the Shapley value, is a research problem that extends beyond data attribution. For instance, it has been extensively studied in the context of *feature attribution* (Lundberg and Lee, 2017), and many of those methods directly transfer to data attribution. Nevertheless, to avoid confusion, we focus our discussion on works that explicitly address data attribution, and refer readers to Li and Yu (2024a,b) for a more comprehensive discussion, which includes methods such as *KernelSHAP* (Lundberg and Lee, 2017; Covert and Lee, 2021), *FastSHAP* (Jethani et al., 2022), *SHAP-IQ* (Fumagalli et al., 2024), *ARM* (Kolpaczki et al., 2024), and others.

**Permutation sampling.** One of the most widely used estimator for Shapley value is based on permutation sampling, proposed by Castro et al. (2009). The key observation is that the Shapley value can be equivalently formulated as

$$ \phi_i = \frac{1}{n!} \sum_{\pi \in \Pi} \Delta(z_i, U, \mathcal{P}^i(\pi)), $$

where $\Pi$ contains all permutations of $[n]$ and $\mathcal{P}^i(\pi)$ is the subset that contains all players preceding $i$ in $\pi$. Based on this equation, the permutation estimator is naturally defined as

$$ \hat{\phi}_i = \frac{1}{T} \sum_{j=1}^{T} \Delta(z_i, U, \mathcal{P}^i(\pi_j)), $$

with $\{\pi_j\}_{j=1}^T$ being uniformly sampled (with replacement) from $\Pi$. Assume that $U$ is bounded, an application of Hoeffding's bound suggests that $O(n \log n)$ samples are required to ensure that $\hat{\phi}_i$ is close to $\phi_i$, implying a total requirement of $T = O(n^2 \log n)$ evaluations.

<footnote n="11">¹¹We include the cost of model training in evaluation.</footnote>
Permutation sampling is less efficient because the marginal contribution $\Delta(z_i, U, S_j)$ (which requires two evaluations) is only used to update $\hat{\phi}_i$, thereby introducing an extra factor of $n$ in the convergence rate. Nevertheless, due to its simplicity, permutation sampling has served as a cornerstone for many subsequent works. For example, Ghorbani and Zou (2019) propose Truncated Monte Carlo Shapley (TMC-Shapley), which truncates the calculation of marginal contributions in a sampled permutation once the utility evaluation is close to that of the full training set. While TMC-Shapley is more efficient than the standard permutation sampling method, it sacrifices unbiasedness. Jia et al. (2019b) leverage permutation sampling to construct the label vector within a compressed sensing framework, achieving a convergence rate of $O(n \log \log n)$. It is important to note that this method relies on the assumption that the ground truth Shapley values are “approximately sparse” (i.e., most $\phi_i$’s are close to their mean).

*Group testing.* Jia et al. (2019b) propose an estimator based on group testing (Du and Hwang, 1999). The key idea is to estimate the Shapley differences between all data pairs instead of the Shapley values directly

$$ \phi_i - \phi_j = \frac{1}{n-1} \sum_{S \subseteq D \setminus \{z_i, z_j\}} \frac{1}{\binom{n-2}{|S|}} [U(S \cup \{z_i\}) - U(S \cup \{z_j\})], $$

followed by solving a feasibility problem with an additional constraint imposed by the efficiency axiom to obtain the $\hat{\phi}_i$'s. This approach achieves a convergence rate of $O(n(\log n)^2)$, saving a factor of $n$ compared to the permutation sampling approach, thanks to increased *sample reuse*: each utility evaluation contributes to the estimation of the Shapley values for all samples. Wang and Jia (2023b) provide a refined analysis of the group testing method and also propose an improved estimator using the dummy player technique (we refer interested readers to their work for details). On the other hand, they make two key observations: 1) group testing does not achieve *maximum sample reuse* (MSR), a concept that will be discussed in detail shortly; 2) empirically, group testing does not significantly outperform permutation sampling, despite having a faster convergence rate.

*(Weighted) sampling lift.* Equation (8) provides a natural probabilistic interpretation. First, the outer summation corresponds to a discrete probability distribution over the *cardinality* of the subset, with the probability of a subset having cardinality $k$ given by $\binom{n-1}{k-1}\alpha_k^{(n)}$. Second, the inner summation represents a uniform distribution over all subsets of cardinality $k$. This idea underlies sampling lift (Moehle et al., 2022), where the estimator is defined as

$$ \hat{\phi}_i = \frac{1}{T} \sum_{j=1}^{T} \Delta(z_i, U, S_j), $$

and the $\{S_j\}_{j=1}^T$ are sampled independently from the two-step process above: first, draw the subset size from the discrete probability distribution, then uniformly sample a subset of that size.

Kwon and Zou (2022) extend this idea with a weighted version of sampling lift. Specifically, they reformulate Equation (8) as

$$ \phi_i = \frac{1}{n} \sum_{k=1}^{n} n \alpha_k^{(n)} \binom{n-1}{k-1} \left[ \frac{1}{\binom{n-1}{k-1}} \sum_{\substack{S \subseteq D \setminus \{z_i\} \\ |S| = k-1}} \Delta(z_i, U, S) \right], $$

and propose the estimator

$$ \hat{\phi}_i = \frac{1}{T} \sum_{j=1}^{T} n \alpha_{|S_j|}^{(n)} \binom{n-1}{|S_j| - 1} \Delta(z_i, U, S_j), $$

where the sampling process for $\{S_j\}_{j=1}^T$ is modified so that, in the first step, the cardinality of $S_j$ is sampled uniformly from $[n]$. For the Shapley value, the coefficients $n\alpha_{|S_j|}^{(n)} \binom{n-1}{|S_j|-1}$ simplify to 1. However, for general probabilistic values, computing these coefficients is numerically unstable for large $n$ due to the combinatorial factors.

*(Weighted) sampling lift* achieves a convergence rate of $T = O(n^2 \log n)$, with the proof strategy mirroring that of the permutation sampling approach (e.g., see Proposition 5 in Li and Yu (2024a)). Importantly, since the estimator is directly derived from Equation (8), it can be applied to any probabilistic values.

*Sparse regression.* Lin et al. (2022) observe that a subfamily of probabilistic values can be formulated as the optimal solution to a least squares problem. This methodology also underlines earlier works in feature attribution such as KernelSHAP (Lundberg and Lee, 2017) and more recent works such as GELS (Li and Yu, 2024a). Specifically, when $\xi$ (as defined in Equation (9)) satisfies $C_\xi =: \int_0^1 \frac{1}{w(1-w)} d\xi(w) < \infty$, we have

$$ \phi = \arg \min_{v \in \mathbb{R}^N} \mathbb{E}_S \left[ (Y(S) - X(S)^\top v)^2 \right]. $$
Here, $Y(S) = U(S)$ and $X(S)_i$ (the $i$-th coordinate of $X(S)$) equals $\frac{1}{C_\xi w}$ if $i \in S$ and otherwise $-\frac{1}{C_\xi (1-w)}$. The random subset $S$ is generated by first sampling $w \in [0, 1]$ from $\xi$, and then including each sample independently with probability $w$. The final estimator is given by

$$ \hat{\phi}_i = \arg\min_{v \in \mathbb{R}^N} \frac{1}{T} \sum_{j=1}^{T} (Y(S_j) - X(S_j)^\top v)^2 + \lambda \|v\|_1, $$

where an $L_1$-regularization term is added to perform sparse regression (i.e., LASSO (Tibshirani, 1996)). Assuming $\|\phi\|_0 = k$, Lin et al. (2022) show that the convergence rate of this estimator is $O(k \log n)$, a significant improvement over previous methods. However, the requirement that $C_\xi$ be finite precludes the application of this method to the Shapley value: in particular, when $\xi = \text{Unif}([0, 1])$, we have $\int_0^1 \frac{1}{w(1-w)} dw = \int_0^1 (\frac{1}{w} + \frac{1}{1-w}) dw = \infty$.

**Complementary contribution.** Another equivalent formulation of the Shapley value is based on the complementary formula (Harsanyi, 1982). Instead of expressing $\phi_i$ as a weighted summation of marginal contributions, it is written as a weighted summation of complementary contributions (CC), which are formally defined as $CC(S) = U(S) - U(D \setminus S)$ for a subset $S$. The complementary formula states that

$$ \phi_i = \frac{1}{n} \sum_{S \subseteq D \setminus \{z_i\}} \left( \frac{n-1}{|S|} \right)^{-1} CC(S \cup \{z_i\}). $$

Based on this, Zhang et al. (2023a) propose the following complement estimator:

$$ \hat{\phi}_i = \frac{1}{n} \sum_{k=1}^{n} \frac{1}{T_{i,k}} \sum_{j=1}^{T} CC(S_j) (1(i \in S_j, |S_j| = k) - 1(i \notin S_j, n - |S_j| = k)), $$

where $T_{i,k} = \sum_{j=1}^T (1(i \in S_j, |S_j| = k) + 1(i \notin S_j, n - |S_j| = k))$, and $\{S_j\}_{j=1}^T$ are sampled independently by first sampling the subset size uniformly from $[n]$, and then uniformly drawing a subset of that size. The main advantage of this estimator is that each sample appears either in a subset $S$ or in its complement, so a single utility evaluation can be leveraged to update the estimates for all samples. In fact, as suggested by Li and Yu (2024b), the complement estimator adheres to the principle of MSR, which will be discussed in the next method. Consequently, although there is no formal convergence guarantee, it is reported to have the best empirical performance for estimating the Shapley value (Li and Yu, 2024b).

**Maximum sample reuse.** The principle of MSR has appeared in prior works, such as group testing and the complement estimator, and was formally introduced in Wang and Jia (2023a) to efficiently approximate the Banzhaf value. Formally, MSR requires that each evaluation is used to update the estimates for all $\phi_i$'s. To achieve this, Wang and Jia (2023a) rewrite the Banzhaf value as

$$ \phi_i = \mathbb{E}_{S|z_i \in S}[U(S)] - \mathbb{E}_{S|z_i \notin S}[U(S)], \quad (10) $$

where $S$ is uniformly distributed over all subsets of $D$. They then propose the following MSR estimator:

$$ \hat{\phi}_i = \frac{1}{T_i^+} \sum_{j=1}^{T} U(S_j) 1(i \in S_j) - \frac{1}{T_i^-} \sum_{j=1}^{T} U(S_j) 1(i \notin S_j), $$

where $T_i^+ = \sum_{j=1}^T 1(i \in S_j)$ and $T_i^- = \sum_{j=1}^T 1(i \notin S_j)$, and $\{S_j\}_{j=1}^T$ are sampled independently from the uniform distribution over all subsets of $D$. Thanks to MSR, this estimator achieves a convergence rate of $O(n \log n)$.

While this estimator can naturally be adapted to the weighted Banzhaf value, it cannot be applied to the Shapley value because it is impossible to construct a suitable distribution over $D$ to express Shapley values in the form of Equation (10) (Wang and Jia, 2023a; Li and Yu, 2024a). We note that this is a limitation of marginal contributions: as demonstrated in the previous method, the complement estimator achieves MSR by using complementary contributions instead.

**One-sample-fits-all.** Most of the estimators discussed thus far, with the exception of (weighted) sampling lift, cannot be applied to all probabilistic values. However, weighted sampling lift not only has a suboptimal convergence rate $O(n^2 \log n)$, but its convergence can be further degraded by the coefficients $n\alpha_{|S_j|}^{(n)}$ (referred to as amplifying scalars). To address this challenge, Li and Yu (2024b) introduce the One-sample-Fits-All (OFA) framework, which simultaneously satisfies the following desirable properties: 1) it requires sampling subsets only once to approximate
all probabilistic values; 2) it adheres to the principle of MSR; and 3) it avoids the amplifying scalars inherent in the weighted sampling lift estimator. The framework is based on the following formulation of probabilistic values:

$$ \phi_i = \sum_{k=1}^{n} \binom{n-1}{k-1} \alpha_k^{(n)} \cdot (\mathbb{E}_{z_i \in S, |S|=k} [U(S)] - \mathbb{E}_{z_i \notin S, |S|=k-1} [U(S)]), $$

where each expectation is taken with respect to the corresponding uniform distribution. The two expectations, denoted as $\hat{\phi}_{i,s}^{+} = \mathbb{E}_{z_i \in S, |S|=k} [U(S)]$ and $\hat{\phi}_{i,s}^{-} = \mathbb{E}_{z_i \notin S, |S|=k-1} [U(S)]$, are estimated using the two-step process of the sampling lift estimator. Furthermore, by optimizing the discrete probability distribution over the cardinality of subsets, OFA achieves a convergence rate of $O(n \log n)$ for all probabilistic values on average.

**Lowering the cost per evaluation.** Despite significant progress along the previous direction, the minimum requirement of utility evaluation remains at least linear in the number of training samples, which is generally impractical for training large-scale models from scratch. In the following, we discuss a complementary approach that aims to address this challenge by lowering the computational cost per evaluation. This reduction can be achieved either by employing techniques to decrease the training cost directly, or by substituting the entire training process with a lighter, more computationally efficient proxy model<fnref n="12" />.

**Reducing training cost.** Building upon the framework of permutation sampling, Ghorbani and Zou (2019) introduce Gradient Shapley (G-Shapley), which approximates a fully trained model by training for a single epoch with a batch size of 1 (this is referred to as an incrementally trainable model in Jia et al. (2019b)) and then computing the marginal contributions as data is sequentially added. Lu et al. (2025) propose Sparsified Fine-Tuning Shapley (Sparsified-FT Shapley) for diffusion models. This method involves pruning a pre-trained diffusion model and fine-tuning it on the full training set to produce a sparsified yet performant model, which is then used to fine-tune on subsets of the training data for utility evaluation. Both approaches significantly reduce the computational cost compared to full retraining.

**KNN as a proxy.** Jia et al. (2019a) derive a closed-form expression for the Shapley value in the context of the K-Nearest Neighbors (KNN) classifier. Specifically, they introduce a KNN utility function that quantifies the impact of a subset S on KNN classification accuracy:

$$ U(S) = \frac{1}{K} \sum_{k=1}^{\min\{K, |S|\}} 1[y_{\alpha_k(S)} = y_{\text{val}}], $$

where $(x_{\text{val}}, y_{\text{val}})$ represents a validation sample and $\alpha_k(S)$ denotes the indices of the training data in S sorted by increasing distance to $x_{\text{val}}$. Building on this utility function, they derive a recursive algorithm that computes the exact Shapley values in $O(n \log n)$ time:

$$ \phi_{\alpha_n} = \frac{1}{n} 1[y_{\alpha_n} = y_{\text{val}}], \\ \phi_{\alpha_i} = \phi_{\alpha_{i+1}} + \frac{1[y_{\alpha_i} = y_{\text{val}}] - 1[y_{\alpha_{i+1}} = y_{\text{val}}]}{K} \min\{K, i\}, i = n-1, \dots, 1. $$

More importantly, Jia et al. (2019a) observe that this framework for calculating Shapley values for KNN classifiers can be extended to deep neural networks. Specifically, they propose to train a KNN classifier based on the embeddings (e.g., extracted from deeper layers of the network or pre-trained feature extractors) of each $x_i$ along with their corresponding labels $y_i$, and then perform data attribution in the latent space. This effectively reduces the problem of data attribution for deep neural networks to that of KNN classification. This approach, which avoids the need to repeatedly train complex models, has been recognized in follow-up works as one of the most practical and scalable techniques for large-scale data valuation (Jia et al., 2021; Pandl et al., 2021; Karlaš et al., 2024).

**NTK as a proxy.** The neural tangent kernel (NTK) (Jacot et al., 2018) is a powerful analytical tool in deep learning theory, used to investigate the optimization and generalization of (deep) neural networks. It posits that with a suitable scale of random initialization, the training of a sufficiently wide neural network remains close to its initialization, effectively reducing the learning dynamics to a kernel regression characterized by NTK. In the context of probabilistic values, NTK enables an approximate evaluation of the utility function $U(S)$ without the need for full training, thereby significantly lowering computational costs. For instance, Wu et al. (2022b) introduce data valuation at initialization (DaVinz), leveraging an off-the-shelf generalization bound from Arora et al. (2019) to compute $U(S)$. Similarly, Wang et al. (2024g) propose fine-tuning-free Shapley value (FreeShap), which employs the empirical NTK (Malladi et al., 2023) to characterize $U(S)$. We comment that NTK is particular useful in the fine-tuning regime where the model is expected to remain close to the pretrained checkpoint.

<footnote n="12">¹²We focus on general-purpose proxy models and omit those with a relatively narrow scope, such as Data-OOB (Kwon and Zou, 2023), which is only applicable to bagging models.</footnote>
### 2.3.3 Key Takeaways and Future Directions

The weighted marginal contribution method builds upon the widely recognized Shapley value framework, with Data Shapley remaining the most popular notion of probabilistic value for training data attribution. Despite significant progress in this area, this approach is generally less scalable than influence functions (as discussed in Section 2.2) because it requires at least a linear number of model trainings relative to the dataset size. In fact, most algorithms discussed in this section, except those based on training-free proxy models, are not directly applicable in the context of generative AI. Nonetheless, its strong theoretical foundations, especially its satisfaction of the four desirable axioms, make it a compelling alternative in scenarios where attribution is needed at a higher level, such as domains (Kang et al., 2024; Liu et al., 2025) or copyright owners (Wang et al., 2024c). This is particularly relevant in emerging applications within the data economy (Sestino et al., 2025) and data markets (Zhang et al., 2024). We envision future work focusing on more efficient approximations of Data Shapley and related concepts tailored to these applications.

**Remark 2.3 (Preservation of axioms under Monte-Carlo approximation).** Among the four classical axioms that characterize the Shapley value, the symmetry axiom is often viewed as the key fairness requirement: equally useful points should receive equal rewards. When the exact Shapley value is approximated with Monte-Carlo (MC) sampling, however, MC estimates satisfy Symmetry only in expectation. With a finite sample budget, two exchangeable points typically obtain slightly different estimates because the Monte-Carlo noise that perturbs their marginal contributions is independent. There is an emerging line of work developing Shapley value approximation methods that provably preserve most of the Shapley axioms. Wang et al. (2024d) developed a deterministic approximation algorithm for weighted KNN-Shapley that preserves the null player and symmetry axioms exactly.

## 2.4 Category: Training Dynamics Methods

This category of methods leverages information collected throughout the training process to assess data influence, an approach we refer to as training dynamics data attribution. The literature in this area generally divides into two groups: tracing gradient descent (TracIn) and trajectory-specific Leave-One-Out (TSL00).

1. Tracing gradient descent. Initiated by Pruthi et al. (2020), this approach, referred to as TracIn, estimates the influence of training data by approximating the cumulative change in model outputs during training.

2. Trajectory-specific LOO. Originally introduced by Hara et al. (2019) under the name SGD-influence (or counterfactual SGD), this method was later reframed as TSL00 in a follow-up study (Wang et al., 2025b) to more accurately reflect its goal of dynamically approximating leave-one-out. Specifically, it quantifies the influence of removing a single training sample on both the model parameters and outputs over the course of training by building a recursive relationship between consecutive iterations.

**Notation.** The loss incurred on a single data point is denoted by $\ell(z; \theta)$, where $\theta \in \mathbb{R}^p$ represents the model parameters and $z \in D$, given a training set $D = \{z_i\}_{i=1}^n$. We also define $f(z_{test}; \theta)$ as the target function and $z_{test}$ as a test sample.

### 2.4.1 Tracing Gradient Descent

**Derivation of tracing gradient descent.** We review the basic derivation of tracing gradient descent method family, particularly highlight the connection between TracIn-Ideal and its first-order approximation TracIn.

**TracIn-Ideal.** The high-level idea of TracIn-Ideal is accumulating the change in model outputs caused by optimizing on a training sample throughout the training process. The idealized influence of a particular training sample $z$ on a given test sample $z_{test}$ is defined as the total model output $f(z_{test}; \theta)$ change induced by the whole training process whenever $z$ is trained on. The definition can be formalized as:

$$ \text{TracIn-Ideal}(z, z_{test}) \triangleq \sum_{t: z_t = z} (f(z_{test}; \theta_t) - f(z_{test}; \theta_{t+1})), \quad (11) $$

where $t: z_t = z$ indicates the iteration where $z$ is taken as the training sample. Equation 11 is for vanilla SGD where parameters optimization is done on one training sample at a time. This ideal definition (Pruthi et al., 2020) is the goal to be approximated by tracing gradient descent method family.

**TracIn.** To make the tracing computationally tractable, TracIn uses first-order Taylor expansion to approximate the idealized influence, and converts the Equation 11 into the following inner-product formulation. The formula is also adapted to mini-batch training.
$$ \mathrm{TracIn}(z, z_{\mathrm{test}}) \triangleq \frac{1}{|B_t|} \sum_{t:z \in B_t} \eta_t \langle \nabla_{\theta} f(z_{\mathrm{test}}; \theta_t), \nabla_{\theta} \ell(z; \theta_t) \rangle, \quad (12) $$

It assumes that parameters are updated by an SGD optimization process with training batch $B_t$ and learning rate $\eta_t$ at iteration $t$. Most of the developments on tracing gradient decent family are based on this first-order approximation.

**Challenges.** Among recent studies, tracing gradient decent family can be primarily categorized into two groups: offline tracing algorithms and online tracing algorithms. Offline tracing algorithms (Pruthi et al., 2020; Lin et al., 2024b; Xie et al., 2024a; Yeh et al., 2022) present the traditional implementations to TracIn, calculating attribution results after model training, typically leveraging some cached training checkpoints. Online tracing algorithms (Wang et al., 2025a, 2024e) compute the attribution results of each training iteration during each training iteration alongside model updates, utilizing access to richer training dynamics such as backpropagation information.

The scalability of Equation 12 is the most critical challenge for tracing gradient descent methods in the generative AI era. Though several studies (Pruthi et al., 2020; Lin et al., 2024b) propose various techniques to reduce the computational complexity and achieve high usability for traditional machine learning models, the core requirement to instantiate and cache gradients for every training sample remains. This creates critical limitations when applying offline algorithms on generative AI models, which have exponentially growing parameter sizes and training datasets. Wang et al. (2025a) proposes ghost dot-product, a technique that calculates the exact inner product of a pair of data samples' gradients without instantiating them. The technique could also be leveraged for offline algorithms as well and serves as the current best practice.

Another challenge is the underperformance stemming from errors introduced by inherent assumptions and approximations. Although different resolutions (Xie et al., 2024a; Yeh et al., 2022) aim to improve gradient tracing performance, the deep, fundamental roots of these algorithms and their impact on downstream application performance require further investigation.

**Offline tracing algorithms.** Offline tracing algorithms are performed during model training. Since computing Equation 12 across every training iteration is computationally expensive, several techniques have been proposed to reduce the cost.

**TracIn-CP.** A straightforward solution is to only select a few checkpoints during the training process to approximate Equation 12. Pruthi et al. (2020) propose to further approximate the accumulation of these inner products using C fixed model checkpoints {$\theta_{t_c}$}$_{c=1}^C$ that are sparsely saved across the whole training process, where $C \ll T$. This leads to the following practical formulation of TracIn-CP:

$$ \mathrm{TracIn-CP}(z, z_{\mathrm{test}}) \triangleq \sum_{c=1}^{C} \eta_{t_c} \langle \nabla_{\theta} f(z_{\mathrm{test}}; \theta_{t_c}), \nabla_{\theta} \ell(z; \theta_{t_c}) \rangle, \quad (13) $$

where $\eta_{t_c}$ is the learning rate at $t_c$. Grad-Dot is a special case when we extend TracIn-CP to use only one checkpoint, i.e., $C = 1$. It coincidences with the similarity in parameterized family defined in Charpiat et al. (2019).

**Gradient vector compression.** Another group of popular techniques are various ways of gradient vector compression, such as layer cherry picking and random projection. Assuming deep neural networks are usually overparameterized, Pruthi et al. (2020) propose to cherry-pick certain layers (e.g., only the last layer) of a deep neural network for gradient calculation to reduce the dimensionality of gradients. Random projection maps the parameter space to a lower-dimension subspace. Both of these techniques are also widely used in Influence Function (Section 2.2.2), thus we will not describe them in detail here.

**Mitigating attribution underperformance.** Previous efficient techniques applied to Equation 12, including checkpoint selection and gradient vector compression, may degrade the performance of offline algorithms across applications. To address this, recent studies identify and try to mitigate some specific failure modes. Yeh et al. (2022) discovers the cancellation effect, where last-layer weights make the attribution result of each sample have a large magnitude that contradicts each other. The authors propose picking the first layer (e.g., the embedding layer of a language model) for the tracing attribution algorithm. Separately, Xie et al. (2024a) finds that in diffusion models (Ho et al., 2020b), training examples with large diffusion timesteps (closer to noise) exhibit higher loss gradient norms, skewing the attribution results. They propose normalizing both training and test gradient terms by their respective norms.

The ease of implementation positions approaches centered around TracIn-CP as a go-to method for various applications (Akyürek et al., 2022). Per-example gradients are calculated (and optionally cached) for each sample at a maximum available batch size $|B|$ constrained by the memory limit $M_{\mathrm{limit}}$. Nevertheless, the per-sample gradient
instantiation is memory-consuming for large generative AI, which limits the maximum available batch size to $\frac{M_{\text{limit}}}{|B|p}$, where $p$ is the parameter size.

**Online tracing algorithms.** Online tracing algorithms aim to estimate the data influence during the model training process. By leveraging available training dynamics (e.g., backpropagation information), Wang et al. (2025a) enables efficient computation of gradient inner products between data sample pairs.

*Ghost dot-product.* For a typical neural net layer without parameter sharing, the per-sample gradient w.r.t. the parameters at this layer can be computed using the input to the layer and the gradient of the loss w.r.t. the output, both of which are available during backpropagation in batch. Furthermore, the inner product of two large gradients can be computed in a layer-by-layer fashion, avoiding the need to instantiate large gradient vectors. The formula for a linear layer can be presented as

$$ \frac{\partial \ell^{(1)}}{\partial \mathbf{W}} \odot \frac{\partial \ell^{(2)}}{\partial \mathbf{W}} = \left( \mathbf{a}^{(1)} \otimes \frac{\partial \ell^{(1)}}{\partial \mathbf{s}^{(1)}} \right) \odot \left( \mathbf{a}^{(2)} \otimes \frac{\partial \ell^{(2)}}{\partial \mathbf{s}^{(2)}} \right) = \left( (\mathbf{a}^{(1)})^{\top} \mathbf{a}^{(2)} \right) \left( \left( \frac{\partial \ell^{(1)}}{\partial \mathbf{s}^{(1)}} \right)^{\top} \left( \frac{\partial \ell^{(2)}}{\partial \mathbf{s}^{(2)}} \right) \right), \quad (14) $$

where $\ell^{(1)}$ and $\ell^{(2)}$ are the losses of the pair of data samples (typically a training sample and a test sample), $\mathbf{W}$ is the weight of the linear layer, $\mathbf{a}^{(1)}$ and $\mathbf{a}^{(2)}$ are the input of the linear layer, $\mathbf{s}^{(1)}$ and $\mathbf{s}^{(2)}$ are the pre-activate output of the linear layer, $\odot$ stands for the inner product and $\otimes$ stands for the outer product.

For online algorithms, $\mathbf{a}^{(1)}$, $\mathbf{a}^{(2)}$, $\frac{\partial \ell^{(1)}}{\partial \mathbf{s}^{(1)}}$ and $\frac{\partial \ell^{(2)}}{\partial \mathbf{s}^{(2)}}$ are ready during back-propagation, so that the overhead caused by ghost dot-product algorithm is relatively small. Nevertheless, the test samples need to be prepared **before** model training for online algorithms, which may be invalid for some use cases such as attributing a generated sample (Deng et al., 2024b). Ghost dot-product can also be applied to offline algorithms and intermediate checkpoints as well, which significantly reduces the space complexity. Detailed analyses are stated in the next paragraph.

**Time and space complexity.** We tabulate the time and space complexity of the methods in Table 7.

For complexity analysis, we split the complexity analysis to “One-time” and “Amortized”. “One-time” stands for the complexity of calculating the data attribution result for one test samples on the full training dataset. “Amortized” stands for the amortized complexity of attributing continually incoming test samples and cache the gradient information if applicable.

We denote the size of training set as $n$, the parameter size as $p$, the number of training iteration as $T$, the number of checkpoints as $C$, the reduced parameter size (e.g., through random projection) as $\hat{p}$, and number of layer as $L$. In most cases $C \ll T$ and $L \ll p$. We also assume one forward pass of the target model costs $F_\theta$, and one backward pass approximately costs two forward passes.

<table><caption>Table 7: Complexity analysis of TracIn and its variants. The complexity analysis here assumes non-sequential data and the attribution for the whole training process.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">One-time</th><th colspan="2">Amortized</th></tr><tr><th>Time</th><th>Space</th><th>Time</th><th>Space</th></tr></thead><tbody><tr><td>TracIn</td><td>O(nT(F<sub>θ</sub> + p))</td><td>O(p)</td><td>O(nTp)</td><td>O(nTp)</td></tr><tr><td>TracIn-CP</td><td>O(nC(F<sub>θ</sub> + p))</td><td>O(p)</td><td>O(nCp)</td><td>O(nCp)</td></tr><tr><td>TracIn-CP + Gradient vector compression</td><td>O(nC(F<sub>θ</sub> + p̂))</td><td>O(p)</td><td>O(nCp̂)</td><td>O(nCp̂)</td></tr><tr><td>TracIn-CP + Ghost dot-product</td><td>O(nC(F<sub>θ</sub> + √pL))</td><td>O(√pL)</td><td>O(nC(√pL))</td><td>O(nC√pL)</td></tr></tbody></table>

1. TracIn and TracIn-CP perform Equation 12 and 13 by directly instantiating per-sample gradients and caching them for multiple test samples. The only difference is the constant scaler $C$ and $T$.

2. Gradient vector compression perform gradient compression right before the caching. Notably, this won't change the space complexity when calculating the gradient since mapping the gradient to a lower dimension happens after the instantiation.

3. Ghost dot-product could avoid explicitly instantiating per-sample gradient and reduce the space cost significantly, which is critical for large generative models. For each layer, caching $\mathbf{a}$ and $\frac{\partial \ell}{\partial \mathbf{s}}$ for one training samples takes $O(d_{in} + d_{out})$ space, which is proportional to the **square root of the original cost** of explicitly caching the per-sample gradient ($O(d_{in} \times d_{out})$).
**Relationship with influence function.** While both TracIn and Influence Functions (IF) are commonly used to quantify the influence of a training sample on a model's prediction for a test sample and both rely on relatively similar closed-form expressions. They differ in how they define influence. TracIn estimates the total change in the model's output during training process caused by a training sample (Equation 11). In contrast, IF approximates the change in the final model parameters and corresponding output when a training sample is removed from the training dataset (Equation 2). Another difference is that IF is agnostic to the optimizer used in the training process, while TracIn is strongly dependent on the optimizer. Stochastic gradient descent(SGD) is assumed to be the optimizer in the derivation of Equation 12, while it's not commonly used for most generative AI applications.

## 2.4.2 Trajectory-Specific LOO

The second type of dynamic method is trajectory-specific LOO (TSL00), which was first proposed by Hara et al. (2019) in the name of SGD-influence. Like influence functions, TSL00 quantifies samples influence by analyzing the counterfactual scenario of LOO. However, unlike influence functions that only consider the counterfactual at the final model checkpoint, TSL00 examines the counterfactual at all intermediate checkpoints. In particular, it uses a first-order approximation to establish a recursive relationship between the LOO influences at consecutive iterations, thereby propagating a sample's influence from its initial appearance to the final model. Importantly, this approach does not assume that the loss function is convex or that the model is trained to convergence, both of which are key assumptions in deriving influence functions, thereby offering broader applicability<fnref n="13" />. Moreover, as we will see shortly, TSL00 eliminates the need to compute the inverse-Hessian-Vector Product (iHVP), which is the most computationally expensive component in influence functions. Instead, it requires only the computation of Hessian-vector products, which can be performed much more efficiently.

**Notation.** We denote the loss function as $\ell(z; \theta)$, where $z$ represents a training sample and $\theta$ denotes the model parameters. We assume $\ell$ is twice-differentiable. The model is trained using SGD for $T$ iterations. At iteration $t$, the model parameters are represented by $\theta^{[t]}$, while $\theta_{-j}^{[t]}$ denotes the model parameters in the counterfactual scenario where sample $j$ is removed. The mini-batch used at this iteration is $S_t$, and the Hessian over the batch is given by $H^{[t]} = \frac{1}{|S_t|} \sum_{i \in S_t} \nabla_{\theta}^2 \ell(z_i; \theta^{[t]})$. Finally, we define $u$ as a query vector, which is typically chosen as the gradient of the loss of a test sample $z_{\text{test}}$ at the final iteration, i.e., $u = \nabla \ell(z_{\text{test}}; \theta^{[T]}).$

**Original formulation of TSL00.** We review the core methodologies of TSL00 established by Hara et al. (2019), including both the derivation of its formula and an efficient implementation. We will also briefly discuss HYDRA, a closely-related method proposed by Chen et al. (2021).

*Derivation of TSL00.* The high level idea of TSL00 is that, in order to obtain the LOO at the final model checkpoint, i.e., $\theta_{-j}^{[T]} - \theta^{[T]}$, it suffices to build a recursive relationship between the LOO at two consecutive iterations, $\theta_{-j}^{[t]} - \theta^{[t]}$ and $\theta_{-j}^{[t-1]} - \theta^{[t-1]}$.

We start with a simple case where the model is trained for only one epoch. Denote $\pi(j)$ as the timestamp where sample $z_j$ is used in training, i.e. $z_j \in S_{\pi(j)}$. For $t > \pi(j) + 1$, we have

$$ \theta_{-j}^{[t]} - \theta^{[t]} = (\theta_{-j}^{[t-1]} - \theta^{[t-1]}) - \frac{\eta_{t-1}}{|S_{t-1}|} \sum_{i \in S_{t-1}} (\nabla_{\theta} \ell(z_i; \theta_{-j}^{[t-1]}) - \nabla_{\theta} \ell(z_i; \theta^{[t-1]})). $$

Using first-order expansion, we have the following approximation:

$$ \frac{1}{|S_{t-1}|} \sum_{i \in S_{t-1}} (\nabla_{\theta} \ell(z_i; \theta_{-j}^{[t-1]}) - \nabla_{\theta} \ell(z_i; \theta^{[t-1]})) \approx H^{[t-1]}(\theta_{-j}^{[t-1]} - \theta^{[t-1]}). $$

Combining the above equations yields the following recursive relationship:

$$ \theta_{-j}^{[t]} - \theta^{[t]} \approx (I - \eta_{t-1} H^{[t-1]})(\theta_{-j}^{[t-1]} - \theta^{[t-1]}), \quad t > \pi(j) + 1. $$

Denote $Z_t = I - \eta_t H^{[t]}$, and note that $\theta_{-j}^{[\pi(j)+1]} - \theta^{[\pi(j)+1]} = \frac{\eta_{\pi(j)}}{|S_{\pi(j)}|} \nabla_{\theta} \ell(z_j; \theta^{[\pi(j)]})$, we have the following approximation of the LOO at the final model:

$$ \theta_{-j}^{[T]} - \theta^{[T]} \approx \frac{\eta_{\pi(j)}}{|S_{\pi(j)}|} Z_{T-1} Z_{T-2} \cdots Z_{\pi(j)+1} \nabla_{\theta} \ell(z_j; \theta^{[\pi(j)]}) =: \Delta \theta_{-j}. \quad (15) $$

<footnote n="13">¹³On the other hand, to date TSL00 supports only the SGD optimizer and does not accommodate alternatives like Adam (Kingma and Ba, 2014), which are essential for training language models.</footnote>
Intuitively, this formula can be interpreted as follows: starting from the first appearance of a sample, propagate its influence (scaled gradient) to the final iteration using second-order information. In data attribution literature, this approximation in Equation 15 first appears in Hara et al. (2019). Notably, similar analysis and expression appear frequently in related domains, including continual learning and deep learning theory (Zou et al., 2021; Evron et al., 2022; Wu et al., 2022a, 2024a; Ding et al., 2024).

For multi-epoch SGD, Hara et al. (2019) propose to sum up the influence of a sample at different epochs. Specifically, suppose SGD is applied for K epochs, and $\pi_k(j)$ is the timestamp of sample j's appearance at k-th epoch. We can similarly define $\Delta\theta_{-j}$ as

$$ \Delta\theta_{-j} = \sum_{k=1}^{K} \left( \prod_{s=1}^{T-\pi_k(j)-1} Z_{T-s} \right) \frac{\eta_{\pi_k(j)}}{|S_{\pi_k(j)}|} \nabla_{\theta} \ell(z_j; \theta^{[\pi_k(j)]}), $$

**Implementation of TSL00.** Hara et al. (2019) propose a backward computation algorithm that computes the TSL00 scores for all training samples in one pass. The key observation is that if the relevant information ($S_t, \eta_t, \theta^{[t]}$) is stored during training, it can later be used during inference to compute the (scaled) per-sample gradient and then its dot product with $u^{[t]}$, which can be sequentially updated from the initial query vector $u$. More specifically, starting from $t = T - 1$ and proceeding backward to $t = 1$, at each iteration the algorithm computes $\langle u, \frac{\eta_t}{|S_t|} \nabla_{\theta} \ell(z_j; \theta^{[t]}) \rangle$ for every $j \in S_t$ and adds it to the TSL00 score of sample j. Then, the query vector is updated as $u \leftarrow Z_t u$. Additionally, Hara et al. (2019) leverage off-the-shelf implementations of the Hessian-vector product to compute $Z_t u$, thereby bypassing the challenges associated with directly computing and storing the Hessian.

**Connection to HYDRA.** Chen et al. (2021) propose HYDRA, which is based on a similar principle to TSL00. Instead of computing the influence at the final model checkpoint directly via implicit differentiation, both methods unroll the differentiation throughout the optimization trajectory by establishing a recursive relationship between consecutive iterations. The key difference is that HYDRA considers the derivative of the upweighted objective as in influence functions, while TSL00 directly deals with LOO. In addition, HYDRA omits the Hessian terms in the recursive relationship for the purpose of computational efficiency, which results in a less precise approximation of the ground truth LOO.

**Scalable variations of TSL00.** While the vanilla TSL00 has been applied to GAN (Terashita et al., 2021), vanilla TSL00 methods face significant challenges when the number of model parameters scales up. Since vanilla TSL00 requires unrolling the training and leveraging information such as the Hessian, it posts two main challenges: high computation cost and large storage demand. In this section, we introduce two scalable variations (Wang et al., 2025b; Bae et al., 2024) of TSL00. We will mainly focus on their design and considerations to scale the method.

**DVEmb.** Wang et al. (2025b) points out that query vector $u$ is unknown at training time, while $\Delta\theta_{-j}$ only requires information along training. The authors then propose data value embedding, noted as DVEmb$^{(t_s)}(z_j)$, defined as

$$ \text{DVEmb}^{(t_s)}(z_j) = \eta_{t_s} \left[ \prod_{k=t_s+1}^{T-1} (I - \eta_k H^{[k]}) \right] \nabla_{\theta} \ell(z_j; \theta^{[t_s]}). $$

Notably, when $t_s = \pi(j)$, the definition of DVEmb$^{(t_s)}(z_j)$ is the same as Equation 15. Authors further leverage the empirical FIM to approximate the Hessian (see Section 2.2)

$$ H^{[k]} = \sum_{z \in S_k} \nabla_{\ell_{\theta}}(z, \theta^{[k]}) \nabla_{\ell_{\theta}}(z; \theta^{[k]})^{\top}, $$

which gives a recursive relationship for the proposed DVEmb$^{(t_s)}(z_j)$

$$ \text{DVEmb}^{(t_s)}(z^*) = \eta_{t_s} \nabla_{\ell_{\theta}}(z^*; \theta^{[t_s]}) - \eta_{t_s} \mathbf{M}^{(t_s)} \nabla_{\ell_{\theta}}(z^*; \theta^{[t_s]}), $$

where the matrix $\mathbf{M}^{(t_s)}$ is defined as:

$$ \mathbf{M}^{(t_s)} = \sum_{t=t_s+1}^{T-1} \left( \sum_{z \in S_t} \left( \text{DVEmb}^{(t)}(z) \nabla_{\ell_{\theta}}(z; \theta^{[t]})^{\top} \right) \right). $$

To compute it recursively:

$$ \mathbf{M}^{(t_s-1)} = \mathbf{M}^{(t_s)} + \sum_{z \in S_{t_s}} \text{DVEmb}^{(t_s)}(z) \nabla_{\ell_{\theta}}(z; \theta^{[t_s]})^{\top}, $$
and for $t = T - 1$, $M^{T-1} = 0$. The above framework allows computing $M$ and DVEmb starting from $t = T - 1$ to $t = 0$.

The method requires $O(TBp)$ storage during training, where $B$ is the batch size and $p$ is number of model parameters. This can be impractical especially for large models. To further improve the storage efficiency, the authors leverage gradient decomposition, storing only decomposed components instead of full gradient vectors, as in the ghost dot product (Equation 14), reducing storage requirements to $O(TB\sqrt{p})$ for non-sequential data. The author further includes dimensional reduction for foundation models so that the storage becomes $O(TB\tilde{p})$ where $\tilde{p}$ is significantly less than $p$.

The paper also includes some other considerations for computational efficiency. The authors propose to calculate data value embedding per layer so that the total FLOPs is reduced from $O(BT\tilde{p}^2)$ to $O(BT\tilde{p}^2/L)$ where $L$ is the number of layers. The authors also propose parallelized computation through checkpointing the models for fixed step intervals, which reduce the iteration numbers.

**SOURCE.** Bae et al. (2024) introduces an averaged notion of SGD-influence that averages across all possible training trajectories rather than focusing on a single specific trajectory. To efficiently approximate this averaged influence, they partition the training process into multiple segments (e.g., corresponding to different training phases) and approximate the Hessian and gradient computations within each segment.

The segment-based TSL00 can be presented as following:

$$ \mathbb{E} \left[ \frac{d\theta^{[T]}}{d\epsilon} \right] = -\mathbb{E} \left[ \sum_{\ell=1}^{L} \left( \prod_{\ell'=L}^{\ell+1} \mathbf{s}_{\ell'} \right) \underbrace{\left( \sum_{k=T_{\ell-1}}^{T_{\ell-1}} \frac{\eta_k}{B} \delta_k Z_{k+1:T_{\ell}} \mathbf{g}_k \right)}_{r_{\ell}} \right] \quad (16) $$

$$ \approx - \sum_{\ell=1}^{L} \prod_{\ell'=L}^{\ell+1} \mathbb{E}[\mathbf{S}_{\ell'}] \mathbb{E}[r_{\ell}], \quad (17) $$

where $Z_k = I - \eta_k H^{[k]}$ and $Z_{k:k'} = Z_{k'-1} \cdots Z_{k+1} Z_k$. $\mathbf{S}_{\ell'} = Z_{T_{\ell'-1}:T_{\ell'}}$ indicates the influence of segment $\ell'$, $r_{\ell}$ indicates the influence of the iterations in segment $\ell$, $\mathbf{g}_k = \nabla_{\theta} \ell(z; \theta^{[k]}).$

The approximation of stable Hessian and gradient within one segment $\ell$ greatly helps to further simplify $\mathbb{E}[\mathbf{S}_{\ell}]$ and $\mathbb{E}[r_{\ell}]$ to be

$$ \mathbb{E}[\mathbf{S}_{\ell}] = \mathbb{E}[Z_{T_{\ell-1}:T_{\ell}}] \approx (I - \bar{\eta}_{\ell} \bar{H}^{[\ell]})^{K_{\ell}} \approx \exp(-\bar{\eta}_{\ell} K_{\ell} \bar{H}^{[\ell]}) =: \bar{\mathbf{S}}_{\ell}, $$

$$ \mathbb{E}[\mathbf{r}_{\ell}] = \mathbb{E} \left[ \sum_{k=T_{\ell-1}}^{T_{\ell}-1} \frac{\eta_k}{B} \delta_k Z_{k+1:T_{\ell}} \mathbf{g}_k \right] \approx \frac{1}{N} (I - (I - \bar{\eta}_{\ell} \bar{H}^{[\ell]})^{K_{\ell}}) (\bar{H}^{[\ell]})^{-1} \bar{\mathbf{g}}_{\ell} \approx \frac{1}{N} \underbrace{(I - \exp(-\bar{\eta}_{\ell} K_{\ell} \bar{H}^{[\ell]})(\bar{H}^{[\ell]})^{-1} \bar{\mathbf{g}}_{\ell})}_{=: \bar{\mathbf{r}}_{\ell}}. $$

Then the whole procedure has the following close form

$$ \mathbb{E} \left[ \frac{d\theta_T}{d\epsilon} \right] \approx -\frac{1}{N} \sum_{\ell=1}^{L} \left( \prod_{\ell'=L}^{\ell+1} \bar{\mathbf{s}}_{\ell'} \right) \bar{r}_{\ell}, $$

Since the calculation of earlier segments requires information from later segments, SOURCE features backward calculation: for later segments, the approximated Hessian product will be passed to earlier segments to assist calculation. Furthermore, SOURCE also proposes several techniques to scale the method. Firstly, by checkpointing, we require fewer calculations of Hessian matrices. Secondly, in order to efficiently calculate the Hessian and its product, SOURCE leverages Eigenvalue-corrected Kronecker-Factored Approximate Curvature (EK-FAC) (Grosse et al., 2023) to assume the Hessian. The authors also present a variant of SOURCE, called Fast-Source, in which the parameters are averages in a segment rather than Hessian and gradient.
### 2.4.3 Key Takeaways and Future Directions

Training dynamic methods are grounded in two seminal studies: TracIN (Pruthi et al., 2020) and SGD-Influence (Hara et al., 2019). Tracing gradient descent methods are popular due to their ease of implementation and relatively low computational cost, making them an active area of research for adaptation to various generative AI models. However, their performance can be compromised by spurious factors and often lags behind other methods in quantitative evaluations. Trajectory-specific leave-one-out (TSL00) methods avoid some of the restrictive assumptions inherent in influence functions, such as convex loss functions and convergence during training, making them better suited for generative AI settings. Nevertheless, their scalability to large models remains a significant computational challenge, limiting their practical applicability. A common future direction for training dynamic methods is enhancing compatibility with modern optimizers. Both tracing gradient descent and Trajectory-specific LOO methods are derived under the assumption that stochastic gradient descent (SGD) is used, whereas SGD is rarely the optimizer of choice for training generative AI models. Xia et al. (2024) have taken initial steps toward addressing this mismatch, but many open questions remain.

## 2.5 Category: Simulator

This category of methods represents a critical line of research that uses surrogate models (i.e., simulators) to estimate the counterfactual effect of model training. Datamodels (Ilyas et al., 2022) is a landmark work in this category, introducing a conceptual framework for analyzing model behavior using a surrogate model and proposing linear datamodel as a practical surrogate function. Other methods in this category have proposed alternative simulators with varying input-output structures (Guu et al., 2023; Chai et al., 2024) and simulator architectures (Yoon et al., 2020; Liu et al., 2025). This section will also highlight the close connection between Linear Datamodel and semivalues (e.g., Data Banzhaf (Wang and Jia, 2023a)) in weighted marginal contribution category as well as how different simulation models extend the conceptual framework of Datamodels.

### 2.5.1 Datamodels

**Notation.** Suppose we have training dataset $D = \{z_i\}_{i=1}^n$. We extend the notation of target function as $f(z_{test}; \hat{\theta}_S)$ indicating the model output $f$ on test data sample $z_{test}$ and the model $\hat{\theta}_S$ trained on a subset of training dataset $S \subseteq D$.

**Datamodeling.** Ilyas et al. (2022) propose a data attribution framework, datamodeling, to quantify the influence of training data on model outcomes through surrogate functions. Specifically, Datamodels for a test sample $z_{test}$ is a parametric function $g_\beta: 2^D \to \mathbb{R}$ optimized to predict $f(z_{test}; \hat{\theta}_S)$ given a training subset $S \subseteq D$:

$$g_{\beta}: \{0, 1\}^n \to \mathbb{R}, \quad \text{where} \quad \beta = \arg \min_{\beta} \mathbb{E}_{S \sim D}^{(M)} \left[ \mathcal{L}_{\text{sim}}(g_{\beta}(\mathbf{1}_S); f(z_{\text{test}}; \hat{\theta}_S)) \right].$$

Here, $D$ defines a sampling distribution of the subset $S$, a fixed training set. $\mathbf{1}_S \in \{0, 1\}^n$ is the characteristic vector of $S$, $\mathcal{L}_{\text{sim}}(\cdot; \cdot)$ is a loss function for simulator training, and $\mathbb{E}^{(M)}$ is an $M$-sample empirical estimate of the expectation.

**Linear Datamodel.** In practice, Ilyas et al. (2022) propose to use a linear function as the surrogate function $g_\beta$, i.e., $g_\beta(S) = \beta_{[1:]}^\top \mathbf{1}_S + \beta_0$ for some $\beta \in \mathbb{R}^{n+1}$. This is referred to as the Linear Datamodel. The additive assumption simplifies the analysis by assuming that the combined effect of data can be presented as the simple summation of independent contributions from each point. Once this assumption is in place, a linear datamodel can be obtained via solving the following LASSO optimization problem (Tibshirani, 1996):

$$\beta = \arg \min_{\beta \in \mathbb{R}^{n+1}} \frac{1}{M} \sum_{k=1}^{M} \left( g_{\beta}(\mathbf{1}_{S_k}) - f(z_{\text{test}}; \hat{\theta}_{S_k}) \right)^2 + \lambda \|\beta\|_1. \quad (18)$$

The subsets $S_1, \dots, S_M$ are sampled from a uniform distribution over $\alpha$-fraction subsets of the training set, i.e., $D = \mathcal{U}(S \subseteq D : |S| = \alpha n)$, which is a simple design choice made by Ilyas et al. (2022). $\lambda \ge 0$ is the regularization coefficient. The sparse prior regularization helps prevent overfitting and reduces the number of training needed. Finally, $\beta_i$ will be the attribution score for the $i$-th training sample in the Linear Datamodel.

**Connection with weighted marginal contribution methods.** It is worth noting that the formulation in equation 18 is similar to the sparse regression estimator proposed in AME (Lin et al., 2022). We will further discuss the connection
between Linear Datamodel and weighted marginal contribution category. In general, Linear Datamodel can be written as linear combinations of utility values, or target function values, with proper designed subset sampling distribution $\mathcal{D}$.

When the $l_1$ regularization term is removed, we can reformulate Equation 18 as a weighted least-squares problem with a subset sampling distribution $\mathcal{D}$:

$$ \arg\min_{\beta} \sum_{s \in [k]} w_s \left[ (\beta^{\top} \mathbf{1}_{S_s} - f(z_{\text{test}}; \theta_{S_s}))^2 \right]. \quad (19) $$

where $k$ subset are sampled from distribution $\mathcal{D}$, $w_s \ge 0$ and $\sum_{s=1}^{k+1} w_s = 1$. For Linear Datamodel with uniform sampling distribution, the weights are defined as

$$ w_s = \begin{cases} \frac{1}{(\alpha n)}, & \text{if } |S| = \alpha n, \\ 0, & \text{otherwise.} \end{cases} $$

Considering the stacked vector $\mathbf{1}_S$ in $X, W$ as the diagonal matrix with each element $w_s$, $f(z_{\text{test}}; \theta) \in \mathbb{R}^n$ as an output vector for all $S$, the optimal solution for this weighted least-squares problem is $\beta^* = (X^{\top}W X)^{-1}X^{\top}W f(z_{\text{test}}; \theta)$ provided $X^{\top}W X$ is invertible. The optimal solution $\beta^*$ of these linear systems satisfies several key properties:

1. Under certain conditions on $w_s$, $\beta^*$ preserves the pairwise differences of certain probabilistic values, that is, it retains the relative order of its components (Li and Yu, 2024b).

2. By placing a suitable distribution over the subsets (i.e., specific choice of $w_s$), $\beta^*$ can recover known probabilistic values. In particular, setting $w_S = \frac{1}{2^n}$ and transforming $X$ into a matrix with elements $\pm 1$, we can recover Data Banzhaf Wang and Jia (2023a) (detailed in Appendix B). More related discussion can be found in (Li and Yu, 2024b).

3. For Linear Datamodel, the derived data values satisfies the linearity and symmetry axioms. Since $\beta^*$ is derived via linear transformations of $f(z_{\text{test}}; \theta)$, linearity holds naturally. Symmetry follows because the weights $w_s$ depend only on the subset size. If two data points contribute identically to $f(z_{\text{test}}; \theta)$ across all subsets, $X^{\top}W X$ and $X^{\top}W f(z_{\text{test}}; \theta)$ will remain symmetric w.r.t $i$ and $j$, leading to $\beta_i^* = \beta_j^*$.

Ilyas et al. (2022) also discusses the connection of Linear Datamodel and empirical influence estimation proposed by Feldman and Zhang (2020). Specifically, taking $(f(z_{\text{test}}; \theta_S)$ as the correctness of predicting test sample $z_{\text{test}}$, the empirical influence estimation of $i$-th training data point on the test sample $z$ can be approximated by:

$$ \mathbb{E}_{S \sim \mathcal{D}}(f(z; \theta_{S \cup \{z_i\}})) - \mathbb{E}_{S \sim \mathcal{D}}(f(z; \theta_S)). $$

Ilyas et al. (2022) shows that with $\alpha = 0.5$ and transforming $X$ into a matrix with elements $\pm 1$, the minimizer of original linear regression in equation 19 can cast the estimated empirical influence in the infinite-sample case.

**Challenges.** There are two notable limitations to this approach of solving linear datamodel. First, the computational cost remains significant since, in equation 18, obtaining the ground truth labels requires training models on $M$ subsets, and the dimension of $\theta$ (i.e., the number of samples) is large. To resolve this problem, some gradient based works (Park et al., 2023) are proposed with promising results on Linear Datamodeling Scores (LDS), an evaluation metric proposed under the concept of Linear Datamodel. The method is introduced in Section 2.2.4 as its close relationship with influence function.

Second, the linear approximation of $g_\beta$ may introduce misspecification errors, as it oversimplifies the potentially complex relationships in the data. These challenges highlight the need for further advancements to improve computational efficiency and reduce approximation errors.

### 2.5.2 Other Simulators

Other than Ilyas et al. (2022), different simulator-based data attribution methods are proposed to have different inputs, outputs and architectures in different studies. Two popular paradigms include 1) predicting the counterfactual target function (Yoon et al., 2020; Liu et al., 2025); 2) predicting the target function trajectory (Guu et al., 2023; Chai et al., 2024) are designed to resolve shortcomings in gradient-based data attribution methods.

**Predict the counterfactual target function output.** A group of methods that directly simulate the counterfactual target function output follow a simple idea that a map between training subset and target function output can be learned
by a neural network. The group of methods is closely related to Linear Datamodel (Ilyas et al., 2022), which can be seen as taken a linear model to predict the counterfactual target function.

DVRL (Data Valuation using Reinforcement Learning) (Yoon et al., 2020) leverages reinforcement learning to train a data value estimator (DVE). The method consists of two components, the target task predictor model and DVE that evaluates each data example. On each training step, a target task predictor model is trained by minimizing the loss on the training set, weighted by the data valuation learned by DVE model. After that, the DVE model is optimized on a small validation set to minimize the loss through adjusting the data weight. The DVE model is trained with RL method policy gradient, receiving the reward signal from the validation loss and its moving average.

RegMix (Liu et al., 2025) is another work that leverages multiple retraining on subset to collect training data for the simulator to directly predict counterfactual model performance. The major difference between Liu et al. (2025) and Ilyas et al. (2022) is that RegMix take the data mixture percentage as simulator input rather than the sample-level information.

CLIMB (Diao et al., 2025) automatically clusters large datasets in semantic space and uses a small proxy model and a predictor to iteratively optimize the data mixture in the pre-training setting of LLMs.

**Predict the training trajectory.** Another group of method predict the training trajectory in each step, such as simfluence (Guu et al., 2023) and GPTfluence (Chai et al., 2024).

Simfluence (Guu et al., 2023) adopt a linear first-order Markov process to model the test loss throughout the training trajectory. The model takes in the curriculum, $C_t$, which is the order of the training batches on time step $t$ fed in, and the target function output before training step $t$ begins and output predicted test loss after this training step. Here we simplify the target function to be the test loss on a test data point $\ell(z_{test}; \theta)$.

Specifically, the simulator is modeled as:

$$ \ell(z_{test}; \theta_t) = \alpha(C_t)\ell(z_{test}; \theta_{t-1}) + \beta(C_t), $$

where $\alpha(C_t) = \sum_{i \in C_t} A_i$, $\beta(C_t) = \sum_{i \in C_t} B_i$, and $A_i, B_i \in \mathbb{R}^n$. The reduction in loss is composed by a relative reduction $\alpha(C_t)$ and an absolute reduction $\beta(C_t)$.

Similar to Simfluence (Guu et al., 2023), GPTfluence (Chai et al., 2024) propose another simulator model to predict the training trajectory. There are two major differences between Simfluence and GPTfluence. First, instead of a first order Markov process, the authors use an $n$-th order Markov process to model the training trajectory, which introduces more parameters but allows for modeling longer dependencies. Specifically, the simulator is modeled as:

$$ \ell(z_{test}; \theta_t) = \sum_{j=1}^{n} \alpha_j(C_t) \ell(z_{test}; \theta_{t-j}) + \beta(C_t). $$

Note that the process is of $n$-order, meaning that the latest target function output is related to the previous $n$ steps' target function output. Moreover, GPTfluence also proposed to a pre-trained encoder such as BERT and GPT to preserve the information within the data examples as well as support unseen training data in $C_t$. The coefficients of the Markov process $\alpha(c_t)$ and $\beta(c_t)$ is then obtained by taking the inner product of these embeddings after some linear layers.

### 2.5.3 Key Takeaways and Future Directions

Training a simulator to directly predict attribution results represents a conceptual framework aimed at learning the relationship between training examples and model behavior. Among studies in this direction, Datamodels stands out as the most popular one. The primary challenge for this category of methods lies in the need to generate training data for the simulator—typically through model retraining. This requirement severely limits the scalability of simulator-based approaches to generative models. Although predicting the training trajectory has been proposed as a way to reduce retraining costs, it introduces another practical concern: error accumulation during trajectory prediction, which can significantly impact the reliability of the results.

## 2.6 Category: Other Approaches

There are also some methods that could not be grouped to previous categories. In this section, we collect the untouched methods and identify two rough groups: unlearning-based methods and similarity-based methods.
## 2.6.1 Unlearning-based Methods

A group of methods follow a general procedure that unlearn certain data points from the fully trained model and directly examine the counterfactual prediction through the unlearned model. Different methods are proposed to unlearn one or multiple data points, which can be broadly categorized into “unlearning training samples” and “unlearning test samples”. Notably, most studies in this category focus on data attribution for diffusion models, likely due to the ambiguous definition and intriguing discover of the target function of diffusion model (Zheng et al., 2023; Lin et al., 2025) in gradient-based data attribution methods.

**Unlearning training samples.** A typical path is to retrain or fine-tuning on different subsets of training data (Lu et al., 2025; Dai and Gifford, 2023) as unlearning. The resulting unlearned models can then be used to estimate the marginal impact of individual training entities.

Dai and Gifford (2023) propose dividing the entire training set into $n$ overlapping subsets, where $n$ diffusion models are individually trained on these subsets. During the generation process, images are generated using an ensemble of these diffusion models. To unlearn a specific training example, all models that have seen the example are removed from the generation process. Similarly, Lu et al. (2025) fine-tune a fully trained model on different contributor-specific data subsets and compute the Shapley value for each contributor.

The path is closely connected to weighted marginal contribution and share some common challenges, where high computational cost of retraining and fine-tuning is the most critical one. To mitigate fine-tuning costs, Lu et al. (2025) first prune the model to approximate the fully trained version before fine-tuning.

**Unlearning test samples.** Another way is leveraging the mirrored influence hypothesis (Ko et al., 2024), which suggests a strong correlation between train-to-test and test-to-train influence. This insight enables the estimation of a training sample's marginal impact by unlearning the corresponding test sample (Wang et al., 2024h).

Similar to influence functions, Wang et al. (2024h) propose an attribution method for diffusion models. The key quantity of interest is the change in the loss of a training sample before and after unlearning a specific generated sample. Since retraining is computationally prohibitive, Wang et al. (2024h) employ unlearning techniques to approximate the removal effect. Specifically, they minimize the elastic weight consolidation (EWC) loss using multiple iterations of Newton update to mitigate the influence of the generated example.

$$ \theta \leftarrow \theta + \frac{\alpha}{N} F^{-1} \nabla \ell(\hat{z}; \theta), \\ \tau(\hat{z}, z) = \ell(z; \theta_{-\hat{z}}) - \ell(z; \theta_0), $$

where $F$ is the Fisher Information Matrix, $\ell$ is the training loss, $\hat{z}$ is the generated sample and $z$ is a training sample.

## 2.6.2 Similarity-based Methods

Similarity-based data attribution methods aim to estimate the influence of data through a defined notion of similarity, following the intuition that influential training data should resemble test data. Based on this principle, several studies have proposed model-agnostic approaches for data attribution. However, these model-agnostic methods do not strictly fall within the scope of this survey, as our definition of data attribution in Section 2 focuses on model behavior. Therefore, we include only a few representative studies to illustrate research efforts in this direction.

Just et al. (2023) propose using multi-class Wasserstein distance to quantify the similarity between two datasets. They frame the data valuation problem as assessing how changes in the distribution of training data affect the similarity between the training and validation sets. Each training sample's contribution is then determined by computing the gradient of this distance with respect to perturbations in probability mass. Compared to LAVA, Yang et al. (2025) adopt a discrete similarity-based approach. Intuitively, they compute the value of a training sample as its calibrated distance to the entire validation set. Representer Point Selection (RPS) (Yeh et al., 2018) also employs weighted feature similarity to represent the importance of the training samples. Specifically, leveraging the representer theorem, it uses pre-activate input of model's last linear classification layer as feature and the resistance for training example feature towards minimizing the norm of the last layer's parameters as weight.
### 2.6.3 Key Takeaways and Future Directions

Machine unlearning is closely related to data attribution in terms of methodology, and thus faces similar challenges, such as the high computational cost of retraining and the non-convexity issues in gradient-based unlearning. The *mirrored influence hypothesis* offers a promising direction for reducing computational cost, particularly in scenarios with a small number of test samples and an extremely large training dataset. Similarity-based methods are generally model-agnostic; although they do not strictly align with the definition of data attribution, they are relatively easy to implement and can perform well in certain downstream tasks, such as mislabeled data detection.

## 3 Applications

In this section, we explore the applications of data attribution, focusing on generative AI settings. We categorize the applications into four primary areas: *data selection, fact tracing and concept probing, adversarial attack and defense*, and finally *data economy*. In addition, we discuss a few emerging but less common applications in Section 3.5.

For each application, we begin by defining the problem setup and objectives, emphasizing how data attribution principles can help address its specific challenges. We then review and organize relevant studies that employ attribution-based methods to tackle these challenges. In applications where literature in a generative AI context is still limited—given the recent emergence of attribution methods for generative models—we focus on conventional machine learning approaches instead, including them with the expectation that their insights can be transferred to generative AI settings. Finally, we conclude with a comprehensive synthesis of the current landscape and future directions. This closing section summarizes existing literature, reflects on unique difficulties and open challenges, considers practical aspects of applying attribution methods, and offers recommendations for future research and development.

### 3.1 Data Selection

Modern generative AI models are typically trained on extensive datasets to improve performance. However, not all data in these large corpora contribute positively to training outcomes: some examples may be irrelevant, mislabeled (Brodley and Friedl, 1999), redundant (Li et al., 2023), or even harmful (Gong et al., 2023). *Data selection* focuses on identifying and retaining the most useful data while filtering out less relevant portions (often referred to as *data cleaning*), with the aim of improving model accuracy and increasing training efficiency by reducing dataset size.

#### 3.1.1 Problem Settings

Data selection aims to construct a curated training dataset $D$ from a larger raw data pool $D_{\text{raw}}$, such that the model $M_D = A(D)$ trained by algorithm $A$ on $D$ achieves desirable properties, e.g., low evaluation loss, better generalization, or safety guarantees. Formally, this corresponds to minimizing an objective function $f_{\text{obj}}$ over possible datasets: $D^* = \arg \min_D f_{\text{obj}}(D)$, where $D$ is drawn from a feasible subset of $D_{\text{raw}}$ (Albalak et al., 2024). Notably, the selected dataset $D$ need not be a strict subset; some approaches allow for weighted sampling or duplication, meaning $D$ may contain repetitions or weightings of examples from $D_{\text{raw}}$. In general, the combinatorial nature of the search space renders exact optimization intractable.

To address this, a variety of approximation methods have been proposed. Recent work, in particular, has leveraged data attribution techniques, which estimate the *counterfactual* influence of individual training examples (or groups thereof) on a specified model behavior. These influence scores serve as tractable surrogates for the objective function and can guide data selection through computationally efficient heuristics. Thus, data attribution offers a scalable and principled approach to approximating the intractable selection problem.

Such approaches have been explored across a wide range of generative AI applications, varying in models and their training setup (e.g., LLM pre-training/fine-tuning, diffusion models), selection granularity (e.g., per-example, per-domain), and target objectives (e.g., model performance, safety). Before discussing the instantiations across these different settings, we first present a general algorithmic abstraction of data-attribution-based selection below.

**General data-attribution-based data selection pipeline.** To understand the core design principles of data-attribution-based selection, we begin by identifying the main computational bottlenecks in directly solving $\arg \min_D f_{\text{obj}}(D)$. First, the combinatorial explosion of candidate datasets $D$ renders exhaustive search intractable. Second, evaluating $f_{\text{obj}}(D)$ even once typically requires full training and evaluation, which is prohibitively expensive in large-scale generative models. In practice, only a limited number of such “oracle” evaluations of $f_{\text{obj}}$ are computationally feasible.
Data attribution techniques address both challenges by replacing exhaustive evaluation with efficient *counterfactual prediction*, followed by a computationally tractable selection strategy. Conceptually, counterfactual prediction provides a localized estimate of each training unit’s influence, and the selection strategy operates on these scores to curate the final dataset *D*. We note that this general framework has been widely explored in recent data selection work, particularly in the context of large language models (Albalak et al., 2024).

In what follows, we focus on the distinctive features of data-attribution-based methods for both counterfactual prediction and the selection strategy. The former offers a more principled and interpretable approach to data selection compared to heuristic or proxy-based alternatives. The latter, while more general and not unique to attribution-based pipelines, remains integral to the overall effectiveness and trade-offs of the framework.

*Counterfactual prediction.* Instead of directly evaluating the objective for every possible dataset *D*, counterfactual prediction considers a localized, fine-grained target function $f(S)$, commonly defined as

$$f(S) := f_{\text{obj}}(D_{\text{raw}} \setminus S),$$

where *S* is a subset of the training data, such as a single point, a group of examples, or an entire domain. This function captures the sensitivity of $f_{\text{obj}}$ to small perturbations of the dataset and reflects the counterfactual change in the objective when *S* is removed. Specifically, the corresponding counterfactual change is referred to as the *influence* of *S*, defined as

$$\tau(S) := \Delta f(S) = f(S) - f(\emptyset) = f_{\text{obj}}(D_{\text{raw}} \setminus S) - f_{\text{obj}}(D_{\text{raw}}),$$

which quantifies how much the objective would degrade (or improve) if *S* were removed. This measures the contribution of training units to the property encoded in the objective. For example, when the goal is to improve model performance, one may define the objective as $f_{\text{obj}}(D) = \ell(D_{\text{eval}}; M_D)$, where $\ell$ is a loss function evaluated on a fixed evaluation set $D_{\text{eval}}$. A natural choice of target function is then $f(S; D_{\text{eval}}) := \ell(D_{\text{eval}}; M_{D_{\text{raw}}\setminus S})$, and the corresponding influence becomes $\tau(S) = \ell(D_{\text{eval}}; M_{D_{\text{raw}}\setminus S}) - \ell(D_{\text{eval}}; M_{D_{\text{raw}}})$, capturing the marginal impact of *S* on evaluation performance.

Crucially, various *data attribution* algorithms exist to approximate $\tau(S)$ without explicitly retraining models for each possible subset *S*, which enables scalable, localized assessment of individual data units’ impact on the objective. We note that the specific choice of data attribution algorithm depends on several factors, including the model training setting (e.g., pre-training, fine-tuning), the form of the target function (e.g., evaluation loss, likelihood), and the granularity of attribution (e.g., individual points, groups, or domains). With a choice of data attribution method, the resulting estimated influence scores $\tau(S)$ then serve as surrogates to approximate the original dataset-level optimization objective in a computationally tractable manner.

*Selection strategy.* The final ingredient in the data-attribution-based data selection pipeline is to act on the estimated influence scores by deciding which data to retain, discard, or re-weight. This is governed by a *selection strategy*, which translates local influence estimates into global modifications of the dataset. Although not unique to attribution-based methods, the selection strategy plays a crucial role in determining how effectively we can utilize the localized estimates $\tau(S)$ to approximately maximize the original dataset-level objective $f_{\text{obj}}(D)$. As the influence scores are defined for every small training unit *S* (e.g., individual examples or groups), several widely used strategies arise naturally by making decisions based on these individual scores:

* Threshold selects data units whose influence exceeds (or falls below) a fixed or adaptive threshold. This avoids combinatorial search by assuming that influence scores provide a meaningful separation between helpful and harmful data.

* Top-k selection fixes the number of retained examples, selecting the k highest-ranked units by influence. This strategy assumes that influence magnitudes provide a globally consistent ordering of utility.

* Importance weighting retains all data but adjusts each unit’s contribution according to its influence. This relaxes the binary nature of selection and assumes that marginal improvements can be linearly composed.

These strategies each reflect specific structural or heuristic approximations to the original dataset-level optimization, and entail trade-offs: threshold is simple but may be brittle under noisy influence estimates; top-k offers explicit control over dataset size but may discard complementary examples; importance weighting preserves more information but depends critically on the fidelity of influence estimation. Crucially, all of these approaches rely on a common modeling assumption: that each training unit’s influence on the objective can be assessed independently of others. That is, we assume $\tau(S \cup S') = \tau(S) + \tau(S')$, or at least that $\tau$ is highly additive or positively correlated across units, so that making local decisions independently (i.e., selecting *S* into *D* or not) yields a near-optimal global outcome $f_{\text{obj}}(D)$. This additive or decomposable view of influence enables tractable approximation of the otherwise combinatorial selection problem. However, this assumption may not always hold: for instance, when the utility of one point depends on the presence of others (Hu et al., 2024a).
Lastly, we also highlight that selection strategies can be deployed in either offline or online settings. In the offline setting, data is selected prior to training based on influence scores computed using a fixed model, typically trained on a subset of the full dataset as a proxy. This setting offers computational efficiency and is well-suited to static selection. In contrast, online selection interleaves influence estimation and model training, often by periodically re-evaluating scores as the model evolves. This enables more adaptive selection that can account for shifting data utility during training, and also helps mitigate the additive limitation discussed above by dynamically updating scores to partially capture higher-order interactions. The two approaches can also be integrated: for example, by initializing with an offline selection and refining it online.

### 3.1.2 Landscape of Data-Attribution-Based Data Selection in Generative AI

In this section, we examine how concrete data selection algorithms instantiate key aspects of data attribution methods— counterfactual prediction and selection strategy—across diverse generative AI training settings, to achieve the specified objective. Specifically, in each setting, we structure the discussion as follows: we first identify the general objectives and key obstacles shared in this setting, and discuss how different data attribution algorithms and selection strategies are utilized to overcome them together with specific considerations that are unique to different works.

**Data selection for LLM pre-training.** In LLM pre-training, the goal is to acquire LLMs with “general knowledge” from the massive text corpora collected from various sources and domains. Due to the massive raw training dataset $D_{\text{raw}}$, the main bottleneck at this stage is the intense computational resources required for a pre-training run. Therefore, data selection becomes a natural technique to reduce the training cost by reducing the dataset size.

In the case of data-attribution-based selection methods, as pre-training focuses on acquiring “general knowledge,” the objective $f_{\text{obj}}(D)$ is often the validation loss $\ell(D_{\text{eval}}; M_D)$, with the natural target function $f(S) = f_{\text{obj}}(D_{\text{raw}} \setminus S) = \ell(D_{\text{eval}}; M_{D_{\text{raw}} \setminus S})$ as discussed above (Zhang et al., 2025a; Wang et al., 2024e; Liu et al., 2025). In some special cases when the downstream tasks are known in advance, the evaluation set $D_{\text{eval}}$ can be drawn specifically from the task domains (Wang et al., 2023c; Yu et al., 2024). On the other hand, Thakkar et al. (2023) approaches the problem with a complementary viewpoint by focusing on identifying noisy and outlier data. In such a scenario, self-loss $f(S) = \ell(S, M_{D_{\text{raw}} \setminus S})$ becomes a natural target function with its induced influence score being the so-called self-influence since it has been observed that self-influence is a strong indicator of outlier and noisy data (Bejan et al., 2023).<fnref n="14" />

Given a target function, various attribution methods can be applied to estimate counterfactual effects, which are then used in conjunction with a selection strategy to curate the final dataset. However, due to the high computational cost of estimating counterfactual influence scores—especially in large-scale generative models with massive raw training dataset—the main consideration in choosing a specific attribution method and selection strategy is minimizing additional overhead. Below, we summarize three common optimization strategies for achieving this goal:

*   **Efficient influence approximation:** The most direct approach is to reduce the computational burden by accelerating the attribution step itself. A wide range of methods have been proposed to scale up standard attribution algorithms (Section 2), including utilizing proxy models, exploiting model architectures, or more aggressive approximations that trade off prediction fidelity for speed. We will expand on these techniques in later paragraphs.

*   **Coarser granularity:** Most prior work adopts per-example attribution (i.e., $|S| = 1$), which requires estimating $|D_{\text{raw}}|$ separate influence scores. A natural way to reduce this cost is to increase the attribution granularity, e.g., performing attribution over groups or domains (i.e., $|S| > 1$), which reduces the total number of training units that must be evaluated and selected.

*   **Reduced frequency of selection:** Another orthogonal strategy is to minimize the number of times attribution and selection are performed. Offline selection methods typically require only a single attribution pass before pre-training, whereas online methods continuously re-estimate influence during training. While online strategies may improve adaptivity, they incur substantially higher overhead due to repeated evaluation.

Recent works explore various combinations of these optimization strategies to balance computational efficiency with attribution fidelity. Table 8 summarizes representative methods applied to LLM pre-training, where bolded entries indicate those that explicitly incorporate one of the above strategies to reduce computational cost. Note that this is omitted in the **Attribution Method** column, as almost all works have a certain level of optimization to one of the standard attribution algorithms. In particular, we highlight several representative works that exemplify distinctive pipelines or introduce algorithmic innovations worth discussing from various angles.

<footnote n="14">¹⁴As an intuitive explanation is that, if removing a data point deteriorates the loss value on itself, then the data should be different enough so that the prediction on which could not be learned from the rest of the data, making it an outlier in most cases.</footnote>
<table><caption>Table 8: Comparison of data-attribution-based methods for LLM pre-training.</caption><thead><tr><th>References</th><th>Attribution Method</th><th>Target Function</th><th>Granularity</th><th>Selection Strategy</th></tr></thead><tbody><tr><td>Thakkar et al. (2023)</td><td>TracIn</td><td>Self-Loss</td><td>Per-example</td><td>Re-weighting &amp; Online</td></tr><tr><td>Wang et al. (2023c)</td><td>Grad-Dot</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Wang et al. (2024e)</td><td>(Higher-order) Grad-Dot</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Online</td></tr><tr><td>Yu et al. (2024)</td><td>BERT-based Datamodel</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Online</td></tr><tr><td>Engstrom et al. (2024)</td><td>TRAK on Proxy Model</td><td>Validation Loss</td><td>Subset</td><td>Top-k &amp; Offline</td></tr><tr><td>Zhang et al. (2025a)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Liu et al. (2025)</td><td>Regression Datamodel</td><td>Validation Loss</td><td>Domain</td><td>Re-weighting &amp; Offline</td></tr></tbody></table>

Firstly, a popular research direction has focused on improving the efficiency of gradient-based attribution methods by exploiting structural properties of gradients, often in the form of *Kronecker factorization*, to speed up influence computation. For example, Wang et al. (2024e) proposes the *ghost dot-product* for linear layers, which enables efficient influence estimation (Equation (14)). Similarly, Zhang et al. (2025a) leverages a Kronecker factorization structure tailored to attention layers, achieving comparable acceleration in transformer-based models. These algorithmic refinements lead to notable speedups compared to their original attribution counterparts.

Next, we revisit the distinction between online and offline selection strategies and how attribution methods adapt in each case. As discussed earlier, online selection has the conceptual advantage of using the current model checkpoint as the reference model for attribution. This makes it natural to apply *dynamic* attribution methods, which are designed to trace the behavior of a specific training run. For example, Thakkar et al. (2023) uses TracIn to compute self-influence scores for individual samples throughout training, enabling an adaptive re-weighting scheme that evolves over time. Along similar lines, Yu et al. (2024) trains a separate influence prediction model to approximate oracle scores in real time. This model is combined with a Gumbel-top-k sampling mechanism to introduce stochasticity into data selection, thereby balancing exploration and exploitation during pre-training.

In contrast, offline strategies face a conceptual challenge: since training has not yet started, no reference model $M_{D_{raw}}$ exists to support counterfactual prediction. A common workaround is to train small proxy models at low cost and use them for attribution. For instance, Liu et al. (2025) frames domain-level selection as a regression task, using small models trained on various data mixtures to estimate domain influence. Similarly, Engstrom et al. (2024) trains multiple GPT2-small variants on diverse subsets and applies TRAK to identify influential subsets for pre-training. An alternative approach is to assume access to a partially trained model: both Zhang et al. (2025a) and Wang et al. (2023c) consider a setting in which a pre-trained model already exists, and the goal is to refine it further. They hypothesize that influence estimates computed with respect to this intermediate model are sufficiently accurate to guide effective data selection.

**Data selection for LLM fine-tuning.** LLM fine-tuning focuses on improving generalization and robustness beyond a pre-trained foundational model by training on task-specific examples. Compared to pre-training, this setting typically involves domain-specific, high-quality *annotated* data and operates at a smaller scale, making it more controlled and selective. As a result, more computational resources can be allocated to the data selection process, enabling more precise and confident decisions to maximize data efficiency.

In the context of data-attribution-based selection for LLM fine-tuning, the choice of target functions is more diverse than in pre-training, as overall task performance is not always the sole objective. While validation loss on a general held-out evaluation set $D_{eval}$ remains a common choice (Jiao et al., 2025; Chhabra et al., 2024; San Joaquin et al., 2024; Xia et al., 2024; Wang et al., 2024g; Cao et al., 2024b; Zhou et al., 2024), other target functions have been explored to reflect different priorities. For example, Zhao et al. (2024) considers federated fine-tuning, where private client data cannot be shared; accordingly, the target function is defined using a *public* validation set. He et al. (2024) focuses on safety-critical applications, where the validation set is curated to emphasize safety-sensitive examples. Han and Tsvetkov (2021), in contrast, targets robustness to spurious correlations: their objective is to minimize loss on a held-out set that reflects only true task-relevant features, thereby penalizing training examples that strongly influence spurious signals. To summarize, although validation loss remains the predominant choice of target function in LLM fine-tuning, the design and curation of the evaluation set can vary significantly depending on task-specific goals, such as safety, privacy, or robustness, leading to different downstream data selection behaviors.

With the target function defined, various attribution methods are applied to conduct the corresponding counterfactual predictions, followed by a selection strategy to curate the final dataset. Interestingly, in contrast to the optimization strategies commonly explored for LLM pre-training—such as *efficient influence approximation* and *coarser granularity*—
data selection for fine-tuning often takes a different direction, emphasizing *precision* over efficiency due to the smaller and more curated nature of $D_{raw}$. Two common design choices reflect this shift:

*   *Precise influence estimation:* Unlike in pre-training, where attribution efficiency is a primary concern, fine-tuning literature places less emphasis on speed and instead adopts computationally expensive methods that align closely with task objectives. Classical influence functions, Shapley-value-based methods, and dynamic approaches are all explored in this setting.

*   *Finest granularity:* Because fine-tuning datasets are typically small and costly to collect, *per-sample* attribution becomes a natural choice, as this granularity enables highly targeted and precise selection decisions to achieve data efficiency.

We remark that offline methods remain widely used: unlike in pre-training, the absence of a reference model is not a bottleneck in fine-tuning since the pre-trained model itself provides a strong initialization. Since the fine-tuned model is often assumed to remain close to its pre-trained counterpart, attribution scores computed *a priori* with respect to the pre-trained model tend to remain valid throughout the fine-tuning process. This makes offline, one-shot estimation of influence both practical and reliable in this regime. In summary, recent work often combines both of the above strategies to enable accurate data selection in data- and resource-constrained environments.

<table><caption>Table 9: Comparison of data-attribution-based methods for LLM fine-tuning.</caption><thead><tr><td>References</td><td>Attribution Method</td><td>Target Function</td><td>Granularity</td><td>Selection Strategy</td></tr></thead><tbody><tr><td>Han and Tsvetkov (2021)</td><td>Grad-Dot</td><td>Validation Loss</td><td>Per-example</td><td>Custom &amp; Online</td></tr><tr><td>He et al. (2024)</td><td>Grad-Dot</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Xia et al. (2024)</td><td>TSL00</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Wang et al. (2024g)</td><td>Shapley</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Cao et al. (2024b)</td><td>Regression Datamodel</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Zhou et al. (2024)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>San Joaquin et al. (2024)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Zhao et al. (2024)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Threshold &amp; Offline</td></tr><tr><td>Chhabra et al. (2024)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Custom &amp; Offline</td></tr><tr><td>Jiao et al. (2025)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Online</td></tr></tbody></table>

Table 9 summarizes representative methods applied to LLM fine-tuning. Bolded entries in the **Target Function** column indicate cases where the validation loss is evaluated on a carefully curated evaluation set to reflect task-specific goals, as discussed earlier. Below, we highlight several notable works that illustrate distinctive algorithmic choices and design decisions in this regime.

Computation-wise, most of the referenced methods rely on relatively expensive attribution algorithms compared to those used in LLM pre-training. Among them, Wang et al. (2024g); Cao et al. (2024b) is particularly notable for employing weighted marginal contribution methods, which generally require multiple rounds of re-training. While computational cost is somewhat less prohibitive in the fine-tuning setting, repeated fine-tuning remains impractical at scale. To address this, Cao et al. (2024b) leverages proxy models techniques as we have seen in LLM pre-training; on the other hand, Wang et al. (2024g) approximates Shapley values using the empirical neural tangent kernel (NTK), which avoids explicit retraining by assuming the model remains close to its pre-trained initialization, which is a reasonable assumption in the fine-tuning regime (see Section 2.3).

Beyond attribution itself, the fine-tuning setting also permits more complex training-time integration of selection signals. For example, Han and Tsvetkov (2021) proposes an *online* and *differentiable* selection strategy by directly incorporating attribution scores into the training objective. The model is guided to prioritize samples aligned with a curated validation set, while penalizing those that reinforce spurious correlations. As the underlying attribution scores are estimated via Grad-Dot, this method requires computing higher-order gradients during training, making it computationally intense compared to more straightforward selection strategies. This broader theme of tailoring attribution pipelines to task-specific goals is also reflected in other designs. Chhabra et al. (2024), for instance, applies outlier detection algorithms directly on attribution scores rather than raw input features to identify and discard anomalous training examples.

**Data selection for large vision-language models.** Large vision-language models (LVLMs) have emerged to meet the growing demand for foundational models capable of processing multi-modal inputs, particularly combining visual data with text. A typical training paradigm, following the foundational work of Liu et al. (2023), adopts a *two-stage* process:
*pre-training* followed by *fine-tuning*. This setup assumes modality-specific backbone models, which typically consist of a pre-trained LLM for text and a pre-trained visual encoder. In the pre-training stage, LVLMs focus on aligning representations across modalities, while the fine-tuning stage adapts the model to specific downstream vision-language tasks, much like in LLM fine-tuning.

A (potentially surprising) key challenge in **both** stages is that the training data is typically large-scale and noisy, similar to LLM pre-training. This is especially true when synthetic vision-language datasets are created by pairing raw modality-specific data, where even modest base datasets (e.g., typical fine-tuning datasets) can lead to a combinatorial explosion in the number of examples. Consequently, data selection again emerges as a natural technique to reduce training cost by prioritizing high-utility samples.

<table><caption>Table 10: Comparison of data-attribution-based methods for LVLM pre-training (Top) and fine-tuning (Bottom).</caption><thead><tr><th>Reference</th><th>Attribution Method</th><th>Target Function</th><th>Granularity</th><th>Selection Strategy</th></tr></thead><tbody><tr><td>Zhou et al. (2024)</td><td>IF</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Liu et al. (2024b)</td><td>Grad-Dot</td><td>Validation &amp; Self-Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr></tbody></table>

As LVLMs are still a relatively new research area, only a few works have explicitly explored data selection in this context. We summarize two representative studies in Table 10: Zhou et al. (2024), which focuses on pre-training, and Liu et al. (2024b), which addresses fine-tuning. Both adopt *validation loss* as a core target function. A distinguishing feature of Liu et al. (2024b) lies in its design of the target function: since the fine-tuning task typically spans multiple subtasks, each with potentially different characteristics, the authors augment the standard validation loss with a measure of *task difficulty*. Specifically, they compute the average *self-loss* of samples within each task as a proxy for task difficulty, and combine this with validation loss to guide more nuanced data selection across subtasks.

To mitigate the potentially high computational cost of counterfactual prediction, both works adopt techniques that echo those used in LLM pre-training—particularly *efficient influence approximation* and *reduced frequency of selection*. In terms of attribution algorithms, Zhou et al. (2024) introduces a highly optimized influence function based on the *hyperpower method*, specifically Schulz’s iterative algorithm, while Liu et al. (2024b) employs Grad-Dot as a scalable approximation.

On the selection side, both studies implement a simple top-k selection strategy in an offline setting. Importantly, this choice is not just a pragmatic trade-off between accuracy and efficiency; it is also conceptually motivated. In the case of LVLM pre-training, where the goal is primarily to align features between backbone encoders (usually through training additional adapters by fixing the pre-trained backbone encoders’ weights), and hence the model is not expected to evolve significantly during training. The same happens in the case of LVLM fine-tuning. As a result, attribution scores computed with respect to the initial model can remain valid throughout training, making offline selection a reasonable and effective approach.

**Data selection in the pre-LLM era.** In contrast to recent trends in large language models, earlier work on generative models focused on narrower NLP tasks such as machine translation with smaller model sizes. In these settings, data selection was primarily geared towards *data cleaning*, where identifying outliers and error-prone instances played a crucial role. Robustness was essential because these approaches typically involved end-to-end training of text-generation models. Moreover, the smaller scale of these tasks made repeated re-training feasible, thus justifying and facilitating offline selection.

<table><caption>Table 11: Comparison of data-attribution-based methods for LM training in the pre-LLM era.</caption><thead><tr><th>Reference</th><th>Attribution Method</th><th>Target Function</th><th>Granularity</th><th>Selection Strategy</th></tr></thead><tbody><tr><td>Lam et al. (2022)</td><td>Influence Function</td><td>Validation Loss</td><td>Per-example</td><td>Top-k &amp; Offline</td></tr><tr><td>Bejan et al. (2023)</td><td>Influence Function and TracIn</td><td>Self-Loss</td><td>Per-example</td><td>Custom &amp; Online</td></tr><tr><td>Li et al. (2024a)</td><td>Custom Heuristics</td><td>Likelihood Loss</td><td>Token-level</td><td>Custom &amp; Offline</td></tr></tbody></table>

We summarize three representative works in Table 11, which collectively exhibit a diverse range of techniques and creative problem formulations. Among these, Bejan et al. (2023) and Li et al. (2024a) stand out for their distinctive, custom approaches. Bejan et al. (2023) investigates the role of self-influence in data selection and proposes an online, *learnable* selection mechanism that leverages automated curriculum learning; by integrating with dynamic attribution methods such as TracIn, their approach naturally supports online decision-making that adapts to evolving model behavior. On the other hand, Li et al. (2024a) takes a fine-grained approach by analyzing the distribution of non-target
tokens during text generation. They observe that a high overall loss for a token may arise either from inherently ambiguous next-token distributions or from actual mispredictions (or, equivalently, wrong labels). To differentiate between these cases, they compute an $l_2$ error norm between the predicted distribution and the one-hot target, thereby considering all non-target tokens. The hypothesis is that training on tokens with excessively high $l_2$ error (indicative of “wrong” label) leads to non-robust performance. Based on this insight, their custom selection strategy truncates tokens with error norms above a preset threshold by setting their loss to zero—effectively filtering out noisy or high-entropy training signals rather than discarding full examples.

Overall, these methods reflect the creative and varied approaches developed during the pre-LLM era, offering valuable insights into how data cleaning and selection strategies can enhance the robustness of generative language models.

### 3.1.3 Discussion

Data-attribution-based methods for data selection aim to directly estimate the counterfactual impact of modifying the training set with respect to a target function, offering actionable insights for dataset refinement. While the objective and constraints vary across application settings, several common challenges persist:

* **Additivity assumption.** As discussed earlier, common selection strategies such as top-k and thresholding implicitly assume that attribution scores are additive across training units. That is, the effect of a subset is approximated by summing per-example scores. However, this assumption is not theoretically justified for many popular attribution methods, including influence functions (Hu et al., 2024a) and data Shapley (Wang et al., 2024f). Accurately attributing the collective effects of data subsets remains computationally challenging. Online selection strategies offer a partial remedy by capturing marginal effects dynamically during training.

* **Dependence on a reference model.** A fundamental challenge in attribution-based selection is its post-hoc nature: attribution scores are computed with respect to a reference model, typically trained on the full raw dataset $D_{raw}$. This presents a dilemma—training on $D_{raw}$ to enable attribution defeats the purpose of data selection in resource-constrained settings. Two common workarounds have been explored: one is using lightweight proxy models as stand-ins for the full model, and another is to adopt online attribution methods that estimate influence during training without requiring a fully trained model upfront.

* **Sensitivity and robustness.** Since data selection aims to minimize human intervention and operate autonomously, the reliability of attribution scores is paramount. If these scores are noisy or biased, there is no clear mechanism for validation or correction. Attribution methods are known to be sensitive to input perturbations (Koh and Liang, 2017) and vulnerable to adversarial manipulation (Wang et al., 2025c). Despite their growing use, robust and trustworthy attribution algorithms remain an underexplored research direction.

Despite these challenges, data-attribution-based methods offer several key advantages over alternative data selection approaches. Unlike classifier-based filtering (Raffel et al., 2020; Longpre et al., 2024), which often lacks transparency, or heuristic filtering (Gan et al., 2024; Albalak et al., 2024), which can introduce unintended biases or degrade performance, attribution-based selection provides a more principled and interpretable foundation for identifying impactful training data. For a broader overview of alternative data selection strategies, we refer readers to the survey by Albalak et al. (2024).

In summary, data-attribution-based methods for data selection offer a flexible and principled framework that can accommodate a wide range of desiderata in trustworthy machine learning, simply by adjusting the choice of target function. Compared to alternative approaches, attribution-based methods provide greater headroom for improvement, particularly when the benefits of heuristic-based selection begin to saturate (Wang et al., 2025a). As such, they can serve as a *complementary* component in broader data selection pipelines. A promising hybrid strategy involves using classical filtering techniques for coarse pruning, followed by data-attribution-based methods for more fine-grained selection and cleaning. This combination balances computational efficiency with selection precision, making it especially well-suited for large-scale generative AI applications.

We also note that while current work on data-attribution-based selection has primarily focused on LLMs, other generative settings, such as image generation via diffusion models, remain largely unexplored. Although recent efforts have advanced data attribution techniques for diffusion models (Lin et al., 2025; Xie et al., 2024a), their integration into data selection pipelines has yet to be studied. Extending attribution-based selection to vision generative models thus represents a promising direction for future research.

## 3.2 Fact Tracing and Concept Probing

Modern generative models have made remarkable progress in both text and image generation. However, their complexity and opaque generation processes pose significant challenges in understanding how the content is generated. A key