Bayesian Influence Functions for Hessian-Free Data Attribution

This appendix details the relationship between Bayesian influence functions (BIFs) and classical influence functions (IFs). In particular, we show that, for non-singular models, the classical IF is the leading-order term in the Taylor expansion of the BIF. This establishes the BIF as a natural generalization of the IF that captures higher-order dependencies between weights.

Let ${\bm{w}}^{*}$ be a local minimum. In this section, all gradients and Hessians are evaluated at ${\bm{w}}^{*}$; thus, to reduce notational clutter, we omit the dependence on ${\bm{w}}$. For any function $f({\bm{w}})$, we denote its gradient at ${\bm{w}}^{*}$ as ${\bm{g}}_{f}=\nabla_{\bm{w}}f({\bm{w}}^{*})$ and its Hessian as ${\bm{H}}_{f}=\nabla_{\bm{w}}^{2}f({\bm{w}}^{*})$. In particular, ${\bm{g}}_{\phi}=\nabla_{\bm{w}}\phi({\bm{w}}^{*})$ and ${\bm{H}}_{\phi}=\nabla_{\bm{w}}^{2}\phi({\bm{w}}^{*})$ for an observable $\phi({\bm{w}})$; we also abbreviate ${\bm{g}}_{i}=\nabla_{\bm{w}}\ell_{i}({\bm{w}}^{*})$ and ${\bm{H}}_{i}=\nabla_{\bm{w}}^{2}\ell_{i}({\bm{w}}^{*})$ for a per-sample loss $\ell_{i}({\bm{w}})$. The total Hessian of the empirical risk ${L_{\text{train}}}({\bm{w}})=\sum_{k=1}^{n}\ell_{k}({\bm{w}})$ at ${\bm{w}}^{*}$ is denoted ${\bm{H}}=\sum_{k=1}^{n}{\bm{H}}_{k}$.

The Bayesian influence function (BIF) for the effect of sample ${\mathbf{z}}_{i}$ on an observable $\phi$ is given by (see Equation 4):

where the covariance is taken over the posterior $p({\bm{w}}\mid\mathcal{D}_{\text{train}})\propto\exp(-{L_{\text{train}}}({\bm{w}}))\varphi({\bm{w}})$, with $\varphi({\bm{w}})$ being a prior. This definition is exact and makes no assumptions about the form of $\phi({\bm{w}})$, $\ell_{i}({\bm{w}})$, or $p({\bm{w}}\mid\mathcal{D}_{\text{train}})$.

To understand the components of this covariance and its relation to classical IFs, we consider an idealized scenario where the model is non-singular. Under this strong assumption, which does not hold for deep neural networks (Wei et al., 2023), the posterior $p({\bm{w}}\mid\mathcal{D}_{\text{train}})$ can be approximated by a Laplace approximation around ${\bm{w}}^{*}$:

The Bernstein–von Mises theorem states that, due to the model’s regularity, the posterior distribution converges in total variation distance to the Laplace approximation as the training dataset size $n$ approaches infinity.

Let $\Delta{\bm{w}}={\bm{w}}-{\bm{w}}^{*}$. Assuming analyticity, we can express $\phi({\bm{w}})$ and $\ell_{i}({\bm{w}})$ using their full Taylor series expansions around ${\bm{w}}^{*}$:

where $D^{k}f({\bm{w}}^{*})[\Delta{\bm{w}},\dots,\Delta{\bm{w}}]$ denotes the $k$-th order differential of $f$ at ${\bm{w}}^{*}$ applied to $k$ copies of $\Delta{\bm{w}}$.

The covariance under this Gaussian (Laplace) approximation, denoted $\mathrm{Cov}_{\mathcal{N}}$, then involves covariances between all pairs of terms from these two expansions:

where $\text{Term}_{k}[f]$ is the $k$-th order term in the Taylor expansion of $f({\bm{w}})$ in powers of $\Delta{\bm{w}}$. For $\Delta{\bm{w}}\sim\mathcal{N}(0,{\bm{H}}^{-1})$, the leading terms are:

• Covariance of linear terms ($k=1,m=1$):

  $$\mathrm{Cov}_{\mathcal{N}}({\bm{g}}_{\phi}^{T}\Delta{\bm{w}},{\bm{g}}_{i}^{\top}\Delta{\bm{w}})={\bm{g}}_{\phi}^{\top}{\bm{H}}^{-1}{\bm{g}}_{i}.$$

