Scalable Multi-Stage Influence Function for Large Language Models via Eigenvalue-Corrected Kronecker-Factored Parameterization

1 Introduction

Understanding the relationship between training data and model behavior is essential for building trustworthy machine learning systems. For example, attributing a model's answer in a closed-book question-answering system to a specific Wikipedia article can enhance user trust. Training data attribution (TDA) techniques quantify the contributions of training instances to model behaviors by addressing a counterfactual question: how would the model's behavior change if an example were removed from the training set? Originating from robust statistics [Hampel, 1974] and introduced into deep learning by Koh and Liang [2017], influence function

∗Corresponding author

(IF) provides an end-to-end, scalar-valued interpretation of a model's high-level behavior.

While several works on IFs have explored their conceptual framework and applications, they often overlook the scalability of these methods to large-scale neural networks trained on massive datasets. Furthermore, prior analyses predominantly focus on classical architectures, such as feed-forward networks, rather than modern architectures like transformers. Findings of Zhou et al. [2024] indicate that most knowledge in large language models (LLMs) is acquired during pretraining. Consequently, explaining the predictions of finetuned models necessitates tracing influence back to the pretraining dataset rather than the fine-tuning dataset. Although Chen et al. [2020] introduced "multi-stage" IFs to address this need, their analysis of NLP models is limited to frozen encoders stacked with linear classifiers, failing to accommodate the full-parameter tuning paradigm prevalent in LLMs.

In this work, we propose an IF-based method to estimate the contribution of pre-training data to the predictions of finetuned models, scaling efficiently to LLMs with billions of parameters. Our approach addresses two key challenges.

First, the original IF framework [Koh and Liang, 2017] does not accommodate multi-stage influence when finetuning tasks require substantial modifications to the model, such as replacing the unembedding layer. Inspired by Chen et al. [2020], we extend IFs to the multi-stage paradigm ("pretrain then fine-tune"), enabling attribution of downstream predictions of LLMs to pre-training examples.

Second, scaling IFs to LLMs involves overcoming computational bottlenecks related to inverse Hessian-Vector Product (iHVP) and training gradients. For iHVPs, we adopt the Eigenvalue-corrected Kronecker-factored Approximate Curvature (EK-FAC) parameterization [George et al., 2018], as suggested by Grosse et al. [2023]. To address the latter bottleneck, we leverage semantic-similarity-based heuristics to narrow the candidate training samples, avoiding iterations over the entire dataset.

We conduct extensive experiments to evaluate our method. For general influence estimation, we demonstrate the superior scalability of EK-FAC approximations compared to various TDA methods. For our multi-stage IF, we evaluate it on our fact-tracing benchmark and show that it outperforms the single-stage IF [Grosse et al., 2023]. Additionally, we analyze the contributions of MLP and multi-head attention (MHA) parameters to influence estimation, finding that MLP parameters play a proportionally greater role. This insight suggests a practical trade-off, allowing analysis to focus on MLP parameters in large models. Finally, we apply our multistage IF on a publicly available instruction-tuned LLM, dollyv2-3b [Conover et al., 2023], qualitatively explaining its generations based on pre-training data.

In summary, we propose a general framework for efficiently estimating multi-stage influence with the help of EKFAC parameterization and semantic-similarity-based candidate selection. Our results provide practical insights into resolving influence estimation trade-offs. Furthermore, we demonstrate how TDA approaches can be applied to calibrate the trustworthiness of generative AI systems.

2 RelatedWork

Training data attribution. Training data attribution (TDA) techniques explain a model's predictions by quantifying the contribution of training data. As noted by Hammoudeh and Lowd [2024], TDA methods can be broadly categorized into retraining-based and gradient-based approaches. Retrainingbased methods estimate the influence of individual examples by retraining models on random subsets of the training dataset. However, these methods incur high computational costs due to the need for multiple retraining rounds, rendering them impractical for large-scale models and datasets. Gradient-based methods, which infer training data influence using gradients, are further divided into dynamic and static approaches. Dynamic estimators, such as TracIn [Pruthi et al., 2020], assess influence by analyzing gradients from intermediate model snapshots captured during training. In contrast, static estimators, including IFs [Koh and Liang, 2017] and representer point selection [Yeh et al., 2018], rely solely on the final model parameters to compute influence.

