Data Attribution

Contents

• Introduction
• 1. Causality and Counterfactuals
  • 1.1 Data attribution as causal analysis
  • 1.2 Counterfactual attribution and its failures
  • 1.3 Beyond leave-one-out: Shapley values
  • 1.4 Notation and goals
• 2. Influence Functions
  • 2.1 The influence function formula
  • 2.2 Translation to modern neural networks
• 3. Bayesian Influence Functions
  • 3.1 Bayesian influence functions
  • 3.2 Exercise: Connecting Bayesian and classical influence functions
  • 3.3 Long Exercise: Influence functions as optimal linear transport
• 4. Unrolling
  • 4.1 The training-dynamics approach to attribution
  • 4.2 The unrolling formula
    • 4.2.1 Materializing the Jacobian
    • 4.2.2 Implicit JVPs and the REPLAY algorithm
  • 4.3 Long Exercise: From unrolling to influence functions
  • 4.4 Exercise: Influence depends on training time
• 5. Practical Considerations and Open Problems
• 6. Further Readings

Introduction

Contents. These notes are structured as follows:

1. Causality and Counterfactuals motivates the problem of data attribution through the lens of causal reasoning, exhibiting the shortcomings of naive counterfactual approaches.
2. Influence Functions introduces the classical theory of influence functions for statistical models and discusses their translation to modern neural networks.
3. Bayesian Influence Functions introduces susceptibilities and Bayesian influence functions, connecting them to the classical theory via a power series expansion.
4. Unrolling presents the training-dynamics approach to attribution, deriving the basic formula and connecting it to influence functions in the appropriate limit.
5. Practical Considerations and Open Problems discusses distributional versus single-model attribution, entangled latent concepts, and the distinction between similarity metrics and genuine attribution.

Prerequisites. We assume the following background, most of it standard for a machine-learning audience.

• Linear algebra: vectors, matrices, inner products, eigenvalues, quadratic forms, Hessians.
• Multivariable calculus: gradients, the chain rule, Taylor expansion, the implicit function theorem.
• Probability: expectations, conditional probability, Bayes' rule, Gaussian distributions.
• Deep learning basics: stochastic gradient descent, empirical risk minimisation, a standard supervised training loop.

Familiarity with the material from the singular learning theory tutorial (degeneracy, the local learning coefficient, posterior broadening) is useful but not required. We flag connections to the SLT day in remarks as they arise.

Fast-track. There is more material here than comfortably fits into a single day. Readers short on time may prefer one of the following routes; each is self-contained.

1. From classical to modern influence functions. Skim Section 1 for motivation, then read Section 2 in full, including the derivation exercise.
2. The Bayesian perspective. Read Section 2.1 for the influence function formula, then Section 3 in full. The connection exercise (Section 3.2) is the conceptual payoff.
3. Training-dynamics attribution. Skim Section 2, then read Section 4 in full. The long exercise in Section 4.3 recovers the influence function as a limit of unrolling.

All three routes depend on the motivational material in Section 1. If reading only one section, read that one.

Acknowledgements. This tutorial was developed for the April 2026 Iliad Intensive and its template was inspired by the SLT day. Exercises and examples adapted from specific sources are credited inline.

1. Causality and Counterfactuals

1.1 Data attribution as causal analysis

At the heart of data attribution, we find a causal question: Which training examples caused a model to behave the way it does? If mechanistic interpretability aims for the causal analysis of a forward pass, data attribution aims for the causal analysis of a training run.¹ Our north star is a decomposition of the training data into "causes" of the model's behaviour.

Two asymmetries between the two settings are worth keeping in mind from the start. First, both problems inherit the same conceptual difficulties of causal reasoning in general (explained in the next section) and there is no reason to expect data attribution to escape them. Second, a forward pass is cheap, while a training run is not. Where mechanistic interpretability can validate a claim by running many ablations, data attribution cannot in general afford even one retraining. That being said, there are also advantages: Generally speaking, the training process is more robust and the ablations we run may be more human understandable than high-dimensional modification of activations or weights.

Now, before developing such data attribution methods, we should be aware of the tension between causality and the ways we will validate a hypothesis. Finding a "true" causal relationship in practice can be much more difficult than it may seem at first glance. In the following section we take the question seriously, and find that the naive counterfactual answer fails in ways that matter.

1.2 Counterfactual attribution and its failures

The classical notion of causation is due to Lewis 1973: $A$ caused $B$ iff, had $A$ not occurred, $B$ would not have occurred. Applied to data attribution, this says that a training example "causes" a behaviour when removing it would have altered the behaviour, also known as the leave-one-out (LOO) counterfactual. We relax this definition from a binary to a gradual one.

Definition 1.1 (Counterfactual attribution). We attribute an action $A$ to a behaviour $B$ iff

1. had $A$ occurred, $B$ would have occurred to a higher degree;
2. had $A$ not occurred, $B$ would have occurred to a lesser degree (or not at all).

Definition 1.1 is the simplest formalisation of cause. But as Mueller 2024 emphasises, it breaks in at least two cases.

Overdetermination: redundant causes are missed. When two training examples are each sufficient to produce a behaviour, removing either one leaves the behaviour intact, and LOO attributes credit to neither.

The same pattern may arise in learning, where once a behaviour is learned, i.e. has hit saturation, additional examples will not affect it and appear as causally irrelevant.

Non-transitivity: counterfactuals do not chain. Direct ablation will not capture how an action causes a behaviour.

In the data-attribution setting this corresponds to chains of mediation. A training example may induce an intermediate feature, and that feature may drive a downstream capability, without the datapoint itself passing the counterfactual test against the capability end-to-end. Methods that look only at marginal effects miss the mediation.

1.3 Beyond leave-one-out: Shapley values

A natural response to overdetermination is to stop looking only at the LOO swap and instead ask how a training example contributes across all subsets of the training set it could be part of. This is the idea behind Shapley values (Shapley 1953).

Definition 1.4 (Shapley value). Let $v:2^{[n]}\to \mathbb{R}$ be a set function (a "game"). The Shapley value of element $i$ is

$${\phi }_i(v):=\frac{1}{n}\sum _{S\subseteq [n]\setminus \{i\}}{\left(\genfrac{}{}{0ex}{}{n-1}{|S|}\right)}^{-1}(v(S\cup \{i\})-v(S)).$$

We can think of $[n]$ as a set of players and $v(S)$ as the value of the coalition. In our previous example, the players are the premium ingredients and the value is the resulting fluffiness. Then $v(S)$ will measure the fluffiness of the croissant when the ingredients in $S$ are premium and the rest are standard.

Shapley values partially address overdetermination. Returning to the baker example: in coalitions $S\subset [n]$ where neither butter nor yeast is premium, upgrading one has a large marginal effect; in coalitions where the other leavener is already premium, the marginal effect is small. Averaging gives each ingredient roughly half the credit.

1.4 Notation and goals

We close the section by fixing notation that will be used throughout.

• A training set is $\mathcal{D}=\{z_i{\}}_{i=1}^n$ with $z_i=(x_i,y_i)$.
• A model is a function $f_w:X\to Y$ parameterised by $w\in W\subset {\mathbb{R}}^{d_p}$.
• The per-example loss is $L_j(w)=L(w,z_i)=\ell (f_w(x_i),y_i)$ for some pointwise loss $\ell$.
• We introduce a vector of data weights $\beta =({\beta }_1,\dots ,{\beta }_n)\in {\mathbb{R}}^n$ and define the weighted training loss

$$L(w;\beta ):=\frac{1}{n}\sum _{i=1}^n{\beta }_{i}L(w,z_i),$$

and we write $L(w):=L(w,1)$ for the unweighted loss.

• The behaviour we want to attribute is a differentiable scalar observable $\varphi :W\to \mathbb{R}$. This could be the loss on a held-out point, an average loss on a task, the logit of a specific completion, or any other scalar function of the trained parameters.

In general we want to know "If I change ${\beta }_i$, how will this affect $\varphi$?" This question is typically split into two parts:

• Parameter influence: "If I change ${\beta }_i$, how will this affect the final parameters?"
• Influence on $\varphi$: "If I change the final parameters, how will this affect $\varphi$?"

In practice we may think of the influence then as a measure of the composition of these maps

where we are intentionally vague about the first map and where the different approaches, influence functions, Bayesian influence functions, and unrolling, differ in how they interpret it.

That being said, all attribution methods in these lecture notes can be phrased as estimators of partial derivatives or finite differences of this map. The influence of an example $i$ is then computed as first-order approximation:

2. Influence Functions

In this section we develop influence functions: a closed-form, first-order approximation to the leave-one-out (LOO) counterfactual. The classical theory (Hampel 1974; Cook 1977) predates machine learning by decades. It was developed in robust statistics, where the question was which observations a parameter estimate is most sensitive to. For simple statistical models the answer is short: a single application of the implicit function theorem yields an inverse-Hessian-times-gradient formula whose factors have a clean operational interpretation.

The formula we arrive at then can be translated to the setup of neural networks, but its interpretation becomes unclear: The classical derivation assumes a unique global minimum and an invertible Hessian, neither of which holds for a neural network which is typically not trained to full convergence. We will see a partial fix for this issue in Section 2.2.

2.1 The influence function formula

Influence functions originate in robust statistics (Hampel 1974): given a statistical estimator, how much does it change when a single observation is added or perturbed?

We now translate this idea into the setting of Section 1.4, whose notation we use throughout. Recall that the weighted training loss $L(w;\beta )={\sum }_i{\beta }_{i}L(w,z_i)$ defines a family of optimization problems parameterized by the weight vector $\beta$. In the case of influence functions, we measure the influence of a data point on the optimal parameters

$$w^{\ast }(\beta ):=\underset{w\in W}{\mathrm{argmin}}L(w;\beta ),$$

and we write $w^{\ast }:=w^{\ast }(1)$ for the trained parameters. We assume throughout this subsection that $w^{\ast }$ is the unique global minimum of $L$ and that the Hessian $H:={\nabla }_w^{2}L(w^{\ast })$ is positive definite, so that the implicit function theorem makes $w^{\ast }(\beta )$ well-defined and smooth in a neighborhood of $\beta =1$. In the notation of Equation 1 we set $\pi (\beta )=w^{\ast }(\beta )$ and Equation 3 yields then the following formula:

Definition 2.2 (Influence function). The parameter influence of training example $z_i$ at the trained parameters $w^{\ast }$ is

$${\mathcal{I}}_{\mathrm{param}}(z_i):=\frac{\partial w^{\ast }(\beta )}{\partial {\beta }_i}{|}_{\beta =1}.$$

Under the assumptions above, the implicit function theorem gives (Exercise 2.1)

$$\boxed{{\mathcal{I}}_{\mathrm{param}}(z_i)=-H^{-1}{\nabla }_{w}L(w^{\ast },z_i).}$$

where $H$ is the Hessian of $L$ at $w^{\ast }$. For any differentiable observable $\varphi :W\to \mathbb{R}$, the chain rule then yields the influence on $\varphi$:

$$\begin{array}{rl}\mathcal{I}(z_i,\varphi )& :=\frac{\partial \varphi (w^{\ast }(\beta ))}{\partial {\beta }_i}{|}_{\beta =1}\\ & ={\nabla }_w\varphi (w^{\ast }{)}^{\top }{\mathcal{I}}_{\mathrm{param}}(z_i)\\ & =-{\nabla }_w\varphi (w^{\ast }{)}^{\top }{H}^{-1}{\nabla }_{w}L(w^{\ast },z_i).\end{array}$$

Thus we get three components

• ${\nabla }_{w}L(w^{\ast },z_i)$ is the direction in parameter space along which $z_i$ pulls the optimizer.
• $H^{-1}$ rescales this direction by local curvature: directions in which the loss is flat (small eigenvalue of $H$) get amplified, since a small force in a flat direction moves the minimum a lot.
• ${\nabla }_w\varphi (w^{\ast })$ projects the parameter change onto the direction that matters for the observable $\varphi$.

In this section, we will focus on ${\mathcal{I}}_{\mathrm{param}}(z_i)$. The influence on any observable $\varphi$ is obtained by computing the dot product ${\mathcal{I}}_{\mathrm{param}}(z_i)\cdot {\nabla }_w\varphi (w^{\ast })$.

Connection to leave-one-out. Setting $\beta =1-e_i$ (i.e., removing $z_i$) and applying the first-order Taylor expansion gives

i.e. the parameter influence is, to first order, the change in $w^{\ast }$ that would result from removing $z_i$ from the training set entirely.

Exercise 2.1 (Derivation of the influence function). Let $w^{\ast }(\beta )={\mathrm{argmin}}_{w\in W}L(w;\beta )$ as above and assume the Hessian $H={\nabla }_w^{2}L(w^{\ast })$ is invertible. Recall that the implicit function theorem states: if $F:W\times B\to {\mathbb{R}}^{\dim W}$ is $C^1$, $F(w^{\ast },1)=0$, and the Jacobian $\partial F/\partial w$ at $(w^{\ast },1)$ is invertible, then there is a unique $C^1$ map $\beta \mapsto w^{\ast }(\beta )$ on a neighborhood of $1$ with $F(w^{\ast }(\beta ),\beta )\equiv 0$, whose derivative is obtained by implicit differentiation of $F$.

(a) Apply the implicit function theorem to $F(w,\beta ):={\nabla }_{w}L(w;\beta )$ to conclude that $\beta \mapsto w^{\ast }(\beta )$ is well-defined and $C^1$ near $1$. Which hypothesis of the theorem uses the invertibility of $H$? Where do we use that $w^{\ast }$ is the unique global minimum (rather than merely a critical point)?

(b) By differentiating the first-order optimality condition ${\nabla }_{w}L(w^{\ast }(\beta );\beta )=0$ with respect to ${\beta }_i$, show that

$$\frac{\partial w^{\ast }(\beta )}{\partial {\beta }_i}{|}_{\beta =1}=-H^{-1}{\nabla }_{w}L(w^{\ast },z_i).$$

(c) Use the chain rule to obtain the formula for $\mathcal{I}(z_i,\varphi )$ in Definition 2.2. Make sure to keep track of the evaluation at $\beta =1$.

Solution

(a) Set $F(w,\beta ):={\nabla }_{w}L(w;\beta )={\sum }_j{\beta }_j{\nabla }_{w}L(w,z_j)$. Then $F$ is $C^1$ in $(w,\beta )$ (assuming $L(\cdot ,z_j)$ is $C^2$), and $F(w^{\ast },1)={\nabla }_{w}L(w^{\ast })=0$ because $w^{\ast }$ is a critical point of the unperturbed loss. The Jacobian $\partial F/\partial w$ at $(w^{\ast },1)$ equals ${\nabla }_w^{2}L(w^{\ast })=H$, and the invertibility hypothesis of the implicit function theorem is exactly the assumption that $H$ is invertible. The theorem then yields a unique $C^1$ map $\beta \mapsto w^{\ast }(\beta )$ on a neighborhood of $1$ with $F(w^{\ast }(\beta ),\beta )\equiv 0$.