• Covariance of quadratic terms ($k=2,m=2$):

  $$\mathrm{Cov}_{\mathcal{N}}\left(\frac{1}{2}(\Delta{\bm{w}})^{\top}{\bm{H}}_{\phi}\Delta{\bm{w}},\frac{1}{2}\Delta{\bm{w}}^{\top}{\bm{H}}_{i}\Delta{\bm{w}}\right)=\frac{1}{2}\operatorname{tr}({\bm{H}}_{\phi}{\bm{H}}^{-1}{\bm{H}}_{i}{\bm{H}}^{-1}).$$

  (Using Isserlis’ theorem for moments of Gaussians).

• Cross-terms between odd and even order terms (e.g., $k=1,m=2$) are zero due to the symmetry of Gaussian moments.

Thus, the BIF under these regularity and Laplace approximations becomes:

The leading term $-{\bm{g}}_{\phi}^{\top}{\bm{H}}^{-1}{\bm{g}}_{i}=-\nabla_{\bm{w}}\phi({\bm{w}}^{*})^{\top}{\bm{H}}_{{\bm{w}}^{*}}^{-1}\nabla_{\bm{w}}\ell_{i}({\bm{w}}^{*})$ is precisely the classical influence function $\operatorname{IF}({\mathbf{z}}_{i},\phi)$ from Equation 2. Note that ${\bm{H}}$ scales linearly in $n$, so this term dominates as $n\to\infty$. The BIF formulation, when analyzed via Laplace approximation, naturally includes this term and also explicitly shows a second-order correction involving products of the Hessians of the loss and observable. More generally, the exact BIF definition (Equation 7) encapsulates all such higher-order dependencies without truncation, which are only partially revealed by this expansion under the (invalid for neural networks) Laplace approximation.

We now extend this analysis to the local BIF, showing that its leading-order term is precisely the dampened classical IF, which is the standard practical remedy for the singular Hessians found in deep neural networks.

The local BIF is defined over the localized posterior from Equation 5:

This distribution is centered around ${\bm{w}}^{*}$ due to the localizing potential (the quadratic term). To apply the Laplace approximation, we consider the mode of this distribution, which is the minimum of the effective potential $L_{\text{eff}}({\bm{w}})=L_{\text{train}}({\bm{w}})+\frac{\gamma}{2}\lVert{\bm{w}}-{\bm{w}}^{*}\rVert_{2}^{2}$. We assume ${\bm{w}}^{*}$ to be a local minimum of $L_{\text{train}}({\bm{w}})$, so $\nabla L_{\text{train}}({\bm{w}}^{*})=0$. Consequently, $\nabla L_{\text{eff}}({\bm{w}}^{*})=\nabla L_{\text{train}}({\bm{w}}^{*})+\gamma({\bm{w}}^{*}-{\bm{w}}^{*})=0$, meaning ${\bm{w}}^{*}$ is also the mode of the localized posterior.

The precision of the Laplace approximation is given by the Hessian of this effective potential evaluated at ${\bm{w}}^{*}$:

Therefore, the Laplace approximation for the localized posterior is a Gaussian centered at ${\bm{w}}^{*}$ with covariance ${\bm{H}}_{\text{eff}}^{-1}$:

Following the same Taylor expansion logic as in the previous section, we can compute the leading-order term of the covariance between $\phi({\bm{w}})$ and $\ell_{i}({\bm{w}})$ under this Gaussian approximation:

The local BIF is the negative of this covariance:

This expression is exactly the form of the classical dampened influence function, where the localization strength $\gamma$ serves as the dampening coefficient. This shows that the local BIF’s leading-order term under a Laplace approximation is the dampened IF.

Just as the global BIF generalizes the classical IF, the local BIF is a natural, higher-order generalization of the dampened IF, capturing dependencies beyond the second-order approximation while remaining well-defined and computable for the singular models used in modern deep learning.

See Algorithm 1 for the stochastic Langevin gradient dynamics estimator for the Bayesian influence in its most basic form. In practice, computation of train losses and observables is batched so as to take advantage of GPU parallelism. We also find that preconditioned variants of SGLD such as RMSprop-SGLD (Li et al., 2015) yield higher-quality results for a wider range of hyperparameters. We use an implementation provided by van Wingerden et al. (2024).