Influence functions. Despite their utility, IFs exhibit several limitations. First, in terms of model architecture, most existing studies focus on traditional architectures such as feed-forward and recurrent networks, with limited exploration of modern architectures like transformers. A recent study by Grosse et al. [2023] extended IF analysis to transformer-based LLMs using EK-FAC approximation. In this work, we also investigate transformer language models.

Second, regarding scalability, most research adopts a "matrix-free" approach to avoid the computational costs of explicitly and inverting the Hessian for large models, but with limited success. One prominent strategy involves parametric approximations of the Hessian to enable efficient inverse Hessian computations [Schioppa et al., 2022; Grosse et al., 2023]. Another approach uses iterative stochastic approximation methods, such as LiSSA and Conjugate Gradient, to approximate iHVPs, as introduced by Koh and Liang [2017]. However, our experiments show that these methods fail to yield usable influence estimates within a practical timeframe.

Third, on the training paradigm, the original IF proposed by Koh and Liang [2017] is unsuitable for analyzing the influence of pre-training data on a fine-tuned model when it has a different output domain. While Grosse et al. [2023] proposed an efficient method for single-stage training scenarios, it does not address multi-stage paradigms. Chen et al. [2020] introduced multi-stage IF as a generalization of the original IF, but their work is limited to the ELMo [Peters et al., 2018] architecture and fine-tuning scenarios where a classification head is added to a frozen pre-trained encoder. In contrast, we focus on the popular full-parameter fine-tuning paradigm.

3 Background and Notations

3.1 Transformer Architecture

The focus of our work is the transformer language model [Vaswani et al., 2017], which starts with a token embedding, followed by a series of L residual blocks, and ends with a token unembeddings [Elhage et al., 2021] (layer normalization is omitted for brevity). For a sequence t with n tokens, the initial embeddings are x0 = Embed(t) ∈ Rn×d , hence the start of the residual stream.

At each residual block, the residual stream is first processed by the MHA module: ˜xl = xl-1 + ∑h∈Hl h(xl-1), where Hl is the set of attention heads and xl-1 is the input for the l-th (1 ≤ l ≤ L) residual block. The attention pattern for head h is obtained via attention mechanism: rh = Attention (qh , kh , vh ), where qh = xl-1 (WhQ)⊤, kh = xl-1 (WhK )⊤, vh = xl-1 (Wh V )⊤. WhQ,WhK ,Wh V ∈ RnHcontext×d are weights for query, key and value, respectively. The attention pattern is then written into the residual stream: h(xl-1) = rh (WO)⊤, WO ∈ Rd×nHcontext .

The MLP module then processes the output of the MHA and writes back to the residual stream: xl = ˜xl + σ (˜xl (Wm I )⊤) (Wm O )⊤, where σ(·) is the element-wise nonlinear activation, Wm I ∈ Rnmcontext×h andWm O ∈ Rh×nmcontext are projection weights. After L residual blocks, the unembeddings produce the final logits: T(t) = xLW⊤ U .

3.2 Influence Functions

In this section, we demonstrate the original single-stage IF [Koh and Liang, 2017]. This formulation quantifies the influence of a training example on both model parameters and a measurement of a query sample.

Consider a training dataset D = {zi}Ni=1 , where each example zi represents a token sequence, and the learning task is self-supervised language modeling. The training objective is to minimize the expectation L of the loss function ℓ(·):

The influence of a training sample z on the optimal parameters θ∗ is Iθ∗(zm), while the influence on a query zq is mediated via a differentiable measure, m(zq, θ), e.g., autoregressive cross-entropy loss of zq. The influence of z with respect to θ∗ can thus be expressed as:

In practice, the Hessian may be singular when the model is not fully converged. To address this, we follow Bae et al. [2022] and replace the Hessian with a damped generalized Gauss-Newton (GGN) matrix [Schraudolph, 2002; Martens, 2020] to ensure positive-definiteness. This modification enables the use of the final model parameters in place of strictly converged parameters θ∗ . Additionally, we use the cross-entropy loss to ensure that the loss function is convex with respect to model outputs.

3.3 EK-FAC for Feed-Forward Networks