The uniqueness the implicit function theorem provides is only local: $w^{\ast }(\beta )$ is uniquely defined as the critical point near $w^{\ast }$. The assumption that $w^{\ast }$ is the unique global minimum is what lets us identify this local branch of critical points with the object we actually care about — the argmin. Without it, ${\mathrm{argmin}}_{w}L(w;\beta )$ could jump to a different basin as $\beta$ varies, and the response function would fail to be continuous (let alone differentiable).

(b) Define the two-slot function

$$G(w,\beta ):={\nabla }_{w}L(w;\beta )=\sum _j{\beta }_j{\nabla }_{w}L(w,z_j),$$

so that the optimality condition reads $G(w^{\ast }(\beta ),\beta )\equiv 0$. The variable $\beta$ enters this identity in two places: once through the first slot via $w^{\ast }(\beta )$, and once directly through the second slot. The multivariate chain rule therefore produces two terms when we differentiate with respect to ${\beta }_i$:

$$\underset{\text{indirect: }{w}^{\ast }\text{ moves with }\beta }{\underbrace{\frac{\partial G}{\partial w}{|}_{(w^{\ast }(\beta ),\beta )}\frac{\partial w^{\ast }(\beta )}{\partial {\beta }_i}}}+\underset{\text{direct: explicit }{\beta }_i\text{ dependence}}{\underbrace{\frac{\partial G}{\partial {\beta }_i}{|}_{(w^{\ast }(\beta ),\beta )}}}=0.$$

The two partials are:

• $\partial G/\partial w={\nabla }_w^{2}L(w;\beta )$, which at $\beta =1$ equals the Hessian $H$.
• $\partial G/\partial {\beta }_i$: differentiating ${\sum }_j{\beta }_j{\nabla }_{w}L(w,z_j)$ with respect to ${\beta }_i$ with $w$ held fixed picks out the single $j=i$ summand, giving ${\nabla }_{w}L(w,z_i)$.

Evaluating at $\beta =1$ and solving,

$$\frac{\partial w^{\ast }(\beta )}{\partial {\beta }_i}{|}_{\beta =1}=-H^{-1}{\nabla }_{w}L(w^{\ast },z_i).$$

(c) For a differentiable observable $\varphi :W\to \mathbb{R}$, the chain rule applied to $\beta \mapsto \varphi (w^{\ast }(\beta ))$ gives

$$\mathcal{I}(z_i,\varphi )=\frac{\partial \varphi (w^{\ast }(\beta ))}{\partial {\beta }_i}{|}_{\beta =1}={\nabla }_w\varphi (w^{\ast }(\beta ){)}^{\top }\frac{\partial w^{\ast }(\beta )}{\partial {\beta }_i}{|}_{\beta =1}=-{\nabla }_w\varphi (w^{\ast }{)}^{\top }{H}^{-1}{\nabla }_{w}L(w^{\ast },z_i),$$

where in the last step both factors are evaluated at $\beta =1$, i.e. at $w^{\ast }$.

Exercise 2.2 (Linear regression). Consider linear regression with training data $\{(x_i,y_i){\}}_{i=1}^n$, $x_i\in {\mathbb{R}}^p$, and squared loss $L(w,(x_i,y_i))=\frac{1}{2}(y_i-w^{\top }{x}_i{)}^2$. Write $A:={\sum }_i{x}_i{x}_i^{\top }$. In this case the minimizer, the influence function, and the exact LOO update all have closed-form solutions. We will see that in general the influence functions can diverge in extreme cases.

(a) Derive the minimizer $w^{\ast }$ by setting ${\nabla }_{w}L(w)=0$. Show that $w^{\ast }=A^{-1}{\sum }_i{x}_i{y}_i$.

(b) Compute the parameter influence ${\mathcal{I}}_{\mathrm{param}}(z_j)$ from Definition 2.2.

(c) Define $w^{\ast }({\beta }_j):={\mathrm{argmin}}_w{\sum }_{i\ne j}L(w,z_i)+{\beta }_{j}L(w,z_j)$, so that ${\beta }_j=1$ is the unperturbed model and ${\beta }_j=0$ corresponds to fully removing $z_j$. Using the Sherman–Morrison formula², show that

$$w^{\ast }({\beta }_j)-w^{\ast }=-\frac{(1-{\beta }_j)A^{-1}{x}_j{r}_j}{1-(1-{\beta }_j)h_j},$$

where $r_j:=y_j-x_j^{\top }{w}^{\ast }$ is the residual at $z_j$ and $h_j:=x_j^{\top }{A}^{-1}{x}_j$ is the leverage of $x_j$. Note that this is a rational function of ${\beta }_j$. This is the ground truth effect of perturbing $z_j$.

(d) Expand the closed form from (c) as a Taylor series in $(1-{\beta }_j)$ around the unperturbed point ${\beta }_j=1$ by recognizing $1/(1-(1-{\beta }_j)h_j)$ as a geometric series, and obtain

$$w^{\ast }({\beta }_j)-w^{\ast }=-A^{-1}{x}_j{r}_j\sum _{k=1}^{\infty }(1-{\beta }_j{)}^k{h}_j^{k-1}.$$

Identify the $k=1$ term with the parameter influence from part (b), and observe that the higher-order terms ($k\ge 2$) are suppressed by powers of the leverage $h_j$.

(e) Specialize to the LOO endpoint ${\beta }_j=0$. Sum the geometric series in $h_j$ to obtain the closed-form deletion update

$$w_{(-j)}^{\ast }-w^{\ast }=-A^{-1}{x}_j{r}_j(1+h_j+h_j^2+h_j^3+\cdots )=-\frac{A^{-1}{x}_j{r}_j}{1-h_j}.$$

The IF prediction is the $h_j^0$ leading term. Discuss the two regimes $h_j\to 0$ and $h_j\to 1$: when does keeping only the linear term suffice, and when do all the higher-order corrections matter?

Solution

(a) The gradient of the training loss is ${\nabla }_{w}L(w)=-{\sum }_i{x}_i(y_i-x_i^{\top }w)=-{\sum }_i{x}_i{y}_i+Aw$. Setting this to zero gives $A{w}^{\ast }={\sum }_i{x}_i{y}_i$, hence $w^{\ast }=A^{-1}{\sum }_i{x}_i{y}_i$.

(b) The Hessian of the unweighted training loss is $H={\sum }_i{x}_i{x}_i^{\top }=A$, and the gradient of the per-example loss at $z_j=(x_j,y_j)$ is ${\nabla }_{w}L(w^{\ast },z_j)=-x_j(y_j-x_j^{\top }{w}^{\ast })=-x_j{r}_j$. Plugging into Definition 2.2 with $\varphi (w)=w$ (so ${\nabla }_w\varphi (w^{\ast })=I$),

$$\mathcal{I}(z_j,w^{\ast })=-A^{-1}(-x_j{r}_j)=A^{-1}{x}_j{r}_j,$$

a residual times a leverage-weighted input.

(c) The normal equations for $w^{\ast }({\beta }_j)$ read

$$(\sum _{i\ne j}{x}_i{x}_i^{\top }+{\beta }_j{x}_j{x}_j^{\top })w=\sum _{i\ne j}{x}_i{y}_i+{\beta }_j{x}_j{y}_j,$$

or equivalently $(A-(1-{\beta }_j)x_j{x}_j^{\top })w={\sum }_i{x}_i{y}_i-(1-{\beta }_j)x_j{y}_j$. Sherman–Morrison applied to the rank-one update on the left gives

$$(A-(1-{\beta }_j)x_j{x}_j^{\top }{)}^{-1}=A^{-1}+\frac{(1-{\beta }_j)p_j{p}_j^{\top }}{1-(1-{\beta }_j)h_j},$$

where we write $p_j:=A^{-1}{x}_j$ and $h_j:=x_j^{\top }{p}_j$. Multiplying through and using the identities $A^{-1}{\sum }_i{x}_i{y}_i=w^{\ast }$ and $p_j^{\top }{\sum }_i{x}_i{y}_i=x_j^{\top }{w}^{\ast }$,

$$\begin{array}{rl}{w}^{\ast }({\beta }_j)& =w^{\ast }-(1-{\beta }_j)p_j{y}_j+\frac{(1-{\beta }_j)p_j(x_j^{\top }{w}^{\ast }-(1-{\beta }_j)h_j{y}_j)}{1-(1-{\beta }_j)h_j}\\ & =w^{\ast }+p_j\cdot \frac{-(1-{\beta }_j)y_j(1-(1-{\beta }_j)h_j)+(1-{\beta }_j)x_j^{\top }{w}^{\ast }-(1-{\beta }_j{)}^2{h}_j{y}_j}{1-(1-{\beta }_j)h_j}\\ & =w^{\ast }+p_j\cdot \frac{(1-{\beta }_j)(x_j^{\top }{w}^{\ast }-y_j)}{1-(1-{\beta }_j)h_j}\\ & =w^{\ast }-\frac{(1-{\beta }_j)p_j{r}_j}{1-(1-{\beta }_j)h_j},\end{array}$$

the claimed expression. The result is rational in ${\beta }_j$ because the matrix being inverted is linear in ${\beta }_j$.

(d) Expanding $1/(1-(1-{\beta }_j)h_j)$ as a geometric series in $(1-{\beta }_j)h_j$ (valid for $|(1-{\beta }_j)h_j|<1$),

$$w^{\ast }({\beta }_j)-w^{\ast }=-(1-{\beta }_j)A^{-1}{x}_j{r}_j\sum _{k=0}^{\infty }(1-{\beta }_j{)}^k{h}_j^k=-A^{-1}{x}_j{r}_j\sum _{k=1}^{\infty }(1-{\beta }_j{)}^k{h}_j^{k-1}.$$

Reading off term by term:

• The $k=1$ term is $-A^{-1}{x}_j{r}_j(1-{\beta }_j)=\mathcal{I}(z_j,w^{\ast })\cdot ({\beta }_j-1)$, exactly the influence-function prediction — as it must be, since the IF is the linear coefficient of $w^{\ast }({\beta }_j)$ at ${\beta }_j=1$.
• The $k=2$ term, $-A^{-1}{x}_j{r}_j(1-{\beta }_j{)}^2{h}_j$, is the second-order Taylor correction. It carries one factor of leverage and is doubly suppressed: small both in $(1-{\beta }_j)$ and in $h_j$.
• Each subsequent term picks up another factor of $(1-{\beta }_j)h_j$. The IF/LOO "error" is precisely this tail, viewed as a series in the perturbation size and the leverage.

(e) At ${\beta }_j=0$ the perturbation factor $(1-{\beta }_j)$ is unity, and the higher-order terms no longer have a small parameter in front of them — only powers of $h_j$ remain. The series collapses to a pure geometric series in $h_j$:

$$w_{(-j)}^{\ast }-w^{\ast }=-A^{-1}{x}_j{r}_j(1+h_j+h_j^2+h_j^3+\cdots )=-\frac{A^{-1}{x}_j{r}_j}{1-h_j}.$$

The IF prediction $-A^{-1}{x}_j{r}_j$ is the $k=0$ leading term; everything else is the Taylor remainder. So the often-quoted "factor of $1/(1-h_j)$" between IF and LOO is not a non-perturbative effect: it is the sum of all the higher-order Taylor terms, which OLS happens to admit in closed form because $w^{\ast }({\beta }_j)$ is rational. For more general M-estimators the response function is no longer rational and the higher-order terms do not sum to such a clean expression — but they are still there, with the same qualitative behavior.

In the two limiting regimes:

• $h_j\to 0$ (low leverage). Each higher-order term carries a factor of $h_j^{k-1}$, so they are individually negligible and the IF approximation is essentially exact.
• $h_j\to 1$ (high leverage). The geometric series barely converges and the higher-order terms together carry an arbitrarily large total weight. The IF underestimates the true LOO change without bound.

This is a clean demonstration that even in a strictly convex, fully-converged setting the linearization in Definition 2.2 loses control over precisely the points one most wants to attribute — the ones the model could not fit without them. The lesson generalizes: high-leverage / outlier examples are exactly where IF approximations should be distrusted, and the modern fixes from Section 2.2 (damping, GNH, PBRF) all address related but more severe versions of the same issue.

2.2 Translation to modern neural networks

The derivation in Section 2.1 relied on two assumptions:

(i) $w^{\ast }$ is the unique global minimum of the training loss;

(ii) the Hessian $H={\nabla }_w^{2}L(w^{\ast })$ is positive definite.

While neither holds for a modern neural network, the quantities in Definition 2.2 are, in principle, computable. Koh & Liang 2020 popularized the influence function formula as a tool for deep learning, and a substantial body of work has tried to make it useful at scale.

(1) Non-uniqueness and a singular Hessian. For a neural network the global minimum is usually not unique. The set of minimizers is generically a positive-dimensional zero locus, and the Hessian at any point on it is singular: it has a nontrivial null space corresponding to directions of parameter movement that leave the loss unchanged. This is the degeneracy phenomenon observed in the SLT day, and it means we cannot invert $H$ as in Definition 2.2.

The standard remedy is to dampen the Hessian: replace $H^{-1}$ with $(H+\lambda I{)}^{-1}$ for some $\lambda >0$. This is the same as adding an $L_2$-regularizer $\frac{\lambda }{2}\Vert w-w^{\ast }{\Vert }^2$ to the training objective $L(w,\beta )$; this yields a new loss function $L_{\mathrm{reg}}(w,\beta )=L(w,\beta )+\frac{\lambda }{2}\cdot \Vert w-w^{\ast }{\Vert }^2$.

(2) Negative curvature. Even after damping away the zero eigenvalues, $H$ at a non-converged checkpoint can have negative eigenvalues, in which case $(H+\lambda I{)}^{-1}$ may still be ill-conditioned or even point in the wrong direction. The standard fix is to replace $H$ with the Gauss–Newton Hessian (GNH)

$$G={\mathbb{E}}_{(x,y)\sim \mathcal{D}}\left[J^{\top }{H}_{y}J\right],$$