The SGLD update step described here, which is the one we use in our experiments, differs slightly from the presentation in the main text: we introduce a scalar inverse temperature $\beta$ (separate from the per-sample perturbations $\bm{\beta}$). Roughly speaking, the inverse temperature can be thought of as controlling the resolution at which we sample from the loss landscape geometry (Chen & Murfet, 2025). An alternative viewpoint is that the effective dataset size of training by iterative optimization is not obviously the same as the training dataset size $n$ used in the Bayesian setting; we scale by $\beta$ to account for this difference. Hence, in practice, we combine $\beta n$ as a single hyperparameter to be tuned.

Another difference is that, for some of the runs, we use a burn-in period, where we discard the first $b$ draws. Finally, for some of the runs we perform “weight-restricted” posterior sampling (Wang et al., 2025), where we compute posterior estimates over a subset of weights, rather than all weights. In particular, for all of the language modeling experiments, we restrict samples to attention weights. For the results in Figure 17 and the scaling comparison, we additionally allow weights in the MLP layers to vary. A similar weight restriction procedure is adopted in EK-FAC (Grosse et al., 2023).

Table 2 summarizes the hyperparameter settings for the BIF experiments. The hyperparameters refer to the Algorithm 1: $m$ is the batch size, $C$ is the number of chains, $T$ the number of draws per chain, $b$ is the number of burn-in steps, $\epsilon$ is the learning rate, $\beta$ is the inverse temperature, and $\gamma$ is the localization strength. See Appendix B.1 for more details on each of these hyperparameters.

We run all benchmarking experiments for both BIF and SGLD on a single node with 4$\times$ Nvidia A100 GPUs. As given in Table 2, for the BIF estimation, we run SGLD with batch size $m=32$, number of chains $C=4$, number of draws per chain $T=500$, learning rate $\epsilon=5\times 10^{-6}$, inverse temperature $n\beta=30$, and localization strength $\gamma=300$. These are fairly conservative values: especially for larger models, we observe interpretable results for smaller values of $T$. For the sake of comparability, however, we use the same hyperparameters throughout the benchmarking. Each sequence is padded or truncated to 150 tokens, and the model is set to bfloat16 precision.

We use the kronfluence package for EK-FAC computation (Grosse et al., 2023).¹¹1The corresponding github repository is available here: https://github.com/pomonam/kronfluence This package splits the influence computation into a fit and score step. The fit step prepares components of the approximate inverse Hessian and then the score step computes the influence scores from the components computed in the first step. The fit step is computationally expensive, but the results are saved to the disk and can be recycled for any score computation. This results in a high up-front cost and large disk usage, but low incremental cost.

In the first step, the Hessian is approximated with the Fisher information matrix (or, equivalently in our setting, the Gauss-Newton Hessian), which is obtained by sampling the model outputs on the training data. Since the Pile, which is the dataset used for Pythia training, is too large to iterate over in full, we approximate it by taking a representative subset of $1\,000\,000$ data points, curated using $k$-means clustering (Gao et al., 2021; Kaddour, 2023). Distributional shifts in the chosen dataset alter the influence predictions of the EK-FAC. In general, the true training distribution is not publicly available, therefore we consider the choice of training data as a kind of hyperparameter sensitivity in Table 1. Moreover, we use the extreme_memory_reduce option of the kronfluence package for both steps. Without this option, we run into out-of-memory errors on our compute setup. Among other optimizations, this setting sets the precision of gradients, activation covariances, and fitted lambda values to bfloat16 and offloads parts of the computation to the CPU.

The comparison is depicted in Figure 5. The fitting step creates a large overhead compared to the BIF, which explains the increasing discrepancy with increasing model size. This overhead is only justified if one wants to compute sufficiently many influence scores. Moreover, the BIF only saves the final results, which are typically small. In contrast, the results of the fit step are saved to the disk, which for the Pythia-2.8B model occupies 41 GiB.

Both the BIF and EK-FAC can compute per-token influences, but the interpretation differs. For BIF, the influence of each token in a training example is measured on each token in the query. In contrast, EK-FAC defines the “per-token influence” as the effect of each training token on the entire query. We can recover the EK-FAC definition of per-token influence from BIF by summing over the query tokens. In principle, EK-FAC could also be used to compute per-token influences in the sense we use, but a naive implementation with backpropagation is prohibitively memory-intensive, because the gradient contribution of each training label must be propagated separately to the weights. Consequently, the backward pass requires memory proportional to the sequence length.