Naive computation of the GGN (G) or its inverse is computationally prohibitive. To address this, EK-FAC [George et al., 2018] was proposed as an efficient approximation method for computing iHVPs. This approach was subsequently adopted by Grosse et al. [2023] to estimate single-stage influence for LLMs. Structurally, the GGN is approximated as a blockdiagonal matrix. Each diagonal block is then approximated using EK-FAC parameterization.

Consider a feed-forward neural network with L layers interleaved with nonlinear activations. Let the parameters be W = (W1,W2, ...,WL) and θ = (θ⊤1 , ..., θ⊤L )⊤, where θl = vec(Wl) denotes the vectorized form of Wl(l = 1, ..., L). The GGN for this network is expressed as:

where ˆQx is the training dataset, Dθ = ∇θ log p(y|fθ(xn), θ), and Gi,j = E [DθiDθ⊤j ] is a block of the GGN. Note that the label y is sampled from the model's predictive distribution rather than the training label.

Approximating the GGN as block-diagonal retains only the diagonal blocks: G ≈ ˜G = diag(G1,1, ...,GL,L). Given a vector v with the same dimensions as θ, the damped iHVP can be computed as ( ˜G+λI)-1 v = diag ((G1,1 + λI)-1 v1, ..., (GL,L + λI)-1 vL), where vl is the slice of v corresponding to θl.

For each block Gl,l, EK-FAC provides a further approximation. Let the inputs and outputs of the l-th layer be al-1 and sl, respectively. The gradient Dθl is expressed as al-1 ⊗ Dsl, where ⊗ denotes Kronecker product. Gl,l can thus be approximated using K-FAC as:

EK-FAC improves upon K-FAC by incorporating the diagonal variance in the eigenbasis of Al-1,l-1 and Sl,l. Denoting these matrices as A and S for brevity, their eigendecompsition yields:

where QA and QS are eigenvectors, and ΛA and ΛS are diagonal matrices of eigenvalues. To account for the diagonal variance, the middle factor is replaced by a new diagonal matrix Λ with its diagonal entries as follows:

Finally, the damped iHVP for the l-th block is computed as:

where Λλ = Λ + λI(λ > 0), ⊘ denotes element-wise division, diag-1 (·) extracts the diagonal elements of a matrix into a vector, and unvec(·) converts a vector into a matrix.

4 Method

In this section, we present the design of our multi-stage IF, which attributes the predictions of a fine-tuned LLM to its pre-training data. Additionally, we describe the use of approximation techniques to make the multi-stage IF computationally tractable for LLMs, addressing the trade-off between efficiency and effectiveness.

The objective of the multi-stage IF is to quantify the influence of a pre-training sample z on a test-time query sample xt. The estimation process consists of two primary steps: candidate selection (Section 4.3) and influence computation (Section 4.2). First, given a query, we filter the pre-training dataset using similarity heuristics for a much smaller subset of training examples as candidates for influence estimation. This step ensures that the selection is focused on training examples that are semantically relevant to the query. In the second step , we compute the influence of the selected candidates on the query, producing a series of influence scores.

4.1 Multi-Stage Influence Function

Before formulating the multi-stage IF, we first review the "pre-train then fine-tune" paradigm. During pre-training, all model parameters are randomly initialized and subsequently updated to fit the distribution of a large-scale corpus.

where ℓpt(·) is the pre-training loss.

During fine-tuning, the model parameters are initialized with θpt and subsequently optimized on a fine-tuning dataset under the cost function ℓft(·).

A naive approach to derive multi-stage influence is as follows:

This formulation is valid as long as the fine-tuned model shares the same output domain as the pre-trained model. For example, it works for pre-trained language models fine-tuned to follow instructions. However, discrepancies may arise when the final unembedding layer is replaced during finetuning. For instance, pre-training typically involves the language modeling task with outputs as logits over a large vocabulary, fθpt (z) ∈ R|V| , whereas fine-tuning may adapt the model for binary sequence classification, resulting in outputs fθft (x) ∈ R2 . Consequently, the gradient ∇θℓpt(z, θft ) is formally invalid for the fine-tuned model.