where $J=\partial f_w/\partial w$ is the parameter-output Jacobian and $H_y$ is the Hessian of the per-example loss with respect to the network outputs (not the parameters). For the cross-entropy and squared losses $H_y$ is positive semidefinite, hence so is $G$. The GNH is a local approximation to the true Hessian that ignores second derivatives of the network outputs with respect to the parameters. It is not hard to show that $G$ and $H$ differ by a term that vanishes at a critical point of $L$, so the GNH is a reasonable proxy for $H$ when $w^{\ast }$ is close to optimal.

(3) Non-converged checkpoints. In practice $w^{\ast }$ is not a critical point of $L$. The classical formula in Definition 2.2 is then evaluated at a point where $\nabla L\ne 0$.

The fix proposed by Bae et al. 2022 is to replace the original training loss by a surrogate for which $w^{\ast }$ is optimal by construction. Concretely, define the proximal Bregman objective (PBO) as

where $D_{L_y}(\hat{y},{\hat{y}}^{\ast })$ is the Bregman divergence of the output-space loss $L_y$ around ${\hat{y}}^{\ast }$. The first term penalizes any change in the network's predictions on the training set, with the current model $w^{\ast }$ as the reference; the second is the perturbation of interest; the third is the damping. By construction, $w^{\ast }$ is the minimizer of $L_{\mathrm{PBO}}(\cdot ;1)$ regardless of whether it is a critical point of the original loss. The minimizer of $L_{\mathrm{PBO}}(\cdot ;\beta )$ for nearby $\beta$ defines the proximal Bregman response function (PBRF).

When one re-derives the influence function formula starting from $L_{\mathrm{PBO}}$ instead of $L$, the result is

the same expression that the modern literature uses. This is effectively the counterfactual it approximates. As Bae et al. 2022 put it: if influence functions are the answer, the question is the PBRF, not LOO retraining. We refer to their paper for a detailed decomposition of the remaining gap between the PBRF and the original LOO counterfactual.

Exercise 2.3 (Failure of the first-order approximation in a degenerate model). Consider the toy two-parameter model from the SLT lecture day in which $w=(a,b)$ enters the loss only through the product $u=ab$, with squared loss $L(w,z)=\frac{1}{2}(z-ab{)}^2$ on training data $\{z_1,\dots ,z_n\}\subset \mathbb{R}$.

(a) Write down the Hessian of the training loss at any minimizer $(a^{\ast },b^{\ast })$ with $a^{\ast }{b}^{\ast }=\bar{z}:=\frac{1}{n}{\sum }_i{z}_i$. Show that it has rank one, and identify the one-dimensional null space geometrically.

(b) Compute the influence function $\mathcal{I}(z_i,\varphi )$ for the observable $\varphi (w)=ab$, treating $H^{-1}$ as a Moore–Penrose pseudoinverse³. How does it compare with the exact LOO change in $\varphi$?

(c) Repeat (b) using the damped inverse $(H+\lambda I{)}^{-1}$ for a small $\lambda >0$. Discuss how the answer depends on $\lambda$, on the choice of minimizer $(a^{\ast },b^{\ast })$ along the singular set, and on $n$. What does this say about the proximity gap?

Solution

(a) The per-example loss is ${\ell }_i=\frac{1}{2}(z_i-ab{)}^2$. Its second derivatives are

$$\frac{{\partial }^2{\ell }_i}{\partial a^2}=b^2,\frac{{\partial }^2{\ell }_i}{\partial b^2}=a^2,\frac{{\partial }^2{\ell }_i}{\partial a\partial b}=2ab-z_i.$$

Summing over $i$ and evaluating at any minimizer $(a^{\ast },b^{\ast })$ with $a^{\ast }{b}^{\ast }=\bar{z}$,

$$H=n\left(\begin{array}{cc}(b^{\ast }{)}^2& a^{\ast }{b}^{\ast }\\ a^{\ast }{b}^{\ast }& (a^{\ast }{)}^2\end{array}\right)=n\mathbf{v}{\mathbf{v}}^{\top },\mathbf{v}:=\left(\begin{array}{c}{b}^{\ast }\\ a^{\ast }\end{array}\right),$$

where we used ${\sum }_i(2a^{\ast }{b}^{\ast }-z_i)=n(2\bar{z}-\bar{z})=n{a}^{\ast }{b}^{\ast }$ for the off-diagonal. This is rank one, with image $\mathrm{span}(\mathbf{v})$.

The null space is $\mathrm{span}((a^{\ast },-b^{\ast }{)}^{\top })$. Geometrically: the minimizers form the hyperbola $\{(a,b):ab=\bar{z}\}$, and the tangent to this hyperbola at $(a^{\ast },b^{\ast })$ is exactly the direction $(a^{\ast },-b^{\ast })$ — the infinitesimal generator of the rescaling symmetry $(a,b)\to (\alpha a,b/\alpha )$ at $\alpha =1$. The null space of $H$ is the tangent to the manifold of minimizers.

(b) We need two gradients at the minimizer. For the observable $\varphi (w)=ab$:

$${\nabla }_w\varphi =(b^{\ast },a^{\ast }{)}^{\top }=\mathbf{v}.$$

For the per-example loss at $z_i$:

$${\nabla }_{w}L(w^{\ast },z_i)=-(z_i-\bar{z})(b^{\ast },a^{\ast }{)}^{\top }=-(z_i-\bar{z})\mathbf{v}.$$

Both point in the $\mathbf{v}$ direction — neither has a component in the null space of $H$. The Moore–Penrose pseudoinverse of $H=n\mathbf{v}{\mathbf{v}}^{\top }$ is

$$H^{+}=\frac{1}{n\Vert \mathbf{v}{\Vert }^4}\mathbf{v}{\mathbf{v}}^{\top },\Vert \mathbf{v}{\Vert }^2=(a^{\ast }{)}^2+(b^{\ast }{)}^2.$$

Plugging into the IF formula,

$$\mathcal{I}(z_i,\varphi )=-{\mathbf{v}}^{\top }{H}^{+}(-(z_i-\bar{z})\mathbf{v})=(z_i-\bar{z})\frac{\Vert \mathbf{v}{\Vert }^4}{n\Vert \mathbf{v}{\Vert }^4}=\frac{z_i-\bar{z}}{n}.$$

The $\Vert \mathbf{v}\Vert$ factors cancel: the pseudoinverse IF is independent of the choice of minimizer $(a^{\ast },b^{\ast })$.

Comparison with exact LOO. Removing $z_i$ changes the sample mean to ${\bar{z}}_{(-i)}=(n\bar{z}-z_i)/(n-1)$, and the observable at any new minimizer is $\varphi ={\bar{z}}_{(-i)}$. So the exact LOO change is

$$\varphi (w_{(-i)}^{\ast })-\varphi (w^{\ast })={\bar{z}}_{(-i)}-\bar{z}=-\frac{z_i-\bar{z}}{n-1}.$$

The IF predicts $-\mathcal{I}(z_i,\varphi )=-(z_i-\bar{z})/n$. The ratio is $n/(n-1)=1/(1-h)$ with leverage $h=1/n$ — every point has equal leverage in this one-effective-parameter model, and we recover the same $1/(1-h)$ correction as in Exercise 2.2.

(c) The eigenvalues of $H=n\mathbf{v}{\mathbf{v}}^{\top }$ are $n\Vert \mathbf{v}{\Vert }^2$ (eigenvector $\hat{\mathbf{v}}=\mathbf{v}/\Vert \mathbf{v}\Vert$) and $0$ (eigenvector ${\hat{\mathbf{v}}}_{\perp }=(a^{\ast },-b^{\ast }{)}^{\top }/\Vert \mathbf{v}\Vert$). So

$$(H+\lambda I{)}^{-1}=\frac{1}{n\Vert \mathbf{v}{\Vert }^2+\lambda }\hat{\mathbf{v}}{\hat{\mathbf{v}}}^{\top }+\frac{1}{\lambda }{\hat{\mathbf{v}}}_{\perp }{\hat{\mathbf{v}}}_{\perp }^{\top }.$$

Since $\nabla \varphi$ and $\nabla L_i$ both lie in the $\hat{\mathbf{v}}$ direction, the $1/\lambda$ term does not contribute, and

$${\mathcal{I}}_{\lambda }(z_i,\varphi )=\frac{(z_i-\bar{z})\Vert \mathbf{v}{\Vert }^2}{n\Vert \mathbf{v}{\Vert }^2+\lambda }.$$

Dependence on the minimizer. For $\lambda >0$ the answer depends on $\Vert \mathbf{v}{\Vert }^2=(a^{\ast }{)}^2+(b^{\ast }{)}^2$, which varies along the orbit $ab=\bar{z}$. For instance, the "balanced" minimizer $(\sqrt{\bar{z}},\sqrt{\bar{z}})$ gives $\Vert \mathbf{v}{\Vert }^2=2\bar{z}$, while the "unbalanced" choice $(\bar{z},1)$ gives $\Vert \mathbf{v}{\Vert }^2={\bar{z}}^2+1$. Different minimizers yield different IFs. This is the proximity gap: the damping term $\frac{\lambda }{2}\Vert w-w^{\ast }{\Vert }^2$ breaks the rescaling symmetry and anchors the computation at a particular point on the orbit.

Limit $\lambda \to 0$. ${\mathcal{I}}_{\lambda }\to (z_i-\bar{z})/n$ regardless of $(a^{\ast },b^{\ast })$, recovering the pseudoinverse answer from (b). The damping-induced dependence on the minimizer vanishes.

Limit $\lambda \to \infty$. ${\mathcal{I}}_{\lambda }\to 0$: heavy damping kills all attribution, because the proximal penalty makes it too costly to move $w$ at all.

The punchline: for this particular observable ($\varphi =ab$, which is constant along the orbit), the pseudoinverse gives a clean, minimizer-independent answer and the damping only introduces a spurious dependence. But for an observable that does vary along the orbit — say $\varphi (w)=a^2$ — $\nabla \varphi$ would have a component in the null direction $(a^{\ast },-b^{\ast })$, the $1/\lambda$ term would contribute, and the IF would diverge as $\lambda \to 0$. In that case damping is not a nuisance but a necessity, and the choice of minimizer genuinely matters. This is the core tension of the proximity gap: it is a distortion for orbit-invariant queries, but a regularizer for orbit-dependent ones.

3. Bayesian Influence Functions

Classical influence functions give a first-order approximation to the change in a point estimate $w^{\ast }$ when the training data is perturbed. The Bayesian reformulation replaces this point estimate by a distribution over parameters: if we perturb the training data, how does our posterior over the parameter space change? This leads naturally to the Bayesian generalization of influence functions introduced in Kreer et al. 2026. In the setting of Equation 1 we set $\pi (\beta )$ not to be a single weight, but a whole distribution over parameters $p_{\beta }$, namely the posterior.

Susceptibilities. The Bayesian influence function is one instance of a more general object studied by Baker et al. 2026. Given a tempered Gibbs posterior $p_{\beta }\propto e^{-\beta L(w)}\phi (w)$ at inverse temperature $\beta$, any pair of observables $\varphi ,\psi :W\to \mathbb{R}$ defines a susceptibility

The differentiation argument we will apply in Exercise 3.1 shows $\chi (\varphi |\psi )=\beta {\mathrm{Cov}}_{p_{\beta }}(\varphi ,\psi )$ in full generality. Taking $\psi =-L(\cdot ,z_i)$ gives the Bayesian influence function developed below. Other choices of $(\varphi ,\psi )$ yield other sensitivity notions, see e.g. the refined learning coefficients (Baker et al. 2026, Appendix D.2). We will not need the general framework again; we mention it only to locate the BIF within the wider susceptibility family.

3.1 Bayesian influence functions

We now specialize the susceptibility story to data attribution. Instead of asking "how does $\varphi (w^{\ast })$ change when we perturb the training weights?" (Definition 2.2) we ask "how does ${\mathbb{E}}_{w\sim p}[\varphi (w)]$ change?" The answer turns out to be remarkably clean: the derivative becomes a covariance, and the Hessian inversion of the classical formula is replaced by posterior sampling.

The tempered posterior. Given training data $\mathcal{D}=\{z_1,\dots ,z_n\}$ with per-sample losses $L(w,z_i)$, sample weights $\beta =({\beta }_1,\dots ,{\beta }_n)$, and a prior $\phi (w)$, define the tempered posterior

When $\beta =1$ and the loss is a negative log-likelihood, this is the standard Bayesian posterior. When the losses are not log-likelihoods, $p_{\beta }$ is a Gibbs measure; the mathematics is identical.

Definition 3.1 (Bayesian influence function). The Bayesian influence function (BIF) of training example $z_i$ on the observable $\varphi$ is the derivative of the posterior expectation of $\varphi$ with respect to the sample weight ${\beta }_i$:

$$\mathrm{BIF}(z_i,\varphi ):=\frac{\partial }{\partial {\beta }_i}{\mathbb{E}}_{w\sim p_{\beta }}[\varphi (w)]{|}_{\beta =1}.$$

The definition is the natural Bayesian analogue of Definition 2.2: replace $\varphi (w^{\ast }(\beta ))$ with ${\mathbb{E}}_{p_{\beta }}[\varphi (w)]$. The following exercise shows that this derivative has a clean closed form. Recall that for two real-valued random variables $X$ and $Y$ defined on the same probability space, the covariance is

$$\mathrm{Cov}(X,Y)=\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y].$$

When $X$ and $Y$ are functions of the random parameter $w\sim p$, we write ${\mathrm{Cov}}_{w\sim p}(X(w),Y(w))$ to indicate that the expectation is taken over the posterior $p$. If $\varphi$ is vector-valued ($\varphi :W\to {\mathbb{R}}^d$), the covariance $\mathrm{Cov}(X,\varphi )$ is the vector whose $j$-th entry is $\mathrm{Cov}(X,{\varphi }_j)$.

Exercise 3.1 (The BIF is a covariance). Let $p_{\beta }(w\mid \mathcal{D})$ be the tempered posterior of Equation 8 and write $p=p_1$ for the unperturbed posterior. By differentiating ${\mathbb{E}}_{p_{\beta }}[\varphi (w)]=\int \varphi (w)p_{\beta }(w\mid \mathcal{D})dw$ with respect to ${\beta }_i$, show that

$$\boxed{\mathrm{BIF}(z_i,\varphi )=-{\mathrm{Cov}}_{w\sim p}(L(w,z_i),\varphi (w)).}$$

Hint: Write $p_{\beta }\propto e^{-{\sum }_j{\beta }_{j}L(w,z_j)}\phi (w)$ and differentiate the ratio $\int \varphi p_{\beta }dw/\int p_{\beta }dw$ using the quotient rule. The key identity is ${\partial }_{{\beta }_i}\log Z(\beta )=-{\mathbb{E}}_{p_{\beta }}[L(w,z_i)]$.