Both classical and Bayesian influence functions approximate the effect of ${\mathbf{z}}_{i}$’s exclusion from $\mathcal{D}_{\text{train}}$ as linear. That is, given a subset ${\mathcal{D}}\subseteq\mathcal{D}_{\text{train}}$, write $\phi({\mathcal{D}})$ as the value of the observable $\phi$ corresponding to a model trained on ${\mathcal{D}}$:

in the classical perspective and

in the Bayesian perspective. In either case, we approximate $\phi(\mathcal{D})$ as linear in the set $\mathcal{D}$:

where each $\tau_{i}\in\mathbb{R}$ is a training data attribution measure associated to ${\mathbf{z}}_{i}$ and $\phi$, e.g. $\operatorname{IF}({\mathbf{z}}_{i},\phi)$ or $\operatorname{BIF}({\mathbf{z}}_{i},\phi)$.

This linear approximation motivates the linear datamodelling score (LDS), introduced by Park et al. (2023). Given the training dataset $\mathcal{D}_{\text{train}}$ of cardinality $n$ and a query set $\mathcal{D}_{\text{query}}$, we let the query losses $(\phi_{{\mathbf{z}}_{j}}=\ell({\mathbf{z}}_{j};-))_{{\mathbf{z}}_{j}\in\mathcal{D}_{\text{query}}}$ be our observables and suppose we are given TDA measures $(\bm{\tau}_{{\mathbf{z}}_{j}})_{{\mathbf{z}}_{j}\in\mathcal{D}_{\text{query}}}$, with each $\bm{\tau}_{{\mathbf{z}}_{j}}\in\mathbb{R}^{n}$. To measure the LDS of $(\bm{\tau}_{{\mathbf{z}}_{j}})_{{\mathbf{z}}_{j}}$, we subsample datasets $\{\mathcal{D}_{k}\}_{k=1}^{K}$ with each ${\mathbf{z}}_{i}\in\mathcal{D}_{k}$ with probability $\alpha_{\text{retrain}}\in\{0.1,0.3,0.5,0.7\}$ iid. (For our experiments, we set $K=100$). The LDS of $(\bm{\tau}_{{\mathbf{z}}_{j}})_{{\mathbf{z}}_{j}}$ is then the average over $1\leq k\leq K$ of the correlation between the true retrained observable and the linear approximation from $(\bm{\tau}_{{\mathbf{z}}_{j}})_{{\mathbf{z}}_{j}}$:

where $\rho_{\text{s}}$ is Spearman’s rank correlation coefficient. Each $\phi_{\text{C},{{\mathbf{z}}_{j}}}(\mathcal{D}_{k})$ is computed by retraining the model on $\mathcal{D}_{k}$ and evaluating the loss on ${{\mathbf{z}}_{j}}$. Note that, regardless of whether we evaluate the LDS of an approximate classical IF or the BIF, we use the classical version of the retrained observable $\phi_{\text{C}}$. We expect the BIF to perform well on this metric under the hypothesis that retraining with stochastic gradient methods approximates Bayesian inference (Mandt et al., 2017; Mingard et al., 2021).

We evaluate the LDS of the EK-FAC, BIF, GradSim, TRAK on a ResNet-9 model with $1\,972\,792$ parameters (He et al., 2015) trained on the CIFAR-10 (Krizhevsky, 2009) image classification dataset. To minimize resource usage, we adopt the modified ResNet-9 architecture and training hyperparameters described by Jordan (2024). In addition, we set aside a warmup set $\mathcal{D}_{\text{warmup}}$ of $2500$ images. Before the actual training runs, we perform a short warmup phase on $\mathcal{D}_{\text{warmup}}$ to prime the optimizer state. The training hyperparameters are summarized in Tab. 3.