To address this, we propose making the pre-training gradient accessible to the fine-tuned model by establishing a connection between their parameter spaces. Fine-tuning implicitly assumes that the model retains its general capabilities after adaptation. Thus it is reasonable to assume that the finetuned parameters do not deviate significantly from pre-trained ones. Under this assumption, fine-tuning can be viewed as minimizing empirical risk on the fine-tuning task in the vicinity of the pre-trained parameters in the parameter space.

Inspired by Chen et al. [2020], we instantiate a quantifiable connection between the pre-trained model and its fine-tuned successor by introducing an additional proximity constraint to the fine-tuning objective in a post-hoc manner. Specifically, we define the proximity as the Euclidean distance between fine-tuned and pre-trained parameters:

where α ∈ R+ is a hyperparameter and || · ||2 is L2 norm.

Based on this reformulated objective above, the influence of a pre-training example z on a test instance xt under the measurement m(·) is given by (full proof in Appendix):

4.2 Practical Implementation for LLMs

A straightforward approach to implementing IFs involves calculating Hθ and H-1 θ , followed by computing iHVPs via H-1 θ v. However, for a model with p parameters and N training samples, O(Np2 + p3 ) operations are required [Koh and Liang, 2017], which is computationally prohibitive when p and N are large. To address this, following Grosse et al. [2023], who approximated single-stage IF with EK-FAC, we adopt a similar approach to approximate the multi-stage IF as described in Section 3.3.

First, we replace Hessians with damped GGNs: ∇2θ Lft (θft) ≈ Gft + λftI and ∇2θ Lpt (θpt ) ≈ Gpt + λptI. For Gft, the term α in Equation (12) is absorbed into the damping term λft. Consequently, the multi-stage IF (Equation (12)) can be interpreted as the inner product of two preconditioned gradients: (Gpt +λptI)-1 ∇θℓpt (z, θpt ) and (Gft +λftI)-1 ∇θm (xt, θft). In practice, we implement multi-stage IF following this interpretation.

Next, we compute EK-FAC factors for Gpt and Gft, focusing on the linear components of the model. The EK-FAC factors are precomputed and stored on disk, to be loaded when needed. These includeWm I andWm O of MLP modules, WhQ, WhK andWh V in each attention head, andWO of MHAmodules. We exclude unembedding parameters, as the unembedding layer of θpt and θft may have different output dimensions. We also exclude layer normalization modules, since their parameter count is maginal and they are usually not considered to encode factual knowledge [Grosse et al., 2023].

We focus on autoregressive decoder-only LLMs, the pretraining loss is the cross entropy loss, following Grosse et al. [2023]. For a sequence z with T tokens:

where pˆ y|x is the pre-trained model's output distribution.

Computational and Spatial Cost

The one-time cost of preparing EK-FAC factors is considered an overhead amortized across future influence analyses. During influence score computation, for a weight matrix with input dimension d and output dimension p, the cost of computing an iHVP is O(d2p + dp2 ), followed by an inner product between query iHVP and candidate iHVP at O(dp). The memory and storage overhead arises from storing eigenvectors Q and the diagonal entries of Λ, resulting in an extra spatial cost of d2 + p2 + dp.

For example, in the GPT-NeoX [Andonian et al., 2023] architecture, where d = p for query, key, value, and output linear projection weights in MHA modules, the additional spatial cost is 3 times the size of the original weights. For MLP modules, where d = 4p or d = 1 4 p, the additional spatial cost is 5.25 times the size of the original weights.

4.3 Selecting Candidates for Influence Estimation

To identify a few positively influential training examples for a query, a naive approach would compute the influence of every training example on the query. However, this requires gradient computations across the entire dataset, a cost equivalent to one epoch of training. To address this, we narrow down the candidate set using efficient similarity-based heuristics inspired by Grosse et al. [2023].

While Grosse et al. [2023] employed TF-IDF [Ramos, 2003], we adopt an unsupervised K-Nearest Neighbors (KNN) approach based on the embeddings of pre-training documents, similar to the approach of Guo et al. [2021]. The embeddings are generated using Sentence Transformers [Reimers and Gurevych, 2019], and a KNN classifier is constructed over these embeddings. This choice reflects the intuition that training examples with similar semantics to the query are more interpretable and relevant than those based solely on textual overlap [Karpukhin et al., 2020]. Moreover, this method aligns with the potential application of multistage IFs in identifying pre-training documents that serve as grounding knowledge sources.