Solution

Write $Z(\beta )=\int \exp (-{\sum }_j{\beta }_{j}L(w,z_j))\phi (w)dw$ for the normalizing constant, so that

$${\mathbb{E}}_{p_{\beta }}[\varphi (w)]=\frac{1}{Z(\beta )}\int \varphi (w)\exp (-\sum _j{\beta }_{j}L(w,z_j))\phi (w)dw.$$

Differentiating with respect to ${\beta }_i$ and applying the quotient rule:

$$\begin{array}{rl}\frac{\partial }{\partial {\beta }_i}{\mathbb{E}}_{p_{\beta }}[\varphi ]& =\frac{1}{Z}\int \varphi (w)(-L(w,z_i))e^{-{\sum }_j{\beta }_j{L}_j}\phi dw-\frac{1}{Z^2}\cdot \frac{\partial Z}{\partial {\beta }_i}\int \varphi e^{-{\sum }_j{\beta }_j{L}_j}\phi dw.\end{array}$$

The first term is $-{\mathbb{E}}_{p_{\beta }}[\varphi (w)L(w,z_i)]$. For the second, ${\partial }_{{\beta }_i}Z=-Z{\mathbb{E}}_{p_{\beta }}[L(w,z_i)]$, so it equals $+{\mathbb{E}}_{p_{\beta }}[L(w,z_i)]{\mathbb{E}}_{p_{\beta }}[\varphi (w)]$. Evaluating at $\beta =1$:

$$\mathrm{BIF}(z_i,\varphi )=-{\mathbb{E}}_p[\varphi \cdot L_i]+{\mathbb{E}}_p[\varphi ]{\mathbb{E}}_p[L_i]=-{\mathrm{Cov}}_{w\sim p}(L(w,z_i),\varphi (w)).$$

Structurally: in the classical IF, the "correlation" between $L_i$ and $\varphi$ is mediated by the inverse Hessian acting on their gradients at a point; in the BIF it is mediated by the full posterior distribution. The inverse Hessian is how a Gaussian posterior would produce a covariance (as we will see in Exercise 3.2), so the BIF strictly generalizes the classical formula.

Localization. Computing expectations over the global posterior $p(w\mid \mathcal{D})$ is generally intractable for neural networks, and we are typically interested in the behavior of a specific trained checkpoint $w^{\ast }$, not a global Bayesian average. Following Kreer et al. 2026, we localize the posterior by adding a Gaussian penalty centered at $w^{\ast }$:

where $\gamma >0$ is a localization strength that controls how tightly the posterior concentrates around $w^{\ast }$. The local BIF is then

where ${\mathrm{Cov}}_{\gamma }$ denotes covariance under $p_{\gamma }$. We will see in Exercise 3.2 that the localization parameter $\gamma$ plays exactly the role of the damping parameter $\lambda$ in the PBRF formula (6).

Practical computation via SGLD. The local BIF requires estimating a covariance under $p_{\gamma }$. This is done by drawing samples from $p_{\gamma }$ using stochastic gradient Langevin dynamics (SGLD), which adds calibrated Gaussian noise to SGD updates:

$$w_{t+1}=w_t-\frac{ϵ}{2}(\frac{n}{m}\sum _{k\in B_t}{\nabla }_{w}L(w_t,z_k)+\gamma (w_t-w^{\ast }))+{\eta }_t,{\eta }_t\sim \mathcal{N}(0,ϵI),$$

where $B_t$ is a mini-batch of size $m$ and $ϵ$ is the step size. After a burn-in period, the iterates $\{w_t\}$ are approximately distributed according to $p_{\gamma }$, and the covariance in (11) is estimated by the sample covariance of $(L(w_t,z_i),\varphi (w_t))$ across draws. Running multiple independent chains from $w^{\ast }$ improves coverage.

3.2 Exercise: Connecting Bayesian and classical influence functions

As we have seen in the SLT day, tools involving the Hessian are usually a low-order approximation of an expansion. This is also true for influence functions: In the non-degenerate case the classical IF is the leading-order term of the BIF under a Laplace approximation. The following exercise makes this precise by computing a power series expansion of the BIF covariance and identifying the classical IF as the leading-order term. The argument follows Kreer et al. 2026, Appendix A.

Exercise 3.2 (Power series expansion of Bayesian influence functions). Let $w^{\ast }$ be a local minimum of the training loss $L(w)={\sum }_{i=1}^{n}L(w,z_i)$ with positive-definite Hessian $H={\nabla }_w^{2}L(w^{\ast })$. Write $\Delta w=w-w^{\ast }$, and abbreviate $g_{\varphi }={\nabla }_w\varphi (w^{\ast })$, $H_{\varphi }={\nabla }_w^2\varphi (w^{\ast })$ for the gradient and Hessian of the observable, and $g_i={\nabla }_{w}L(w^{\ast },z_i)$, $H_i={\nabla }_w^{2}L(w^{\ast },z_i)$ for the per-sample loss.

(a) (Laplace approximation.) Consider the posterior $p(w\mid \mathcal{D})\propto e^{-L(w)}\phi (w)$. Taylor-expand $L(w)$ to second order around $w^{\ast }$, using the fact that $\nabla L(w^{\ast })=0$ at a minimum. Argue that for large $n$ the quadratic term in $L$ dominates the prior $\phi$, and conclude that

$$p(w\mid \mathcal{D})\approx \mathcal{N}(w^{\ast },H^{-1}).$$

Under this approximation, $\Delta w\sim \mathcal{N}(0,H^{-1})$. (This is called the Laplace approximation; the Bernstein–von Mises theorem guarantees it is asymptotically exact for non-singular models.)

(b) (Taylor expansion of the covariance.) Expand $\varphi (w)$ and $L(w,z_i)$ in Taylor series around $w^{\ast }$:

$$\begin{array}{rl}\varphi (w)& =\varphi (w^{\ast })+g_{\varphi }^{\top }\Delta w+\frac{1}{2}\Delta w^{\top }{H}_{\varphi }\Delta w+\cdots ,\\ L(w,z_i)& =L(w^{\ast },z_i)+g_i^{\top }\Delta w+\frac{1}{2}\Delta w^{\top }{H}_i\Delta w+\cdots .\end{array}$$

Using the fact that constant terms drop out of covariances and that $\mathrm{Cov}(X,Y)={\sum }_{k,m\ge 1}\mathrm{Cov}(T_k[\varphi ],T_m[L_i])$ where $T_k[f]$ denotes the $k$-th order term in the Taylor expansion of $f$, show that the BIF decomposes as

$$\mathrm{BIF}(z_i,\varphi )=-\sum _{\begin{array}{c}k,m\ge 1\\ k+m\text{ even}\end{array}}{\mathrm{Cov}}_{\mathcal{N}}(T_k[\varphi ],T_m[L_i]).$$

Why do terms with $k+m$ odd vanish?

(c) (Leading order: classical IF.) Compute the $(k,m)=(1,1)$ term:

$${\mathrm{Cov}}_{\mathcal{N}}(g_{\varphi }^{\top }\Delta w,g_i^{\top }\Delta w)=g_{\varphi }^{\top }{H}^{-1}{g}_i.$$

Conclude that $-g_{\varphi }^{\top }{H}^{-1}{g}_i=-{\nabla }_w\varphi (w^{\ast }{)}^{\top }{H}^{-1}{\nabla }_{w}L(w^{\ast },z_i)$ is the classical influence function $\mathcal{I}(z_i,\varphi )$ from Definition 2.2. This is the leading-order term of the BIF.

(d) (Second-order correction.) Compute the $(k,m)=(2,2)$ term using Isserlis' theorem (the Gaussian moment identity $\mathbb{E}[\Delta w_a\Delta w_b\Delta w_c\Delta w_d]={\Sigma }_{ab}{\Sigma }_{cd}+{\Sigma }_{ac}{\Sigma }_{bd}+{\Sigma }_{ad}{\Sigma }_{bc}$ with $\Sigma =H^{-1}$). Show that

$${\mathrm{Cov}}_{\mathcal{N}}\left(\frac{1}{2}\Delta w^{\top }{H}_{\varphi }\Delta w,\frac{1}{2}\Delta w^{\top }{H}_i\Delta w\right)=\frac{1}{2}\mathrm{tr}\left(H_{\varphi }{H}^{-1}{H}_i{H}^{-1}\right).$$

This correction involves the Hessians of the observable and per-sample loss. It captures second-order curvature interactions that the classical IF misses entirely. In the linear regression setting of Example 3.2, why does this correction vanish?

(e) (Localized version: damped IF.) Now consider the local BIF of Equation 11 with localization strength $\gamma$. The localized posterior is approximately $\mathcal{N}(w^{\ast },(H+\gamma I{)}^{-1})$. Repeat the leading-order computation of part (c) to show that

$${\mathrm{BIF}}_{\gamma }(z_i,\varphi )\approx -{\nabla }_w\varphi (w^{\ast }{)}^{\top }(H+\gamma I{)}^{-1}{\nabla }_{w}L(w^{\ast },z_i).$$

This is precisely the damped influence function of Equation 6, with the localization strength $\gamma$ playing the role of the damping parameter $\lambda$. The local BIF is thus a natural, higher-order generalization of the damped IF: it agrees at leading order and includes all the corrections from parts (b)–(d) with $H^{-1}$ replaced by $(H+\gamma I{)}^{-1}$.

Solution

(a) The log-posterior is $\log p(w\mid \mathcal{D})=-L(w)+\log \phi (w)+\mathrm{const}$. Taylor-expanding $L(w)$ around the minimum $w^{\ast }$:

$$L(w)=L(w^{\ast })+\underset{=0}{\underbrace{\nabla L(w^{\ast }{)}^{\top }}}\Delta w+\frac{1}{2}\Delta w^{\top }H\Delta w+O(\Vert \Delta w{\Vert }^3).$$

The linear term vanishes because $w^{\ast }$ is a critical point. Substituting into $\log p$:

$$\log p(w\mid \mathcal{D})\approx -L(w^{\ast })-\frac{1}{2}\Delta w^{\top }H\Delta w+\log \phi (w)+\mathrm{const}.$$

Since $L={\sum }_{i=1}^n{L}_i$, the Hessian $H={\sum }_i{\nabla }^2{L}_i$ scales as $O(n)$, while the prior contributes $O(1)$ to the log-density. For large $n$ the quadratic term dominates, giving $p(w\mid \mathcal{D})\approx \mathcal{N}(w^{\ast },H^{-1})$. This is the Laplace approximation. (The Bernstein–von Mises theorem makes this rigorous: under regularity, the posterior converges to this Gaussian in total variation as $n\to \infty$.)

(b) Since $\mathrm{Cov}(X+c,Y)=\mathrm{Cov}(X,Y)$ for any constant $c$, the constant terms $\varphi (w^{\ast })$ and $L(w^{\ast },z_i)$ drop out. The covariance of the two Taylor series is bilinear, so it distributes over the sum of terms:

$$\mathrm{Cov}(\varphi ,L_i)=\sum _{k\ge 1}\sum _{m\ge 1}\mathrm{Cov}(T_k[\varphi ],T_m[L_i]).$$

Under $\Delta w\sim \mathcal{N}(0,H^{-1})$, $T_k[\varphi ]$ is a degree-$k$ polynomial in $\Delta w$. The covariance $\mathrm{Cov}(T_k,T_m)=\mathbb{E}[T_k{T}_m]-\mathbb{E}[T_k]\mathbb{E}[T_m]$ involves moments of $\Delta w$ of degree $k+m$. For a centered Gaussian, odd moments vanish: $\mathbb{E}[\Delta w_{a_1}\cdots \Delta w_{a_{\ell }}]=0$ when $\ell$ is odd. If $k+m$ is odd, then all moments in $\mathbb{E}[T_k{T}_m]$ and $\mathbb{E}[T_k]\mathbb{E}[T_m]$ involve an odd total degree, so both vanish and $\mathrm{Cov}(T_k,T_m)=0$.

(c) $T_1[\varphi ]=g_{\varphi }^{\top }\Delta w$ and $T_1[L_i]=g_i^{\top }\Delta w$ are linear in $\Delta w$. Their covariance is

$$\mathrm{Cov}(g_{\varphi }^{\top }\Delta w,g_i^{\top }\Delta w)=g_{\varphi }^{\top }\mathbb{E}[\Delta w\Delta w^{\top }]g_i=g_{\varphi }^{\top }{H}^{-1}{g}_i.$$

Negating gives $-g_{\varphi }^{\top }{H}^{-1}{g}_i=-\nabla \varphi (w^{\ast }{)}^{\top }{H}^{-1}\nabla L(w^{\ast },z_i)$, which is the classical influence function $\mathcal{I}(z_i,\varphi )$ from Definition 2.2.

(d) Write $Q_{\varphi }=\frac{1}{2}\Delta w^{\top }{H}_{\varphi }\Delta w$ and $Q_i=\frac{1}{2}\Delta w^{\top }{H}_i\Delta w$. We need $\mathrm{Cov}(Q_{\varphi },Q_i)=\mathbb{E}[Q_{\varphi }{Q}_i]-\mathbb{E}[Q_{\varphi }]\mathbb{E}[Q_i]$ with $\Delta w\sim \mathcal{N}(0,\Sigma )$, $\Sigma =H^{-1}$.

For the expectation: $\mathbb{E}[Q_{\varphi }]=\frac{1}{2}\mathrm{tr}(H_{\varphi }\Sigma )$ and $\mathbb{E}[Q_i]=\frac{1}{2}\mathrm{tr}(H_i\Sigma )$.

For the cross-moment: $\mathbb{E}[Q_{\varphi }{Q}_i]=\frac{1}{4}{\sum }_{a,b,c,d}(H_{\varphi }{)}_{ab}(H_i{)}_{cd}\mathbb{E}[\Delta w_a\Delta w_b\Delta w_c\Delta w_d]$. By Isserlis' theorem:

$$\mathbb{E}[\Delta w_a\Delta w_b\Delta w_c\Delta w_d]={\Sigma }_{ab}{\Sigma }_{cd}+{\Sigma }_{ac}{\Sigma }_{bd}+{\Sigma }_{ad}{\Sigma }_{bc}.$$

The first pairing gives $\frac{1}{4}\mathrm{tr}(H_{\varphi }\Sigma )\mathrm{tr}(H_i\Sigma )=\mathbb{E}[Q_{\varphi }]\mathbb{E}[Q_i]$, which cancels in the covariance. The other two pairings each give $\frac{1}{4}\mathrm{tr}(H_{\varphi }\Sigma H_i\Sigma )$ (by cyclicity of the trace and symmetry of $H_{\varphi }$, $H_i$, and $\Sigma$). So