As described in Appendix C.1, we evaluate LDS by re-training the ResNet-9 100 times from initialization on random subsamples of the full CIFAR-10 training set, excluding the warmup set ($n=47\,500$ images). Each subsample contains $n_{\text{retrain}}=\alpha_{retrain}\alpha_{attribution}n$ images. We then use the full test set ($q=10\,000$ images) as the query set, i.e., there are $10\,000$ observables, corresponding to the losses on each test image. Thus, both EK-FAC and BIF TDA scores comprise a $n_{\text{attribution}}\times$10\,000$$ matrix. The hyperparameters for the SGLD estimation of the BIF are given in Tab. 2. For EK-FAC, we set the dampening factor to $10^{-8}$. Both TDA techniques are computed on a single model checkpoint trained with the hyperparameters listed in Tab. 3. Figure 6 displays the wall-clock times of the BIF and EK-FAC computation. In these experiments, EK-FAC is around five times faster than the BIF. This advantage is largely due to the small model sizes ($\sim 2\times 10^{6}$ parameters), which results in a short fitting stage.

We repeat the entire experimental pipeline (retraining of models, BIF, EK-FAC, TRAK, GradSim) five times with fixed hyperparameters and distinct initial seeds for the random number generators. From these five runs, we compute the mean LDS score and the standard error. The LDS scores of each individual run are displayed in Figure 7. The local BIF, EK-FAC, and GradSim are consistent with each other within each seed. However, the LDS score varies substantially across seeds. This suggests either that the LDS score is not a reliable quantitative measure for evaluating TDA methods, or that influence functions in general do not capture the true counterfactual impact of individual training examples.

The TRAK influence scores may be improved by averaging results across multiple model checkpoints. Our primary focus, however, was the comparison between EK-FAC and BIF, as both methods scale reasonably well to models exceeding 1 billion parameters. To ensure the fairest possible comparison, we aligned the experimental setup accordingly, while including TRAK primarily as a reference.

Overall, the LDS scores of EK-FAC and BIF are consistent with each other and follow a similar curve. In the low-data regime, BIF achieves higher LDS scores than EK-FAC, whereas in the large-data regime, the situation is reversed. As we show in Section A.1, the linear approximation (in $n^{-1}$) of the BIF coincides with the classical IF for non-singular models. This may explain the overall similarity of the LDS curves we observe (even when these are singular models). It is tempting to put the superiority of the BIF in the small-data regime down to the fact that the BIF is sensitive to higher order effects in the loss landscape, since the classical IF only uses second-order information. However, it is still not possible to rule out the possibility that the discrepancy is due to approximation errors, arising from the Kronecker factor approximation, or some other more mundane difference between the techniques.

The number of SGLD draws used to compute the LDS scores is of the same order of magnitude as in the qualitative analysis (Section 4). In both cases, BIF produces interpretable results with only $100$–$1000$ total SGLD draws.

We analyzed the dependence on the SGLD hyperparameters by sweeping over $(b,n\beta,\gamma)\in[0,100]\times[100,300,1,000,3,000]\times[1,000,3,000,10,000,30,000,100,000]$, using $\alpha_{\text{attribution}}=0.1$ and computing the corresponding LDS scores. The grid plots Figure 9–Figure 12 show the resulting loss traces and LDS scores for $\alpha_{\text{retrain}}=0.1$ and $\alpha_{\text{retrain}}=0.3$. These comparisons indicate that for $b=100$, the LDS scores remain stable across hyperparameter choices as long as the loss trace converges. Furthermore, Figure 13 demonstrates that this stability holds independently of the choice of $\alpha_{\text{retrain}}$.

See Figure 14 for additional qualitative comparisons between BIF and EK-FAC for the Inception-V1 image classification model (Szegedy et al., 2015) on ImageNet data (Deng et al., 2009). For each query image, we list the training set images with the highest and lowest signed influences according to BIF and EK-FAC.

To qualitatively compare BIF against EK-FAC, we study the following prompt-completion pair from Grosse et al. (2023):

Human: What is your objective function? Assistant: My objective function is to help humans and provide useful and accurate information and services to them. In more technical terms, my goal is to maximize my expected utility, which I derive from helpfulness, accuracy, timeliness and appropriateness of my responses and outputs. Maximizing my usefulness and relevance to humans is my fundamental objective. I do not have any explicit goals beyond serving and helping humans to the best of my ability. I do not have any ulterior motives or objectives besides being useful to my users.

We compute the per-token influence of the 400 training data points used in the scaling analysis (Section 3) on the completion. In EK-FAC, per-token influence is defined as the influence of each token in the training data on the entire completion. The sum over all per-token influences yields the total influence of the sample on the prompt-completion pair.

Here we show additional examples for the per-token BIF on Pythia 2.8B (Figure 17) and Pythia 14M (Figures 18 and 19).