5 Experiments

This section presents a series of experiments to evaluate our proposed method. First, we assess the scalability of singlestage IFs approximated using EK-FAC compared to various TDA methods in terms of both estimation accuracy and wall-clock runtime on language modeling tasks (Section 5.1). Next, we investigate the validity of the underlying assumption of our multi-stage IF (Section 5.2). Subsequently, we evaluate the effectiveness of the multi-stage IF on a factual knowledge retrieval task (Section 5.3). Finally, we showcase a qualitative case study using the multi-stage IF on an instructionfollowing LLM (Section 5.4). Our code is public at https: //github.com/colored-dye/multi stage influence function.

5.1 Scalability Validation of EK-FAC

This experiment evaluates the effectiveness and efficiency of EK-FAC parameterization in producing influence estimates. For simplicity, we focus exclusively on single-stage language modeling rather than the "pre-train then fine-tune" scenario.

Dataset and Model. We train a custom GPT-NeoX model on The Penn Treebank dataset [Marcus et al., 1993]. The model consists of decoder layers, two attention heads, and a hidden size of 256, comprising approximately 1.18M MLP parameters and 0.39M MHA parameters.

Baselines. We evaluate various TDAmethods as baselines. For influence estimation methods relying on iterative iHVP estimation, we include Conjugate Gradient (CG) and LiSSA following Koh and Liang [2017]. To reduce computational costs, we use CG as the ground truth instead of methods like training under the proximal Bregman response objective [Bae et al., 2022] or linear datamodeling scores [Park et al., 2023].

For IF estimation methods based on estimating diagonal entries of the Hessian, we employ IF with Arnoldi iteration [Schioppa et al., 2022].

TRAK [Park et al., 2023], a retraining-based baseline, uses a set of models trained on random subsets. We also choose Gradient dot product (GDP) [Charpiat et al., 2019; Grosse et al., 2023] and linear Centered Kernel Alignment (CKA) [Kornblith et al., 2019] as gradient-similarity-based baselines. GDP simply computes the dot product between query gradients and candidate gradients. CKA is specifically designed for measuring representation similarity. Here, we use it on query and training gradients.

Metrics. We randomly select ten test samples, each paired with 500 candidates sampled from the training split. The primary metric for influence estimation accuracy is the Spearman correlation coefficient (ρ), as ranking quality is our main focus. "Pair-wise runtime" refers to the cost of computing the influence of 500 candidates on a query. "Overhead runtimes" are specific to each baseline: fitting EK-FAC factors for EK-FAC, calculating dominant eigenpairs for Arnoldi, or training and gradient featurization for TRAK. Spearman correlations and pair-wise runtimes are averaged over ten trials, while overhead runtime is measured once per baseline.

Results. Figure 1 visualizes results, with the best results shown in Table 1. Ideally, a method should achieve highquality estimates within minimal runtime, corresponding to markers closer to the upper left corner of the figure. Key findings include:

1. EK-FAC achieves the best trade-off between approximation quality and computation cost, with its marker closest to the upper left corner. Furthermore, it resides on a more optimal Pareto frontier than CG and LiSSA.
2. Ablating influence analysis for MHA parameters does not result in substantial degradation in approximation quality compared to analyzing both MHA and MLP parameters. Despite MHA parameters accounting for 25% of analyzed parameters, they contribute only 14.5% of the total influence, indicating that MLP parameters have a more significant impact.
3. The approximation quality of CG and LiSSA scales loglinearly with computation time, with CG offering a superior Pareto frontier.

(a) BLOOM-560m v.s. BLOOMZ-560m.

(b) Pythia-2.8b v.s. dolly-v2-3b.

Figure 2: Distribution of ||θft - θpt ||2 and ||θft - θpt ||2/||θpt ||2.

1. Gradient similarity-based methods (GDP and CKA) are the most efficient baselines but yield low-quality influence estimates.
2. TRAK, a representative of retraining-based methods, has the lowest estimation quality and the highest computational cost.
3. Arnoldi yields poor estimates at ∼120× the pair-wise compute cost of EK-FAC.