$$\mathrm{Cov}(Q_{\varphi },Q_i)=2\cdot \frac{1}{4}\mathrm{tr}(H_{\varphi }\Sigma H_i\Sigma )=\frac{1}{2}\mathrm{tr}(H_{\varphi }{H}^{-1}{H}_i{H}^{-1}).$$

In the linear regression setting with $\varphi (w)=w_k$ (a coordinate function), $H_{\varphi }={\nabla }^2{w}_k=0$: the observable is linear in $w$, so its Hessian vanishes and the entire $(2,2)$ correction is zero. This is why the BIF equals the (damped) IF exactly for linear regression — all corrections beyond leading order involve $H_{\varphi }$ or higher derivatives of $\varphi$, which vanish for a linear observable.

(e) The localized posterior of Equation 10 has effective Hessian $H_{\mathrm{eff}}=H+\gamma I$, so the Laplace approximation gives $\Delta w\sim \mathcal{N}(0,(H+\gamma I{)}^{-1})$. The leading-order $(1,1)$ computation from part (c) becomes

$${\mathrm{BIF}}_{\gamma }(z_i,\varphi )\approx -g_{\varphi }^{\top }(H+\gamma I{)}^{-1}{g}_i=-\nabla \varphi (w^{\ast }{)}^{\top }(H+\gamma I{)}^{-1}\nabla L(w^{\ast },z_i),$$

which is the damped IF of Equation 6 with $\gamma$ in the role of $\lambda$. The higher-order corrections from parts (b)–(d) carry through with $H^{-1}$ replaced by $(H+\gamma I{)}^{-1}$ throughout.

(f) Summary: For regular (non-singular) models, the posterior is approximately Gaussian (Bernstein–von Mises), the Taylor expansion converges, and the $(1,1)$ term dominates (it scales as $O(1/n)$ while higher terms scale as $O(1/n^2)$ and beyond). The classical IF is a good approximation because it is the leading term. For singular models (neural networks), the posterior is concentrated on a positive-dimensional variety, not at an isolated point. The Laplace approximation fails: the Hessian has a large null space, the posterior is non-Gaussian, and the Taylor series does not converge around $w^{\ast }$. In this regime, the BIF — defined as an exact covariance under the true posterior — captures the full geometry, while the classical IF is at best the leading term of an expansion that does not converge. The BIF is the fundamental object; the IF is the Gaussian shadow it casts.

3.3 Long Exercise: Influence functions as optimal linear transport

The previous exercise showed that the classical IF is the leading-order term of the BIF under a Laplace approximation. We now take a different perspective, due to Mlodozeniec et al. 2025: instead of asking "what does the BIF reduce to?" we ask "what is the best way to shift the parameters to approximate the effect of a data perturbation?" The answer turns out to be: the influence function, and its optimality does not require the Laplace approximation to hold.

Setup. Let $L(w)={\sum }_{i=1}^{n}L(w,z_i)$ denote the total training loss, and abbreviate $L_k(w):=L(w,z_k)$ for the per-sample loss at example $k$. Define the Boltzmann distributions at inverse temperature $\beta$:

$$\begin{array}{rlrl}{p}_0^{\beta }(w)& \propto e^{-\beta L(w)},& p_{\epsilon }^{\beta }(w)& \propto e^{-\beta (L(w)-\epsilon L_k(w))},\end{array}$$

so that $p_{\epsilon }^{\beta }$ is the Boltzmann distribution after downweighting sample $k$ by $\epsilon$. The low-temperature limit $\beta \to \infty$ concentrates $p_0^{\beta }$ on the set of global minima ${\mathcal{S}}_L=\{w:L(w)=\inf L\}$.

Now suppose we want to approximate the effect of the perturbation without resampling from $p_{\epsilon }^{\beta }$. The simplest approach is to take each parameter $w$ drawn from $p_0^{\beta }$ and shift it to $w+\epsilon r(w)$ for some smooth vector field $r:{\mathbb{R}}^d\to {\mathbb{R}}^d$. If $r$ is chosen well, the distribution of the shifted parameters should be close to $p_{\epsilon }^{\beta }$. The question is: what is the optimal $r$?

To make this precise, we need to know the density of the shifted parameters. If $w\sim p_0^{\beta }$ and we set $\tilde{w}=w+\epsilon r(w)$, then by the multivariate change-of-variables formula, $\tilde{w}$ has density

where $\nabla r$ is the $d\times d$ Jacobian matrix of $r$ and the denominator is the usual volume-change factor.⁴ We measure how well $q_{\epsilon }$ approximates $p_{\epsilon }^{\beta }$ by the KL divergence $D_{\mathrm{KL}}(q_{\epsilon }\Vert p_{\epsilon }^{\beta })$. Rather than computing this directly, it is cleaner to write it as an expectation over the original distribution $p_0^{\beta }$:

This identity follows by substituting (12) into the definition of KL and changing variables back to $w$. We define the asymptotic KL cost

The $1/\beta$ normalization kills the determinant term (which is $O(1)$) and the entropy terms, isolating the energy contribution (which is $O(\beta )$).

Exercise 3.3 (Influence functions as optimal parameter shifts). Let $w^{\ast }$ be a non-degenerate minimum of $L$ with positive-definite Hessian $H={\nabla }^{2}L(w^{\ast })$. In this exercise we work at the unique minimum, so the $\beta \to \infty$ limit simply evaluates everything at $w^{\ast }$.

(a) (KL in terms of energies.) Starting from (13), substitute $p_0^{\beta }(w)=e^{-\beta L(w)}/Z(\beta ,0)$ and $p_{\epsilon }^{\beta }(w)=e^{-\beta (L(w)-\epsilon L_k(w))}/Z(\beta ,\epsilon )$, and show that

$$\begin{array}{rl}\frac{1}{\beta }{D}_{\mathrm{KL}}=& {\mathbb{E}}_{p_0^{\beta }}[L(w+\epsilon r)-\epsilon L_k(w+\epsilon r)-L(w)]\\ & +\frac{1}{\beta }\log \frac{Z(\beta ,\epsilon )}{Z(\beta ,0)}-\frac{1}{\beta }{\mathbb{E}}_{p_0^{\beta }}[\log |\det (I+\epsilon \nabla r)|],\end{array}$$

where $Z(\beta ,\epsilon )=\int e^{-\beta (L(w)-\epsilon L_k(w))}dw$. Explain why the last two terms vanish in the $\beta \to \infty$ limit. Hint: $\log |\det (I+\epsilon \nabla r)|=O(1)$, and by Laplace's method $\frac{1}{\beta }\log Z(\beta ,\epsilon )\to -{\inf }_w[L(w)-\epsilon L_k(w)]$.

(b) (Taylor expansion at $w^{\ast }$.) Since the expectation under $p_0^{\beta }$ concentrates at $w^{\ast }$ as $\beta \to \infty$, expand $L(w^{\ast }+\epsilon r)$ and $\epsilon L_k(w^{\ast }+\epsilon r)$ using $\nabla L(w^{\ast })=0$:

$$\begin{array}{rl}L(w^{\ast }+\epsilon r)& =L(w^{\ast })+\frac{{\epsilon }^2}{2}{r}^{\top }Hr+O({\epsilon }^3),\\ \epsilon L_k(w^{\ast }+\epsilon r)& =\epsilon L_k(w^{\ast })+{\epsilon }^2\nabla L_k(w^{\ast }{)}^{\top }r+O({\epsilon }^3).\end{array}$$

Similarly, show that ${\inf }_w[L(w)-\epsilon L_k(w)]=L(w^{\ast })-\epsilon L_k(w^{\ast })-\frac{{\epsilon }^2}{2}\nabla L_k^{\top }{H}^{-1}\nabla L_k+O({\epsilon }^3)$, by minimizing the Taylor expansion of $L-\epsilon L_k$ around $w^{\ast }$. Combine everything to obtain

$$\mathcal{F}(r,\epsilon )=\frac{{\epsilon }^2}{2}(r-H^{-1}\nabla L_k{)}^{\top }H(r-H^{-1}\nabla L_k)+O({\epsilon }^3).$$

Hint: After substitution, the terms that do not depend on $r$ should assemble into $\frac{{\epsilon }^2}{2}\nabla L_k^{\top }{H}^{-1}\nabla L_k$ and cancel with the partition-function contribution, leaving a perfect square.

(c) (IFs are optimal.) The cost (15) is a squared Mahalanobis distance (in the Hessian metric) between the shift $r$ and the influence function $r_{\mathrm{IF}}:=H^{-1}\nabla L_k(w^{\ast })$. Conclude that the IF minimizes $\mathcal{F}$ to $O({\epsilon }^2)$ among all smooth shift maps $w\mapsto w+\epsilon r(w)$, and that the minimum cost is $O({\epsilon }^3)$.

State the result in words: the influence function is the approximately KL-optimal way to shift the parameters of the unperturbed Boltzmann distribution to approximate the perturbed one.

(d) (Degenerate case.) Now suppose $w^{\ast }$ lies on a minimum manifold ${\mathcal{S}}_L$ and $H$ is singular with null space $N$. Argue (without detailed proof) that:

• The cost (15) generalizes to $\mathcal{F}(r,\epsilon )=\frac{{\epsilon }^2}{2}(r-H^{+}\nabla L_k{)}^{\top }H(r-H^{+}\nabla L_k)+O({\epsilon }^3)$, where $H^{+}$ is the Moore–Penrose pseudoinverse.
• The minimizer is $r=H^{+}\nabla L_k+r_{\mathrm{NS}}$ for any $r_{\mathrm{NS}}\in N$: directions in the null space of $H$ have zero cost, because the Hessian assigns zero curvature to movement along the minimum manifold.
• The IF (with pseudoinverse) is therefore optimal for the directions that matter, but the component of the shift along the flat directions is unconstrained.

Connect this to the SLT perspective: for singular models, the minimum is a variety, and the IF tells you how to move off the variety but is silent about movement along it.

(e) (Linear regression check.) Specialize to the Bayesian linear regression setting of Example 3.2 with prior $w\sim \mathcal{N}(0,{\tau }^{2}I)$. The Boltzmann distribution at inverse temperature $\beta =1$ is the Bayesian posterior $\mathcal{N}(\mu ,V)$, which is already Gaussian, no need to take $\beta \to \infty$. A constant shift $w\mapsto w+\epsilon r$ sends $\mathcal{N}(\mu ,V)$ to $\mathcal{N}(\mu +\epsilon r,V)$. The true perturbed posterior (downweighting $z_k$ by $\epsilon$) has mean $\mu +\epsilon \mathrm{BIF}(z_k,w)+O({\epsilon }^2)$ and covariance $V+O(\epsilon )$. Using the KL formula for Gaussians with the same covariance, show that

$$D_{\mathrm{KL}}(\mathcal{N}(\mu +\epsilon r,V)\Vert \mathcal{N}(\mu +\epsilon b,V))=\frac{{\epsilon }^2}{2}(r-b{)}^{\top }{V}^{-1}(r-b),$$

which is minimized at $r=b=\mathrm{BIF}(z_k,w)=(A+\lambda I{)}^{-1}{x}_k{r}_k$. Verify that this is consistent with (15): the "Hessian of the energy" (training loss plus prior) is $V^{-1}={\sigma }^{-2}A+{\tau }^{-2}I$, and $V^{-1}\cdot \mathrm{BIF}={\sigma }^{-2}{x}_k{r}_k={\nabla }_{w}L(w,z_k){|}_{w=\mu }$.

(f) (Comparison with Exercise 3.2.) The Laplace expansion exercise and this exercise both relate the IF to the BIF, but they answer different questions:

• Exercise 3.2 is an algebraic identity: BIF $=$ IF $+$ higher-order corrections. It requires the Laplace approximation (hence non-singular $H$, Gaussian posterior) and tells you what happens when you truncate the BIF.
• This exercise is an optimality result: the IF shift minimizes the KL divergence between the shifted and true perturbed distributions. The $\beta \to \infty$ result (parts a–d) needs only mild regularity of $\mathcal{L}$ and works even for singular $H$ (via the pseudoinverse).

Summarize: the Laplace expansion tells you what the IF is (the leading Gaussian term of the BIF). The optimality perspective tells you why it works (it is the best first-order correction to the parameter distribution). The second viewpoint does not require the posterior to be Gaussian, which is why Mlodozeniec et al. 2025 argue that it provides a better explanation for the empirical success of IFs in deep learning.

Solution

(a) Using $p_0^{\beta }(w)=e^{-\beta L(w)}/Z(\beta ,0)$ and $p_{\epsilon }^{\beta }(w)=e^{-\beta (L(w)-\epsilon L_k(w))}/Z(\beta ,\epsilon )$, substituting into (13) gives

$$\begin{array}{rl}{D}_{\mathrm{KL}}(q_{\epsilon }\Vert p_{\epsilon }^{\beta })& ={\mathbb{E}}_{w\sim p_0^{\beta }}[\log p_0^{\beta }(w)-\log p_{\epsilon }^{\beta }(w+\epsilon r)-\log |\det (I+\epsilon \nabla r)|]\\ & ={\mathbb{E}}_{p_0^{\beta }}[-\beta L(w)-\log Z(\beta ,0)+\beta (L(w+\epsilon r)-\epsilon L_k(w+\epsilon r))\\ & +\log Z(\beta ,\epsilon )-\log |\det (I+\epsilon \nabla r)|].\end{array}$$

Dividing by $\beta$ and rearranging gives the claimed expression. As $\beta \to \infty$: (i) $\frac{1}{\beta }\log |\det (I+\epsilon \nabla r)|\to 0$ since the log-determinant is $O(1)$ (bounded for smooth $r$ and small $\epsilon$); (ii) by Laplace's method, $\frac{1}{\beta }\log Z(\beta ,\epsilon )\to -{\inf }_w[L(w)-\epsilon L_k(w)]$, so $\frac{1}{\beta }\log \frac{Z(\beta ,\epsilon )}{Z(\beta ,0)}\to \inf L-\inf (L-\epsilon L_k)$, which is a constant independent of $r$; (iii) the expectation concentrates on $w^{\ast }$.

(b) At $w^{\ast }$, $\nabla L(w^{\ast })=0$, so the Taylor expansion of $L(w^{\ast }+\epsilon r)$ starts at the quadratic term: $L(w^{\ast })+\frac{{\epsilon }^2}{2}{r}^{\top }Hr+O({\epsilon }^3)$. For $\epsilon L_k$: $\epsilon L_k(w^{\ast }+\epsilon r)=\epsilon L_k(w^{\ast })+{\epsilon }^2\nabla L_k^{\top }r+O({\epsilon }^3)$. Combined:

