Influence functions for large language models
Influence functions promise to say which training examples made a model behave the way it does. The formula is one line long and contains an inverse Hessian. This page derives that line, explains why the inverse is not optional, and then prices every step of using it on a model the size of Llama‑3‑8B.
Companion: Influence Functions, the Kaczmarz Way, the same formulas in two dimensions where everything can be drawn.
A model with parameters $\theta\in\mathbb R^P$ is trained on examples $z_1,\dots,z_N$ by minimizing an average loss, $$\theta^\star=\arg\min_\theta J(\theta),\qquad J(\theta)=\frac1N\sum_{i=1}^N L(z_i,\theta).$$ Pick one training example $z_m$ and give it a little extra weight $\varepsilon$: $$\theta^\star(\varepsilon)=\arg\min_\theta\; J(\theta)+\varepsilon\,L(z_m,\theta).$$ This curve through parameter space is the response function. At $\varepsilon=0$ it is the model we have; at $\varepsilon=-1/N$ it is, to first order, the model trained without $z_m$.
Pick also a measurement $f(\theta)$, a scalar you care about. For a language model the standard choice is the log‑probability of a completion given a prompt, $f(\theta)=\log p_\theta(\text{completion}\mid\text{prompt})$. The influence of $z_m$ on $f$ is the slope of the measurement along the response curve: $$\mathcal I_f(z_m)=\frac{d}{d\varepsilon}\,f\big(\theta^\star(\varepsilon)\big)\Big|_{\varepsilon=0}.$$ Positive influence means more of $z_m$ would push the completion's log‑probability up. Ranking the training set by this number is what "finding the examples responsible for a behavior" means in practice.
Nothing so far mentions a Hessian. It appears the moment we try to compute the slope without retraining.
Derivation
The minimizer $\theta^\star(\varepsilon)$ is defined implicitly by the stationarity condition, which holds for every $\varepsilon$: $$\nabla J\big(\theta^\star(\varepsilon)\big)+\varepsilon\,\nabla L\big(z_m,\theta^\star(\varepsilon)\big)=0 .$$ Differentiate both sides with respect to $\varepsilon$ and evaluate at $\varepsilon=0$. The chain rule turns $\nabla J$ into its Jacobian, the Hessian $H=\nabla^2 J(\theta^\star)$, times the velocity of the curve: $$H\,\frac{d\theta^\star}{d\varepsilon}\Big|_0+\nabla L(z_m,\theta^\star)=0 \quad\Longrightarrow\quad \frac{d\theta^\star}{d\varepsilon}\Big|_0=-\,H^{-1}\,\nabla L(z_m,\theta^\star).$$ One more chain rule gives the influence on the measurement: $$\boxed{\;\mathcal I_f(z_m)=-\,\nabla f(\theta^\star)^{\!\top}\,H^{-1}\,\nabla L(z_m,\theta^\star)\;}$$ That is the whole formula (Cook 1977 for statistics; Koh and Liang 2017 for deep networks). The inverse Hessian is the Jacobian of the implicit function theorem: it converts a change in the objective into the change in the minimizer.
Four ways to see what $H^{-1}$ is doing
$J$ is a potential energy well and $\theta^\star$ rests at its bottom. The rest of the training data acts like springs holding $\theta$ in place, and $H$ is their stiffness matrix. Upweighting $z_m$ applies a force $-\nabla L_m$. Hooke's law with a matrix says the equilibrium moves by force divided by stiffness, $H^{-1}(-\nabla L_m)$. A plain gradient dot product pretends every direction has stiffness one.
Write $H^{-1}=H^{-1/2}H^{-1/2}$. Then $\mathcal I_f=-\langle H^{-1/2}\nabla f,\;H^{-1/2}\nabla L_m\rangle$: the ordinary dot product after both gradients are rescaled by the curvature. Directions in which nearly every training example pushes (frequent tokens, the bulk of the gradient distribution) have huge curvature and get shrunk. What survives is the part of $\nabla L_m$ that is distinctive. Without the whitening, the ranking is dominated by generic, high‑gradient‑norm examples that look influential for everything.
An SGD step on $z_m$ moves $\theta$ by $-\eta\nabla L_m$, straight along the new gradient. The influence step moves by $-H^{-1}\nabla L_m$, which is where the equilibrium lands after all the other data has pushed back. In the 2‑D companion page these are two different projections onto the same line: orthogonal in the Euclidean metric versus orthogonal in the metric $H$.
The flattest directions of $H$ are the ones the data barely constrains, so a small push there moves the model a lot. Wang et al. (2025) bin the eigenvalues of the curvature and show that attribution quality keeps improving as lower and lower curvature bins are solved accurately. Methods that keep only a top subspace throw that signal away, and so does an iterative solver that stops early.
Three things stand between the boxed formula and a working system.
With $P$ parameters the Hessian has $P^2$ entries. For $P=8\times10^9$ that is $6.4\times10^{19}$ numbers, about 250 exabytes in single precision. Its eigenvalues also include zeros and negatives at any real checkpoint, so $H^{-1}$ does not exist as written. Practice replaces $H$ by the Gauss‑Newton Hessian $G$ (the Fisher information for log‑loss with labels sampled from the model), which is positive semidefinite, and adds damping: $$\mathcal I_f(z_m)\approx-\,\nabla f^{\top}(G+\lambda I)^{-1}\nabla L(z_m).$$ Bae et al. (2022) show this is not a hack but the exact answer to a slightly different question: it is the derivative of the proximal Bregman response function, which stays well defined for non‑convex, non‑converged models. Damping caps the amplification of any direction at $1/\lambda$. The toy above shows the effect: raise $\lambda$ past $\sigma_2$ and candidate B's advantage collapses. Choosing $\lambda$ is therefore choosing how far into the flat tail you trust the linearization. Grosse et al. (2023) used one tenth of the mean eigenvalue; Wang et al. (2025) derive it from the training run as $1/(\eta K)$, where $\eta$ is the learning rate and $K$ the number of steps, because $K$ steps at rate $\eta$ can only have converged the directions with curvature above that scale.
Kronecker‑factored approximate curvature (K‑FAC) treats layers as independent and, within a linear layer $y=Wx$, assumes the input activations $a$ and output backpropagated gradients $g$ are independent, so the layer's block factors as $G_\ell\approx A\otimes B$ with $A=\mathbb E[aa^\top]$ and $B=\mathbb E[gg^\top]$. Eigenvalue‑corrected K‑FAC (EK‑FAC) keeps the Kronecker eigenbasis but fits the eigenvalues $\Lambda$ directly from data, which fixes most of the bias in the spectrum. The inverse is then a rotate, an elementwise divide, and a rotate back: $$(G_\ell+\lambda I)^{-1}\operatorname{vec}(V)=\operatorname{vec}\!\Big(Q_B\big[(Q_B^{\top}VQ_A)\oslash(\Lambda+\lambda)\big]Q_A^{\top}\Big).$$ For a layer of shape $b\times a$ that costs about $4ab(a+b)$ multiply‑adds per vector, comparable to a forward pass through the layer. The factors are computed once per checkpoint and shared by every query.
Alternatively, never form $G$: solve $(G+\lambda I)\,v=\nabla f$ with an iterative method that only needs matrix‑vector products, each of which costs about one forward‑plus‑backward pass on a mini‑batch. Conjugate gradient does not tolerate stochastic products, so the field uses Neumann iterations (LiSSA). Unpreconditioned, the flattest directions shrink by only about $1-\lambda/\sigma_{\max}$ per step, so convergence takes on the order of $\sigma_{\max}/\lambda$ iterations, routinely $10^3$ to $10^6$. Stopping after $J$ steps at step size $\alpha$ is equivalent to raising the damping by $1/(\alpha J)$, silently, and exactly in the directions that carry signal.
ASTRA (Wang et al. 2025) combines the two: initialize at the EK‑FAC answer and precondition the Neumann iteration with the EK‑FAC inverse. One step reproduces EK‑FAC exactly; a few hundred more remove its structural bias. The fixed point is the true damped inverse‑Hessian‑vector product regardless of how crude the preconditioner is. Measured against real retraining, that bias turns out to matter: on a convolutional network the linear datamodeling score rose from about 0.25 to about 0.6.
Everything below is priced in one unit, the cost of computing a gradient for one training sequence of $T$ tokens through a model with $P$ parameters: about $2PT$ floating‑point operations forward and $4PT$ backward, so $$C_g\approx 6\,P\,T\ \text{FLOPs}.$$ Hardware is expressed in H100‑hours at a peak of 989 TFLOP/s in bf16 with 40% utilization, roughly $4\times10^{14}$ FLOP/s sustained. The worked example is a Llama‑3‑8B‑shaped model: 32 blocks, width 4096, MLP hidden width 14336, influence computed over the MLP weights as in Grosse et al. (2023). Every number in this section is recomputed from the calculator, so change the assumptions and the tables follow.
One forward and backward pass per query to get $\nabla f(\theta^\star)$. Negligible. The only real question is what $f$ should be: the log‑probability of the whole completion, or of one token, and whether to condition on the prompt. The choice changes which examples come back.
Run $S$ sequences with labels sampled from the model, accumulating $A=\mathbb E[aa^\top]$ and $B=\mathbb E[gg^\top]$ for every MLP matrix. Eigendecompose each factor ( FLOPs in total, negligible), then make a second pass to fit the corrected eigenvalues $\Lambda$. The multiply‑accumulates for the covariances cost about as much as the gradients themselves, hence the factor three.
Note the memory. A layer of shape $b\times a$ has $ab$ weights but its factors have $a^2+b^2$ entries. With $b=3.5a$ that is about 3.8 times the weights, so the EK‑FAC factors for the MLP layers are larger than the model. This is a fixed cost per checkpoint, shared by all queries.
With EK‑FAC alone the solve is a few dense matrix products per layer, about sequence‑gradient equivalents. This is why Grosse et al. (2023) could treat the inverse as essentially free and point at the scan as the bottleneck.
ASTRA pays for accuracy with iterations. Each one is a Gauss‑Newton‑vector product on a mini‑batch, roughly two gradient computations per sequence, plus one EK‑FAC solve. At the default settings that is times the EK‑FAC‑only cost. Still small next to the scan, but it becomes the dominant per‑query cost once the scan is replaced by a lookup (step 4). Memory is a handful of parameter‑sized vectors: the iterate, its momentum, the query gradient, and the preconditioned residual.
For every candidate sequence: a forward and backward pass to get its gradient, then a dot product with each query's inverse‑Hessian‑vector product. The dot products are of the total even with $Q$ queries sharing the pass. The gradient is the cost, and it is the same gradient the model computed during training.
That gives the cleanest statement of the problem. Scanning the whole pretraining corpus costs $6PD$ FLOPs, which is exactly the cost of pretraining: FLOPs, . Nobody does that per batch of queries, so the candidate set is shrunk first, by TF‑IDF prefiltering in the 2023 paper, or by scanning a random slice.
The obvious amortization is to store every training gradient once and answer each query with dot products. At for the candidate set alone, that is not storable. Random projection to $k$ dimensions (TRAK, LoGra) brings it to , and a query becomes a $k$‑dimensional lookup over $N$ rows. The catch is that a dense projection of a $P$‑vector costs $2Pk$, about gradient‑equivalents per example, so the projection has to exploit the low‑rank structure of per‑layer gradients, and every published projection so far pays in retrieval quality. Grosse's team's current work is sketches with recall guarantees for the top influential examples; if that works, step 3 becomes the whole per‑query bill.
The linear datamodeling score retrains the model on random half‑subsets of the data, several seeds each, and correlates the measured change in $f$ with the sum of predicted influences over the subset. At the scale above that is only affordable for a fine‑tuning run, and even then it costs as much as the scan. Grosse et al. (2023) instead validated against the proximal Bregman response function, the quantity the estimate approximates; Wang et al. (2025) used about 10,000 retrained small models and 1,000 GPT‑2 fine‑tunes. Every claim of the form "better curvature estimates give more sensible retrievals" rests on this step.
| Step | FLOPs | H100‑hours | Paid | Scales with |
|---|---|---|---|---|
| Query gradients | per query | Q | ||
| EK‑FAC factors | per checkpoint | S | ||
| iHVP, EK‑FAC only | per query | Σ ab(a+b) | ||
| iHVP, ASTRA | per query | J · B | ||
| Scan of N candidates | per batch of Q queries | N | ||
| Scan of the full corpus | = pretraining | D | ||
| Validation retrains | per method evaluation | M · Dft | ||
| Per query, amortized over the batch (ASTRA + scan/Q) |
With EK‑FAC the inverse‑Hessian‑vector product costs a few gradient computations. The training‑set scan costs one gradient computation per candidate, and a full scan equals pretraining.
Storing gradients to reuse across queries needs petabytes for a modest candidate set. Compression to a few thousand dimensions is the only route to per‑query lookups, and quality is the price.
ASTRA multiplies the per‑query inverse cost by a few thousand but leaves it far below the scan. Once retrieval is a lookup, that iteration count is the per‑query bill, which is why preconditioner quality matters.
Ground truth means retraining, hundreds of times. It is only possible at fine‑tuning scale, and everything known about which approximations matter comes from there.