Observation (2) highlights the potential for further approximations in influence estimation. As noted by Grosse et al. [2023], factual associations are primarily localized within MLP modules [Meng et al., 2022], thus MLP modules are also likely to contribute significantly to influence estimation. Our findings support this claim, suggesting that for largescale models, focusing solely on MLP parameters can significantly reduce computational costs while maintaining useful influence estimates.

5.2 Euclidean Proximity in Practice

As described in Section 4.1, the multi-stage IF relies on the assumption that the fine-tuned parameters are geometrically close to its pre-trained predecessor in the parameter space. To validate this assumption, we conduct an experiment analyzing two pairs of models: BLOOM-560m versus BLOOMZ560m, and Pythia-2.8b versus dolly-v2-3b. We focus on the linear weights of both the MLP and MHA modules, ignoring bias terms, as the bias parameters only constitute a minor portion of the total parameters.

The results, presented in Figure 2, indicate that the distance between fine-tuned and pre-trained weights (||θft -θpt ||2) accounts for no more than 8% of the L2 norm of the pre-trained weights, ||θpt ||2. Specifically, for the BLOOM-560m and BLOOMZ-560m pair, ||θft - θpt ||22 = 368.7, while for the Pythia-2.8b and dolly-v2-3b pair, ||θft - θpt ||22 = 94.5.

In practice, the damping term of the GGNs is typically very small; in subsequent experiments, we set α ≤ λft = 10-4 .

For the two model pairs, the additional Euclidean proximity term introduces an extra fine-tuning loss of 0.0184 and 0.0047, respectively. These additional losses are negligible compared to the final training losses of BLOOMZ-560m and dolly-v2-3b, which are 1.6132 and 0.8209, respectively.

5.3 Effectiveness of Multi-Stage Influence Function

In this experiment, we evaluate the effectiveness of our multistage IF in identifying pre-training data that significantly influence predictions of a fine-tuned model. To facilitate this evaluation, we establish a benchmark with ground-truth labels for attributing test-set instances to corresponding samples in an attribution set. To enable comparison with the single-stage IF, the chosen task does not require modifying the unembedding layer.

Benchmark Setup. We construct a fact-tracing benchmark based on that of Aky¨urek et al. [2022], consisting of an attribution set and a test set. Each test set instance corresponds to a fact and is associated with several sentences in the attribution set. The evaluated methods are tasked with correctly retrieving the knowledge sources for each fact.

The attribution set is derived from the T-REx dataset [Elsahar et al., 2018], which aligns knowledge base triples with DBpedia abstracts. Instead of directly using the original dataset, we create a sentence-level subset of T-REx that includes all knowledge base triples represented in the LAMA dataset [Petroni et al., 2019].

The test set is derived from the T-REx split of LAMA. Aky¨urek et al. [2022] use the test set as cloze-style language modeling samples, as their target models are mT5based [Xue, 2020] models, which are trained under the masked language modeling "span-corruption" objective. But our target model is an instruction-following autoregressive LLM, so we manually convert the original text completion templates into equivalent question answering formats. For example, a knowledge base triple (X, born in, Y) originally uses the template "X was born in," which is a text-completion task, while our version reformulates it as "Where was X born?" to align with a question-answering task, where the model is required to predict the object Y.

Metrics. We evaluate fact retrieval performance using standard information retrieval metrics, including Mean Reciprocal Rank (MRR) and Recall@10. The MRR is defined as 1Q ∑q∈Q 1 rankq , where Q is the test set and rankq denotes the rank of the first correct knowledge source for query q. Results are averaged over three trials, with each trial sampling 200 test instances.

Baselines. We compare the proposed multi-stage IF (MSIF) against several baseline methods. BM25 [Robertson et al., 1995], a model-agnostic baseline, retrieves relevant samples based on word-level overlap. Additionally, we include two similarity-based methods: Representation Similarity (RepSim) [Caruana et al., 1999], which computes cosine similarity between the hidden states of pre-trained and fine-tuned models, and GDP, a simplified multi-stage IF where the Hessian matrices are identity matrices.