$$L(w^{\ast }+\epsilon r)-\epsilon L_k(w^{\ast }+\epsilon r)-L(w^{\ast })=-\epsilon L_k(w^{\ast })+\frac{{\epsilon }^2}{2}{r}^{\top }Hr-{\epsilon }^2\nabla L_k^{\top }r+O({\epsilon }^3).$$

For the infimum: expanding $L(w^{\ast }+\delta )-\epsilon L_k(w^{\ast }+\delta )$ to second order in $\delta$ and minimizing, the optimal perturbation is ${\delta }^{\ast }=\epsilon H^{-1}\nabla L_k$, giving

$$\underset{w}{\inf }[L-\epsilon L_k]=L(w^{\ast })-\epsilon L_k(w^{\ast })-\frac{{\epsilon }^2}{2}\nabla L_k^{\top }{H}^{-1}\nabla L_k+O({\epsilon }^3).$$

So $\inf L-\inf (L-\epsilon L_k)=\epsilon L_k(w^{\ast })+\frac{{\epsilon }^2}{2}\nabla L_k^{\top }{H}^{-1}\nabla L_k+O({\epsilon }^3)$.

Adding the two contributions:

$$\begin{array}{rl}\mathcal{F}(r,\epsilon )& =\frac{{\epsilon }^2}{2}{r}^{\top }Hr-{\epsilon }^2\nabla L_k^{\top }r+\frac{{\epsilon }^2}{2}\nabla L_k^{\top }{H}^{-1}\nabla L_k+O({\epsilon }^3)\\ & =\frac{{\epsilon }^2}{2}(r^{\top }Hr-2\nabla L_k^{\top }r+\nabla L_k^{\top }{H}^{-1}\nabla L_k)+O({\epsilon }^3).\end{array}$$

Completing the square: $r^{\top }Hr-2\nabla L_k^{\top }r+\nabla L_k^{\top }{H}^{-1}\nabla L_k=(r-H^{-1}\nabla L_k{)}^{\top }H(r-H^{-1}\nabla L_k)$, since expanding the right side gives $r^{\top }Hr-2\nabla L_k^{\top }{H}^{-1}\cdot Hr+\nabla L_k^{\top }{H}^{-1}H{H}^{-1}\nabla L_k=r^{\top }Hr-2\nabla L_k^{\top }r+\nabla L_k^{\top }{H}^{-1}\nabla L_k$.

(c) Since $H$ is positive definite, the squared Mahalanobis form $(r-r_{\mathrm{IF}}{)}^{\top }H(r-r_{\mathrm{IF}})\ge 0$ with equality iff $r=r_{\mathrm{IF}}=H^{-1}\nabla L_k$. At the minimum, $\mathcal{F}(r_{\mathrm{IF}},\epsilon )=O({\epsilon }^3)$: the IF shifts the distribution with only a third-order KL residual. Any other shift $r\ne r_{\mathrm{IF}}$ incurs a strictly positive $O({\epsilon }^2)$ cost proportional to its Mahalanobis distance from the IF. The IF is the unique optimal shift among all smooth vector fields $r$, not just among constant maps.

(d) When $H$ is singular, the quadratic form $r^{\top }Hr$ is zero for $r\in N=\mathrm{Null}(H)$: movement in the null space has zero Hessian cost because the loss is flat along the minimum manifold. The pseudoinverse $H^{+}$ inverts only the non-degenerate directions, and the completed square becomes $(r-H^{+}\nabla L_k{)}^{\top }H(r-H^{+}\nabla L_k)$. This vanishes whenever $r-H^{+}\nabla L_k\in N$, i.e., the minimizer is $r=H^{+}\nabla L_k+r_{\mathrm{NS}}$ for any $r_{\mathrm{NS}}\in N$.

Geometrically: the IF determines the optimal shift off the minimum manifold (perpendicular to ${\mathcal{S}}_L$), but movement along ${\mathcal{S}}_L$ is invisible to the energy at this order. In the SLT language: the degeneracy of the singular set means that many shifts are equally good at leading order — the IF picks one, but the null-space ambiguity is real.

(e) For linear regression with Gaussian prior, the posterior is $\mathcal{N}(\mu ,V)$ with $V^{-1}={\sigma }^{-2}A+{\tau }^{-2}I$. A constant shift $w\mapsto w+\epsilon r$ sends $\mathcal{N}(\mu ,V)$ to $\mathcal{N}(\mu +\epsilon r,V)$. The true perturbed mean is $\mu +\epsilon b+O({\epsilon }^2)$ with $b=\mathrm{BIF}(z_k,w)=(A+\lambda I{)}^{-1}{x}_k{r}_k$. For two Gaussians with the same covariance:

$$D_{\mathrm{KL}}(\mathcal{N}(\mu +\epsilon r,V)\Vert \mathcal{N}(\mu +\epsilon b,V))=\frac{{\epsilon }^2}{2}(r-b{)}^{\top }{V}^{-1}(r-b).$$

This is minimized at $r=b=(A+\lambda I{)}^{-1}{x}_k{r}_k$. To check consistency with (15): the "Hessian of the energy" in the Boltzmann distribution $p\propto e^{-(\text{loss}+\text{prior})}$ is $H_{\mathrm{eff}}={\sigma }^{-2}A+{\tau }^{-2}I=V^{-1}$, and $H_{\mathrm{eff}}^{-1}{\nabla }_{w}L(\mu ,z_k)=V\cdot {\sigma }^{-2}{x}_k{r}_k=(A+\lambda I{)}^{-1}{x}_k{r}_k=\mathrm{BIF}$. The general formula (15) reduces to $\frac{{\epsilon }^2}{2}(r-b{)}^{\top }{V}^{-1}(r-b)$, exactly as computed directly.

(f) The two exercises are complementary:

• The Laplace expansion (Exercise 3.2) tells you what the IF is: the leading-order term of the BIF covariance under a Gaussian approximation to the posterior. The higher-order corrections (Hessian products, etc.) are visible but require the Laplace approximation to hold.
• The transport perspective (this exercise) tells you why the IF works even when the Laplace approximation fails: it minimizes the KL cost of transporting the parameter distribution, and this optimality holds in the $\beta \to \infty$ limit under mild regularity assumptions — no Gaussianity required. For singular models, the pseudoinverse handles the degenerate directions naturally (part d), and the result still says: the IF is the best you can do with a smooth first-order map.

This provides what Mlodozeniec et al. 2025 call a "distributional" justification for IFs in deep learning. The original justification (implicit function theorem, convexity) does not hold for neural networks. The Laplace/BIF justification (Exercise 3.2) requires Gaussianity of the posterior, which also fails. But the transport optimality result asks only for bounded derivatives and a well-defined minimum — conditions much closer to what actually holds in practice.

4. Unrolling

The influence functions of Section 2 and the Bayesian influence functions of Section 3 both have something in common: they analyze the endpoint of training, not the dynamics that led to it. Under the assumptions of the unique global minimum, this is not a problem: the training trajectory is irrelevant, since all roads lead to the same outcome. But in practice the training path will matter a lot. Ultimately, we are studying generalization and this is ultimately a question about the training dynamics. Similarly, the BIF looks at the posterior expectation shifts. To give a concrete failure mode: a training point that appears in the first mini-batch of the first epoch will have a very different effect compared to the same point appearing towards the end of training.

If we write down the training trajectory as a composition of $\beta$-dependent maps:

$$\begin{array}{ccccccccc}{\mathbb{R}}^p& \stackrel{h_0(\cdot ;\beta )}{\to }& {\mathbb{R}}^p& \stackrel{h_1(\cdot ;\beta )}{\to }& \cdots & \stackrel{h_{T-1}(\cdot ;\beta )}{\to }& {\mathbb{R}}^p& \stackrel{\varphi }{\to }& \mathbb{R}\\ w_0& ⟼& w_1& ⟼& \cdots & ⟼& w_T& ⟼& \varphi (w_T)\end{array}$$

where each $h_t$ is one optimizer step (a gradient step, an Adam update, etc.): $w_{t+1}=h_t(w_t;\beta )$, we can observe an analogy to a forward pass through a neural network:

$$\begin{array}{ccccccccc}{\mathbb{R}}^{d_0}& \stackrel{f_1(\cdot ;w)}{\to }& {\mathbb{R}}^{d_1}& \stackrel{f_2(\cdot ;w)}{\to }& \cdots & \stackrel{f_L(\cdot ;w)}{\to }& {\mathbb{R}}^{d_L}& \stackrel{\ell }{\to }& \mathbb{R}\\ a_0& ⟼& a_1& ⟼& \cdots & ⟼& a_L& ⟼& \ell (a_L)\end{array}$$

where each $f_{\ell }$ is one layer: $a_{\ell }=f_{\ell }(a_{\ell -1};w)$.

	**Forward pass **	Training trajectory
Inputs	Model weights $w$	Data weights $\beta$
Intermediate states	Activations $a_1,\dots ,a_L$	Parameters $w_0,\dots ,w_T$
Output	Loss $\ell (\text{logits})$	Observable $\varphi (w_T)$
Differentiation method	Backpropagation	Unrolling

Data attribution via unrolling is to a training run what backprop is to a forward pass: differentiating a composition of maps to identify which inputs matter most.

This invites the following question: Can we backprop through the training trajectory? The gradients with respect to the ${\beta }_i$ then indicate which data points would need to be up- or downweighted to change the final observable $\varphi (w_T)$. And in principle, we can! We can differentiate through the training computation, step by step, treating the final parameters as an explicit function of every gradient update along the way. The result is a chain-rule expression that decomposes data attribution into per-step contributions and naturally incorporates the optimizer, learning rate schedule, mini-batch ordering, and training duration. We will derive the unrolling formula similar to the backprop formula for a forward pass, and then discuss two approaches to making it practical: one that materializes Jacobians and one that computes Jacobian-vector products implicitly. In the setup of Equation 1 we take $\pi (\beta )=w_T(\beta )$, the final checkpoint of the training run.

4.1 The training-dynamics approach to attribution

We fix notation for the training process. Given a dataset $\mathcal{D}=\{z_1,\dots ,z_n\}$ with sample weights $\beta =({\beta }_1,\dots ,{\beta }_n)$, the optimizer runs for $T$ steps. At step $t$, a mini-batch $B_t\subseteq \{1,\dots ,n\}$ of size $b$ is drawn, and the parameters are updated by

where ${\eta }_t$ is the learning rate at step $t$. At $\beta =1$ this is ordinary SGD on the unweighted loss. The initial parameters $w_0$ are independent of $\beta$. For simplicity we can assume that we are training for one epoch only, so that each example appears in exactly one mini-batch, but the formulas hold more generally.

The counterfactual. We want

$$\frac{\partial \varphi (w_T(\beta ))}{\partial {\beta }_i}{|}_{\beta =1},$$

the sensitivity of a final-model observable $\varphi$ to an infinitesimal change in the weight of training example $z_i$.

4.2 The unrolling formula

Central to computing the data weight gradient is the Jacobian of one SGD step with respect to the current parameters. Define

where ${\hat{H}}_t:=\frac{1}{b}{\sum }_{i\in B_t}{\nabla }_w^{2}L(w_t,z_i)$ is the mini-batch Hessian at step $t$. The Jacobian of the map from step $t$ to step $t^{\prime }$ is the ordered product

Similarly, the direct effect of ${\beta }_i$ at step $t$ (when $i\in B_t$) is the per-example gradient:

The full derivative is obtained by chaining direct effects through remaining steps:

Each term in the sum corresponds to one training step in which $z_i$ appeared in the mini-batch: it is the gradient of $z_i$ at the parameters of that step, propagated forward through the Jacobians of all subsequent steps.

The influence on an observable $\varphi$ follows by the chain rule:

Exercise 4.1 (Derivation of the unrolling formula). (a) Starting from the SGD update (16), apply the chain rule to write $\frac{\partial w_{t+1}}{\partial {\beta }_i}{|}_{\beta =1}$ as a sum of two terms: one from the explicit dependence on ${\beta }_i$ (the direct effect) and one from the dependence of $w_t$ on ${\beta }_i$ (the indirect effect). Show that the result is the recursion

$$\frac{\partial w_{t+1}}{\partial {\beta }_i}=J_t\frac{\partial w_t}{\partial {\beta }_i}-\frac{{\eta }_t}{b}1[i\in B_t]{\nabla }_{w}L(w_t,z_i),$$

with initial condition $\frac{\partial w_0}{\partial {\beta }_i}=0$. Compute the Jacobian $J_t$ explicitly in terms of the mini-batch Hessian ${\hat{H}}_t$ which leads to (17).

(b) Solve the recursion by "unrolling" it (i.e. by substituting repeatedly and using (18)) to obtain (20).

(c) (Preconditioned SGD.) Suppose instead that the optimizer uses a preconditioner $P_t\succ 0$:

$$w_{t+1}(\beta )=w_t(\beta )-\frac{{\eta }_t}{b}{P}_t\sum _{i\in B_t}{\beta }_i{\nabla }_{w}L(w_t(\beta ),z_i).$$

This includes momentum-free Adam and natural gradient methods as special cases (with appropriate $P_t$). Show that the step Jacobian becomes $J_t^{(P)}=I-{\eta }_t{P}_t{\hat{H}}_t$ and that the unrolling formula generalizes to

$$\frac{\partial w_T(\beta )}{\partial {\beta }_i}{|}_{\beta =1}=-\sum _{t=0}^{T-1}\frac{{\eta }_t}{b}1[i\in B_t]J_{(t+1):T}^{(P)}{P}_t{\nabla }_{w}L(w_t,z_i),$$

where $J_{(t+1):T}^{(P)}={\prod }_{s=t+1}^{T-1}(I-{\eta }_s{P}_s{\hat{H}}_s)$. The only change is that each gradient is premultiplied by the preconditioner at the step where it appears. What does this say about the effect of using Adam vs. SGD on the attribution of a data point that appears early in training?

Solution

(a) Differentiating (16) with respect to ${\beta }_i$ at $\beta =1$ via the chain rule:

$$\begin{array}{rl}\frac{\partial w_{t+1}}{\partial {\beta }_i}& =\frac{\partial w_t}{\partial {\beta }_i}-\frac{{\eta }_t}{b}\sum _{j\in B_t}[\underset{\text{indirect}}{\underbrace{{\nabla }_w^{2}L(w_t,z_j)\frac{\partial w_t}{\partial {\beta }_i}}}+\underset{\text{direct}}{\underbrace{{\delta }_{ij}{\nabla }_{w}L(w_t,z_j)}}]\\ & =(I-\frac{{\eta }_t}{b}\sum _{j\in B_t}{\nabla }_w^{2}L(w_t,z_j))\frac{\partial w_t}{\partial {\beta }_i}-\frac{{\eta }_t}{b}1[i\in B_t]{\nabla }_{w}L(w_t,z_i)\\ & =J_t\frac{\partial w_t}{\partial {\beta }_i}-\frac{{\eta }_t}{b}1[i\in B_t]{\nabla }_{w}L(w_t,z_i).\end{array}$$

The initial condition is $\partial w_0/\partial {\beta }_i=0$ since $w_0$ does not depend on $\beta$.

(b) Substituting the recursion repeatedly:

$$\begin{array}{rl}\frac{\partial w_T}{\partial {\beta }_i}& =J_{T-1}\frac{\partial w_{T-1}}{\partial {\beta }_i}-\frac{{\eta }_{T-1}}{b}1[i\in B_{T-1}]{\nabla }_{w}L(w_{T-1},z_i)\\ & =J_{T-1}(J_{T-2}\frac{\partial w_{T-2}}{\partial {\beta }_i}-\frac{{\eta }_{T-2}}{b}1[i\in B_{T-2}]{\nabla }_{w}L(w_{T-2},z_i))-\frac{{\eta }_{T-1}}{b}1[i\in B_{T-1}]{\nabla }_{w}L(w_{T-1},z_i)\\ & =\cdots =-\sum _{t=0}^{T-1}\frac{{\eta }_t}{b}1[i\in B_t]J_{(t+1):T}{\nabla }_{w}L(w_t,z_i),\end{array}$$

where the last step follows from $\partial w_0/\partial {\beta }_i=0$ killing the $J_{0:T}$ term. This is (20).

(c) With the preconditioned update (22), the indirect effect picks up the preconditioner in the Hessian term: $J_t^{(P)}=I-({\eta }_t/b){\sum }_{j\in B_t}{P}_t{\nabla }_w^{2}L(w_t,z_j)=I-{\eta }_t{P}_t{\hat{H}}_t$. The direct effect becomes $-({\eta }_t/b)1[i\in B_t]P_t{\nabla }_{w}L(w_t,z_i)$. Unrolling the recursion as in (b) gives (23). The preconditioner $P_t$ at step $t$ acts as a local rescaling of the gradient: Adam's adaptive scaling amplifies gradients in directions with historically small second moments. A datum appearing early in training (when Adam has not yet accumulated accurate statistics) will have its gradient rescaled differently than the same datum appearing later, when $P_t$ has stabilized. This is a concrete mechanism by which optimizer choice affects per-datum attribution.

There are two ways to make (20) practical. The first is to compute the Jacobian products $J_{(t+1):T}$ explicitly, or rather an approximation thereof. The second is to avoid materializing Jacobians altogether and instead compute Jacobian-vector products implicitly by reverse-mode automatic differentiation through the training loop. We will discuss both approaches in turn.

4.2.1 Materializing the Jacobian

In this subsection we will derive the SOURCE method (Bae et al. 2024), which introduces a small number of approximations of Jacobians that reduce the full unrolling computation to something resembling an influence-function calculation at each of several checkpoints.

Stationary segments. Divide the $T$ training steps into $L$ segments at checkpoints $0=T_0<T_1<\cdots <T_L=T$, and let $K_{\ell }=T_{\ell }-T_{\ell -1}$ be the number of steps in segment $\ell$. Within each segment, make the stationarity approximation: the expected mini-batch Hessian and per-example gradient are approximately constant,

Similarly, write ${\bar{\eta }}_{\ell }$ for the average learning rate in segment $\ell$.

Segment Jacobian. Under the stationarity approximation, the product of $K_{\ell }$ step Jacobians within segment $\ell$ simplifies to a matrix power:

The matrix exponential is well approximated by the matrix power when ${\bar{\eta }}_{\ell }{\bar{H}}_{\ell }$ has small spectral norm. In an eigenbasis of ${\bar{H}}_{\ell }$ with eigenvalue $\sigma$, the segment Jacobian acts as the scalar filter $F_S(\sigma )=e^{-{\bar{\eta }}_{\ell }{K}_{\ell }\sigma }$: high-curvature directions ($\sigma$ large) are exponentially damped, while low-curvature directions ($\sigma$ small) are preserved.

Segment response. Return to the unrolling formula (20) and restrict the sum to the steps inside segment $\ell$. The contribution of datum $z_i$ during this segment is

where the Jacobian $J_{(t+1):T_{\ell }}$ propagates each gradient contribution forward to the segment boundary $T_{\ell }$ (propagation from $T_{\ell }$ to the end of training is handled by the inter-segment Jacobians ${\bar{S}}_{{\ell }^{\prime }}$ in the full formula below). Now apply the stationarity approximation: replace each ${\nabla }_{w}L(w_t,z_i)$ by ${\bar{g}}_{i,\ell }$, each learning rate by ${\bar{\eta }}_{\ell }$, replace the indicator $1[i\in B_t]$ by its expectation $b/n$ (datum $z_i$ appears in a random batch with probability $b/n$), and use the segment Jacobian approximation $J_{(t+1):T_{\ell }}\approx (I-{\bar{\eta }}_{\ell }{\bar{H}}_{\ell }{)}^{T_{\ell }-1-t}$. This gives

$$R_{i,\ell }\approx \frac{{\bar{\eta }}_{\ell }}{n}\sum _{k=0}^{K_{\ell }-1}(I-{\bar{\eta }}_{\ell }{\bar{H}}_{\ell }{)}^k{\bar{g}}_{i,\ell }.$$

The sum is a matrix geometric series. Using ${\sum }_{k=0}^{K-1}{M}^k=(I-M^K)(I-M{)}^{-1}$ with $M=I-{\bar{\eta }}_{\ell }{\bar{H}}_{\ell }$ and passing to the matrix exponential:

To understand this, work in an eigenbasis of ${\bar{H}}_{\ell }$. On an eigenvalue $\sigma$, the matrix $(I-e^{-{\bar{\eta }}_{\ell }{K}_{\ell }{\bar{H}}_{\ell }}){\bar{H}}_{\ell }^{-1}$ acts as the scalar

Note the similarity to the classical influence function formula, which applies $1/\sigma$ (the inverse Hessian eigenvalue), and with the damped IF (6), which applies $1/(\sigma +\lambda )$. Thus $F_r$ can be seen as a principled damping constant:

• Large $\sigma$ (directions the optimizer has fully converged along): $e^{-{\bar{\eta }}_{\ell }{K}_{\ell }\sigma }\approx 0$, so $F_r(\sigma )\approx 1/\sigma$. Agrees with the classical IF.
• Small $\sigma$ (directions the optimizer has barely moved along): $e^{-{\bar{\eta }}_{\ell }{K}_{\ell }\sigma }\approx 1-{\bar{\eta }}_{\ell }{K}_{\ell }\sigma$, so $F_r(\sigma )\approx {\bar{\eta }}_{\ell }{K}_{\ell }$. A finite constant, not $1/\sigma \to \infty$.
• Crossover at ${\sigma }_{\ast }\approx 1/({\bar{\eta }}_{\ell }{K}_{\ell })$: directions with curvature below this threshold are automatically suppressed.

Full formula. Chaining the segment Jacobians and responses across all $L$ segments:

where the product is ordered with ${\ell }^{\prime }=L$ on the left. With a single segment ($L=1$), the formula reduces to ${\nabla }_w{\varphi }^{\top }{\bar{r}}_{i,1}$, which is the classical IF with $H^{-1}$ replaced by the filter $F_r$.

4.2.2 Implicit JVPs and the REPLAY algorithm

The second approach to evaluating (20) avoids materializing Jacobians altogether. Instead, it uses reverse-mode automatic differentiation through the training loop, computing Jacobian-vector products (JVPs) implicitly. This is the strategy of the MAGIC method (Ilyas & Engstrom 2025), building on the metagradient framework of Engstrom et al. 2025.

Setup. Model the entire training process as a composition of $T$ differentiable update functions:

where $s_t$ is the full optimizer state at step $t$ (parameters, momentum buffers, etc.), $g_t(s_t,\beta )={\sum }_{i\in B_t}{\beta }_i\nabla L(s_t,z_i)$ is the weighted mini-batch gradient, and $h_t$ is the optimizer's update rule (SGD, Adam, etc.). The attribution target is $\varphi (s_T)$, which depends on $\beta$ through the entire chain of updates.

The metagradient decomposition. Applying the chain rule through this composition:

The vector $A_{t+1}:=\partial \varphi (s_T)/\partial s_{t+1}$ is the sensitivity of the final observable to the state at step $t+1$; it satisfies the backward recursion

This is just backpropagation through the training loop. The vector $B_t$ is the direct effect of $\beta$ at step $t$: it captures how the gradient at step $t$ depends on the data weights. For vanilla SGD, $B_t$ reduces to $-({\eta }_t/b)[{\nabla }_{w}L(w_t,z_i)\cdot 1[i\in B_t]{]}_{i=1}^n$.

REPLAY: efficient reverse-mode through training. The backward recursion (33) requires the optimizer state $s_t$ at every step. Naively this means storing all $T$ states, which is prohibitive for long training runs. The REPLAY algorithm (Engstrom et al. 2025) solves this with a hierarchical checkpointing scheme:

1. Save $k$ evenly-spaced checkpoints along the training trajectory.
2. When the backward pass needs a state between two checkpoints, replay the training forward from the nearest saved checkpoint to regenerate it.
3. Apply this strategy recursively (a $k$-ary tree of depth ${\log }_{k}T$).

The result is exact (to floating-point precision) differentiation through all $T$ training steps, using $O(k{\log }_{k}T)$ stored states and $O(T{\log }_{k}T)$ total forward-pass work, a logarithmic overhead over the cost of training itself.

Comparison with explicit unrolling. The implicit approach has two main advantages:

• Exactness. No stationarity or matrix-exponential approximations are needed; the Jacobian-vector products are computed exactly by autodiff.
• Generality. Any differentiable optimizer (Adam, LAMB, etc.) and any training pipeline (multi-stage, curriculum learning) are handled transparently, since $h_t$ is just a function.

The price is computational: REPLAY requires $O(T\log T)$ training-equivalent steps and careful engineering of the checkpointing schedule, whereas the explicit approach requires only a handful of checkpoints and an EK-FAC-style Hessian approximation at each. The choice depends on the scale: for moderate-size models where replaying training is feasible, MAGIC gives exact answers; for large-scale models where even one extra training pass is expensive, the SOURCE approximation with $L=3$–$6$ segments may be the only practical option.

4.3 Long Exercise: From unrolling to influence functions

The unrolling formula (20) and the influence function (6) appear to be very different objects: one is a sum over training steps of Jacobian-propagated gradients, the other is a single inverse-Hessian-times-gradient computation at the endpoint. The following exercise shows that, in the right limit, unrolling converges to the classical influence function. The key observation, due to Mlodozeniec et al. 2025, is that the joint process $(w_t,r_t)$, where $r_t:=\partial w_t/\partial \epsilon {|}_{\epsilon =0}$ is the unrolling response, forms a Markov chain, and its limiting behavior can be analyzed with standard stochastic-approximation tools. This yields a convergence result that requires only that SGD reaches a local minimum, not that the loss is convex.

Exercise 4.2 (Convergence of unrolling to influence functions). Consider the SGD update with an interpolated loss that downweights example $z_k$ by $\epsilon$:

$$w_{t+1}(\epsilon )=w_t(\epsilon )-\frac{{\eta }_t}{b}\sum _{i\in B_t}(1-\epsilon 1[i=k]){\nabla }_{w}L(w_t(\epsilon ),z_i),$$

so $\epsilon =0$ is the original training run and $\epsilon =1/n$ corresponds to removing $z_k$. Define the response $r_t:=\partial w_t(\epsilon )/\partial \epsilon {|}_{\epsilon =0}$.

$$\left(\begin{array}{c}{w}_{t+1}\\ r_{t+1}\end{array}\right)=\underset{\text{propagation}}{\underbrace{\left(\begin{array}{cc}I& 0\\ 0& J_t\end{array}\right)}}\left(\begin{array}{c}{w}_t\\ r_t\end{array}\right)+\underset{\text{driving terms}}{\underbrace{\left(\begin{array}{c}-\frac{{\eta }_t}{b}{\sum }_{i\in B_t}{\nabla }_{w}L(w_t,z_i)\\ \frac{{\eta }_t}{b}1[k\in B_t]{\nabla }_{w}L(w_t,z_k)\end{array}\right)}}$$

(a) (The response recursion.) By differentiating (34) with respect to $\epsilon$ at $\epsilon =0$, derive (35), where $J_t=I-\frac{{\eta }_t}{b}{\sum }_{i\in B_t}{\nabla }_w^{2}L(w_t,z_i)$ is the step Jacobian from (17). The top row is ordinary SGD; verify that the bottom row is the same recursion as Exercise 4.1(a), but now with stochastic batches. Since the driving terms depend only on $(w_t,r_t)$ and the i.i.d. batch selection $B_t$, the joint process is a Markov chain. Why is this observation useful?

(b) (Deterministic warm-up: full-batch GD.) As a warm-up, consider full-batch gradient descent with constant learning rate $\eta$ (i.e. $B_t=\{1,\dots ,n\}$ for all $t$). Assume training has converged: $w_t\approx w^{\ast }$ and ${\nabla }_w^{2}L(w_t)\approx H$ for all $t$ in a window of length $K$. Show that the response over this window satisfies

$$r_T\approx -[I-(I-\eta H{)}^K]H^{-1}{\nabla }_{w}L(w^{\ast },z_k).$$

Hint: Sum the geometric series ${\sum }_{t=0}^{K-1}(I-\eta H{)}^t$ using the identity ${\sum }_{t=0}^{K-1}{M}^t=(I-M^K)(I-M{)}^{-1}$.

Take $K\to \infty$ (assuming $\eta <2/{\lambda }_{\max }(H)$) and recover the influence function $r_{\infty }=-H^{-1}{\nabla }_{w}L(w^{\ast },z_k)=r_{\mathrm{IF}}$.

(c) (Stochastic case: convergence to IF.) Now return to SGD with i.i.d. batches and a decaying learning rate satisfying ${\sum }_t{\eta }_t=\infty$, ${\sum }_t{\eta }_t^2<\infty$ (the Robbins–Monro conditions). Assume SGD converges to a local minimum $w^{\ast }$ with $H={\nabla }_w^{2}L(w^{\ast })$ positive semidefinite. The continuous-time ODE that the response tracks is