The single-stage IF (SS-IF) is also evaluated, using only the fine-tuned model. Both MS-IF and SS-IF adopt the same measurement m(x, y; θ) = -log pθ(y|x), where x and y represent the question and answer tokens, respectively. Both use EK-FAC approximation and analyze both the MLP and MHA modules. Furthermore, as both IFs rely on damping terms, we report their performance under different damping values, ranging from 10-8 to 10-4 .

Models. We utilize BLOOM-560m [Workshop et al., 2022] and BLOOMZ-560m [Muennighoff et al., 2022], the latter of which has undergone multitask fine-tuning. We choose these models as we are able to verify that their pre-training data encompass most of the knowledge of the fact-tracing benchmark. Specifically, as each knowledge instance is represented as (subject, relation, object) triples, we confirm that 96.81% of the objects are included by the pre-training data of BLOOM(Z)-560m.

Results. As is shown in Figure 3, the multi-stage IF outperforms all baseline methods in terms of both MRR and Recall under smaller damping terms (10-6 and 10-8 ), demonstrating its superior ability to assign higher influence scores to ground-truth knowledge sources. However, the performance gaps between MS-IF, SS-IF, and GDP are relatively small under a large damping term (10-4 ). This indicates that the choice of the damping term is essential to the effectiveness of IFs, and that both SS-IF and MS-IF degenerate to GDP when the damping terms are large. Moreover, even the bestperforming method yields results that are far from perfect retrieval. This discrepancy may stem from approximation errors, the design of our multi-stage IF, or limitations in the language models' suitability for knowledge retrieval tasks.

5.4 Case Study

In addition to the quantitative experiments above, we qualitatively demonstrate the interpretive power of the proposed multi-stage IF through a case study. The analysis is conducted on dolly-v2-3b for the task of factual knowledge attribution. The motivation for this case study stems from possible user concerns regarding whether an LLM-based interactive system's responses are grounded in reliable knowledge sources or are mere hallucinations. To address this, it is reasonable to attribute model-generated outputs to the pre-training data. By inspecting whether the top-influential pre-training texts contain relevant and accurate information, we are able to assess whether the model's responses are properly grounded.

For the factual knowledge attribution task, the model is prompted with a question, and a response is generated using greedy decoding. The article with the highest multi-stage influence is identified from the Wikipedia subset of the pretraining corpus using the pipeline described in Section 4. Table 2 presents the results, showing an excerpt of the retrieved article. While the response itself is incorrect, the retrieved article is highly relevant to the queried topic.

6 Conclusions and Limitations

This paper introduces a generalization of the IF to enable the attribution of predictions made by a fine-tuned model to its pre-training data. To enhance the scalability of influence computation, we employ EK-FAC parameterization and a nearest-neighbor-based candidate selection strategy. Experimental results confirm the effectiveness and efficiency of the proposed multi-stage IF, demonstrating its applicability to LLMs with three billion parameters.

While our work enhances the scalability and generality of influence functions, several limitations remain. First, we only analyze backbone MLP and MHA components, excluding contributions from other components such as embeddings, unembeddings, and layer normalization. Extending influence analysis to these components may improve the quality of influence estimates. However, the potential benefit is likely marginal, as they constitute a small proportion of the overall parameters and are not considered to encode knowledge.

Second, our analyses are limited to decoder-only transformer architectures. Extending scalable influence analysis to other architectures, such as encoder-decoder models or diffusion models, could unlock valuable new applications.

Finally, SOURCE is proposed by Bae et al. [2024] as an effective TDA approach, leveraging approximate unrolled differentiation and inherently suited for multi-stage scenario. Their results demonstrate the superior performance of SOURCE compared to influence functions on models like BERT [Devlin, 2018] and GPT-2 [Radford et al., 2019]. Unfortunately, due to the unavailability of its implementation, we were unable to directly compare our multi-stage influence function with SOURCE in this work.

Acknowledgments

This work was partly supported by the National Key Research and Development Program of China under No. 2024YFB3900105, NSFC under No. 62402418, Zhejiang Province's 2025 "Leading Goose + X" Science and Technology Plan under grant No.2025C02034, the Key R&D Program of Ningbo under No. 2024Z115, and the Open Project of Key Laboratory of General Quality Technology and Application of Intelligent Manufacturing Equipment, Ministry of Industry and Information Technology (HK202403532).