$$\dot{r}(t)=-Hr(t)+{\nabla }_{w}L(w^{\ast },z_k).$$

(You do not need to prove that SGD tracks this ODE, this follows from standard stochastic approximation theory.)

i. Show that the equilibrium of (36) in the column space of $H$ is $r_{\mathrm{IF}}=H^{+}{\nabla }_{w}L(w^{\ast },z_k)$, where $H^{+}$ is the pseudoinverse. This is the influence function, with $H^{+}$ in place of $H^{-1}$ because the Hessian may be singular at a local minimum of an overparameterized model.

ii. Show that the component of $r(t)$ in the null space of $H$ grows linearly: if $P_0$ is the projector onto $\ker (H)$, then $P_{0}r(t)=P_{0}r(0)+t{P}_0{\nabla }_{w}L(w^{\ast },z_k)$. Why does this component not converge? Under what condition on ${\nabla }_{w}L(w^{\ast },z_k)$ does this runaway term vanish?

The full result (Mlodozeniec et al. 2025, Theorem 2) is: on the set of SGD trajectories that converge to a local minimum, $r_t\to r_{\mathrm{IF}}+r_{\mathrm{NS}}$ almost surely, where $r_{\mathrm{NS}}\in \ker (H)$. The influence function is the limiting response, up to a component in the flat directions of the loss.

Solution

(a) Differentiating (34) with respect to $\epsilon$ at $\epsilon =0$: the first row is just the SGD update, independent of $\epsilon$. For the response, the chain rule gives two contributions: the indirect effect through $w_t(\epsilon )$ picks up the mini-batch Hessian (giving $J_t{r}_t$ as in Exercise 4.1), and the direct effect from the $(1-\epsilon 1[i=k])$ factor contributes $+({\eta }_t/b)1[k\in B_t]{\nabla }_{w}L(w_t,z_k)$. This is (35). Given the current state $(w_t,r_t)$ and the batch selection $B_t$ (which is i.i.d. and independent of history), the next state $(w_{t+1},r_{t+1})$ depends only on $(w_t,r_t)$ — no earlier states are needed. This is the Markov property.

(b) With full-batch GD, $B_t=\{1,\dots ,n\}$, the response recursion becomes deterministic: $r_{t+1}=(I-\eta H)r_t+\eta {\nabla }_{w}L(w^{\ast },z_k)$. With $r_0=0$, unrolling gives $r_K=\eta {\sum }_{t=0}^{K-1}(I-\eta H{)}^t{\nabla }_{w}L(w^{\ast },z_k)$. Using ${\sum }_{t=0}^{K-1}{M}^t=(I-M^K)(I-M{)}^{-1}$ with $M=I-\eta H$:

$$r_K=[I-(I-\eta H{)}^K]H^{-1}{\nabla }_{w}L(w^{\ast },z_k).$$

(Note the sign: the interpolation $1-\epsilon 1[i=k]$ downweights $z_k$, so the direct effect has the opposite sign from the $\beta$-upweighting in Exercise 4.1.) Since $\Vert I-\eta H\Vert <1$ when $\eta <2/{\lambda }_{\max }(H)$, $(I-\eta H{)}^K\to 0$ and $r_{\infty }=H^{-1}{\nabla }_{w}L(w^{\ast },z_k)=r_{\mathrm{IF}}$.

(c) i. At equilibrium $\dot{r}=0$, the ODE (36) gives $H{r}_{\infty }={\nabla }_{w}L(w^{\ast },z_k)$. In the column space of $H$, this has the unique solution $r_{\mathrm{IF}}=H^{+}{\nabla }_{w}L(w^{\ast },z_k)$. (If $H$ is invertible, $H^{+}=H^{-1}$ and this is the classical IF.)

ii. Project the ODE onto $\ker (H)$: $P_0\dot{r}(t)=-P_{0}Hr(t)+P_0{\nabla }_{w}L(w^{\ast },z_k)=P_0{\nabla }_{w}L(w^{\ast },z_k)$, since $P_{0}H=0$. Integrating: $P_{0}r(t)=P_{0}r(0)+t{P}_0{\nabla }_{w}L(w^{\ast },z_k)$. This grows linearly unless $P_0{\nabla }_{w}L(w^{\ast },z_k)=0$, i.e. unless the per-example gradient has no component in the Hessian null space. The null space of $H$ at a local minimum corresponds to flat directions along the minimum manifold; the response diverges if the perturbation "pushes" along these flat directions, because the optimizer has no restoring force.

4.4 Exercise: Influence depends on training time

The previous sections analyzed unrolling's asymptotics in the converged limit. But the whole point of unrolling is that the training path matters. The following exercise makes this concrete in a setting where the influence function admits a closed form at every point during training: a two-layer deep linear network (DLN), the same class of models from the SLT lecture day. The setup and analytic formula follow Lee et al. 2025.

Exercise 4.3 (Stagewise influence in a deep linear network). A two-layer linear network $f_W(x)=W_2{W}_{1}x$ with $W_1,W_2\in {\mathbb{R}}^{d\times d}$ is trained on $\{(x_i,y_i){\}}_{i=1}^n$ with squared loss. Assume whitened inputs ($\frac{1}{n}{\sum }_i{x}_i{x}_i^{\top }=I$) and, for simplicity, that the input–output cross-covariance ${\Sigma }_{xy}=\frac{1}{n}{\sum }_i{x}_i{y}_i^{\top }$ is already diagonal with entries $s_1>s_2>\cdots >s_d>0$.⁵

Dynamics. Under gradient flow with small balanced initialization, the network's modes evolve independently. The effective weight matrix at time $t$ is $W(t)=\mathrm{diag}({\mathcal{G}}_1(t),\dots ,{\mathcal{G}}_d(t))$, where each mode strength follows

$${\mathcal{G}}_k(t)=\frac{s_k{e}^{2s_{k}t/\tau }}{e^{2s_{k}t/\tau }-1+s_k/{\mathcal{G}}_k(0)},$$

with ${\mathcal{G}}_k(0)\ll s_k$ and time constant $\tau$. Mode $k$ transitions from $\approx 0$ to $\approx s_k$ around time $t_k^{\ast }\approx \frac{\tau }{2s_k}\log (s_k/{\mathcal{G}}_k(0))$. Larger singular values saturate first: the network learns dominant structure before fine structure.

Analytic influence. Now perturb the weight of training example $z_p=(x_p,y_p)$ by $\epsilon$, so the perturbed cross-covariance is ${\Sigma }_{xy}(\epsilon )={\Sigma }_{xy}+\epsilon y_p{x}_p^{\top }$. Since ${\Sigma }_{xy}$ is diagonal, the perturbation shifts its eigenvalues and rotates its eigenvectors. The influence of $z_p$ on the weight matrix at time $t$ decomposes as

$${\mathcal{I}}_p(t):=\frac{\partial W(t,\epsilon )}{\partial \epsilon }{|}_{\epsilon =0}=A\mathcal{G}(t)+{\mathcal{G}}_{\epsilon }^{\prime }(t)+\mathcal{G}(t)B,$$

with three terms:

• $A$ and $B$ are skew-symmetric matrices describing how the perturbation rotates the left/right singular bases. Their off-diagonal entries are $A_{jk}=B_{kj}=(y_p{)}_j(x_p{)}_k/(s_k-s_j)$ for $j\ne k$.
• ${\mathcal{G}}_{\epsilon }^{\prime }(t)=\mathrm{diag}(g_t^{\prime }(s_1)(y_p{)}_1(x_p{)}_1,\dots ,g_t^{\prime }(s_d)(y_p{)}_d(x_p{)}_d)$ encodes the change in mode strengths, where

$$g_t^{\prime }(s_k):=\frac{\partial {\mathcal{G}}_k(t)}{\partial s_k}$$

is the sensitivity of the $k$-th mode strength to a change in the corresponding singular value at time $t$.

The influence on the loss of a test example $z_q=(x_q,y_q)$, with residual $r_q(t)=y_q-W(t)x_q$, is

$$\mathcal{I}(z_p,{\ell }_q)(t)=-r_q(t{)}^{\top }(A\mathcal{G}(t)+{\mathcal{G}}_{\epsilon }^{\prime }(t)+\mathcal{G}(t)B)x_q.$$

(a) (Three sources of influence.) Interpret the three terms in (38):

• $A\mathcal{G}(t)$ and $\mathcal{G}(t)B$: the perturbation rotates the singular bases, mixing already-learned modes into each other. Why are these terms proportional to the current mode strengths $\mathcal{G}(t)$ rather than to their derivatives?
• ${\mathcal{G}}_{\epsilon }^{\prime }(t)$: the perturbation shifts the singular values, changing how fast each mode is learned. Why does this term involve $g_t^{\prime }(s_k)$, sensitivity of the dynamics to the singular value, rather than ${\mathcal{G}}_k(t)$ itself?

(b) (When does influence peak?) Using the dynamics (37), argue that $g_t^{\prime }(s_k)$ is peaked around $t\approx t_k^{\ast }$ (the transition time of mode $k$). Conclude that the ${\mathcal{G}}_{\epsilon }^{\prime }$ term, the part of influence that acts through the learning speed of each mode, is concentrated in time around the moment the mode is being learned. Before and after, this contribution is negligible.

Hint: Consider the limits $t\ll t_k^{\ast }$ (mode not yet learning, ${\mathcal{G}}_k\approx {\mathcal{G}}_k(0)$) and $t\gg t_k^{\ast }$ (mode saturated, ${\mathcal{G}}_k\approx s_k$). In both cases, how sensitive is ${\mathcal{G}}_k$ to a small change in $s_k$?

(c) (Sign flips.) Consider $d=2$ with $s_1\gg s_2$: mode 1 captures a coarse distinction ("animal vs. plant") and mode 2 a fine distinction ("dog vs. cat" within animals). A dog example $z_{\mathrm{dog}}$ and a cat example $z_{\mathrm{cat}}$ share the same mode-1 coordinate but have opposite mode-2 coordinates: $(x_{\mathrm{dog}}{)}_2=+(x_{\mathrm{cat}}{)}_2$. Using (40), argue that the influence of $z_{\mathrm{dog}}$ on a cat test loss can change sign during training: positive while mode 1 is being learned (shared structure), negative after mode 2 is learned (competing structure). At what time is the sign flip sharpest?

Solution

(a) The $A\mathcal{G}(t)$ and $\mathcal{G}(t)B$ terms describe how the perturbation rotates the singular basis. A rotation mixes already-learned modes into each other, so its effect on the weight matrix is proportional to the current mode strengths: if mode $k$ has not yet been learned (${\mathcal{G}}_k\approx 0$), rotating its basis direction contributes nothing to $W(t)$. The ${\mathcal{G}}_{\epsilon }^{\prime }(t)$ term, by contrast, acts through the learning dynamics: perturbing $s_k$ changes the speed at which mode $k$ is learned, and this matters precisely when the dynamics are sensitive to $s_k$ — not when ${\mathcal{G}}_k(t)$ is large or small per se, but when ${\mathcal{G}}_k(t)$ is changing rapidly as a function of $s_k$. This is why the derivative $g_t^{\prime }(s_k)=\partial {\mathcal{G}}_k(t)/\partial s_k$ appears rather than ${\mathcal{G}}_k(t)$ itself.

(b) From (37), ${\mathcal{G}}_k(t)=s_k{e}^{2s_{k}t/\tau }/(e^{2s_{k}t/\tau }-1+s_k/{\mathcal{G}}_k(0))$. For $t\ll t_k^{\ast }$: $e^{2s_{k}t/\tau }\approx 1$, so ${\mathcal{G}}_k\approx {\mathcal{G}}_k(0)$ regardless of $s_k$ — the mode hasn't started learning yet, and a small change in $s_k$ makes no difference, so $g_t^{\prime }(s_k)\approx 0$. For $t\gg t_k^{\ast }$: $e^{2s_{k}t/\tau }\gg s_k/{\mathcal{G}}_k(0)$, so ${\mathcal{G}}_k\approx s_k$ and $g_t^{\prime }(s_k)\approx 1$ — but the residual in mode $k$ is also $\approx 0$, so this term's contribution to the loss influence (40) is negligible. The product $g_t^{\prime }(s_k)\times (\text{residual in mode }k)$ is peaked around $t\approx t_k^{\ast }$, when the mode is in transition: the dynamics are maximally sensitive to $s_k$ and the residual is still nonzero.

(c) Write $x_{\mathrm{dog}}=(\alpha ,+\delta )$ and $x_{\mathrm{cat}}=(\alpha ,-\delta )$, where $\alpha$ is the shared mode-1 component and $\pm \delta$ are opposite mode-2 components (since ${\Sigma }_{xy}$ is diagonal, the standard basis is the singular basis). During mode-1 learning ($t\approx t_1^{\ast }$, mode 2 not yet active): the residual of the cat test example has a large mode-1 component, and $z_{\mathrm{dog}}$'s perturbation increases $s_1$ (via $(y_{\mathrm{dog}}{)}_1(x_{\mathrm{dog}}{)}_1>0$), accelerating mode-1 learning. This reduces the cat test loss — positive influence. After mode-2 learning ($t\approx t_2^{\ast }$): the cat's residual is now dominated by mode 2, and $z_{\mathrm{dog}}$'s mode-2 perturbation has the wrong sign for the cat (because $(x_{\mathrm{dog}}{)}_2$ and $(x_{\mathrm{cat}}{)}_2$ have opposite signs). This increases the cat test loss — negative influence. The sign flip is sharpest around $t_2^{\ast }$, when $g_t^{\prime }(s_2)$ peaks and the mode-2 residual transitions from large to small.

5. Practical Considerations and Open Problems

Under construction.

6. Further Readings

Under construction.

References

Footnotes

1. In Section 4 we make this analogy explicit. ↩

2. For an invertible matrix $M$ and vectors $u,v$: $(M+u{v}^{\top }{)}^{-1}=M^{-1}-\frac{M^{-1}u{v}^{\top }{M}^{-1}}{1+v^{\top }{M}^{-1}u}$. ↩

3. Concretely you can take the eigendecomposition of $H$ and only invert the non-zero eigenvalues. ↩

4. For small $\epsilon$ the map $w\mapsto w+\epsilon r(w)$ is a diffeomorphism, so the inverse is well-defined. The exact formula requires evaluating $r$ at the pre-image; the expression above is correct to the order we need. ↩

5. The general case replaces the standard basis with the left/right singular vectors $U,V$ of ${\Sigma }_{xy}$ throughout; see Remark 4.4. ↩