Influence Functions for Scalable Data Attribution in Diffusion Models

## Abstract

Abstract Diffusion models have led to significant advancements in generative modelling.
Yet their widespread adoption poses challenges regarding data attribution and interpretability.
In this paper, we aim to help address such challenges in diffusion models by developing an influence function framework.
Influence function-based data attribution methods approximate how a model’s output would have changed if some training data were removed.
In supervised learning, this is usually used for predicting how the loss on a particular example would change.
For diffusion models, we focus on predicting the change in the probability of generating a particular example via several proxy measurements.
We show how to formulate influence functions for such quantities and how previously proposed methods can be interpreted as particular design choices in our framework.
To ensure scalability of the Hessian computations in influence functions, we systematically develop K-FAC approximations based on generalised Gauss-Newton matrices specifically tailored to diffusion models.
We recast previously proposed methods as specific design choices in our framework, and show that our recommended method outperforms previous data attribution approaches on common evaluations, such as the Linear Data-modelling Score (LDS) or retraining without top influences, without the need for method-specific hyperparameter tuning.

## 1 Introduction

Generative modelling for continuous data modalities — like images, video, and audio — has advanced rapidly propelled by improvements in diffusion-based approaches.
Many companies now offer easy access to AI-generated bespoke image content.
However, the use of these models for commercial purposes creates a need for understanding how the training data influences their outputs.
In cases where the model’s outputs are undesirable, it is useful to be able to identify, and possibly remove, the training data instances responsible for those outputs.
Furthermore, as copyrighted works often make up a significant part of the training corpora of these models (Schuhmann et al., ), concerns about the extent to which individual copyright owners’ works influence the generated samples arise.
Some already characterise what these companies offer as “copyright infringement as a service” (Saveri & Butterick, ), which has caused a flurry of high-profile lawsuits Saveri & Butterick (; ).
This motivates exploring tools for data attribution that might be able to quantify how each group of training data points influences the models’ outputs.
Influence functions (Koh & Liang, ; Bae et al., ) offer precisely such a tool.
By approximating the answer to the question, “If the model was trained with some of the data excluded, what would its output be?”, they can help finding data points most responsible for a low loss on an example, or a high probability of generating a particular example.
However, they have yet to be scalably adapted to the general diffusion modelling setting.

Influence functions work by locally approximating how the loss landscape would change if some of the training data points were down-weighted in the training loss (illustrated in Figure 5).
Consequently, this enables prediction for how the (local) optimum of the training loss would change, and how that change in the parameters would affect a measurement of interest (e.g., loss on a particular example).
By extrapolating this prediction, one can estimate what would happen if the data points were fully removed from the training set.
However, to locally approximate the shape of the loss landscape, influence functions require computing and inverting the Hessian of the training loss, which is computationally expensive.
One common approximation of the training loss’s Hessian is the generalised Gauss-Newton matrix (GGN, Schraudolph, ; Martens, ).
The GGN has not been clearly formulated for the diffusion modelling objective before and cannot be uniquely determined based on its general definition.
Moreover, to compute and store a GGN for large neural networks further approximations are necessary.
We propose using Kronecker-Factored Approximate Curvature (K-FAC, Heskes, ; Martens & Grosse, ) to approximate the GGN.
It is not commonly known how to apply it to neural network architectures used in diffusion models;
for example, Kwon et al. () resort to alternative Hessian approximation methods because “[K-FAC] might not be applicable to general deep neural network models as it highly depends on the model architecture”.
However, based on recent work, it is indeed clear that it can be applied to architectures used in diffusion models (Grosse & Martens, ; Eschenhagen et al., ), which typically combine linear layers, convolutions, and attention (Ho et al., ).

In this work, we describe a scalable approach to influence function-based approximations for data attribution in diffusion models, using a K-FAC approximation of GGNs as Hessian approximations.
We articulate a design space based on influence functions, unify previous methods for data attribution in diffusion models (Georgiev et al., ; Zheng et al., ) through our framework, and argue for the design choices that distinguish our method from previous ones.
One important design choice is the GGN used as the Hessian approximation.
We formulate different GGN matrices for the diffusion modelling objective and discuss their implicit assumptions.
We empirically ablate variations of the GGN and other design choices in our framework and show that our proposed method outperforms the existing data attribution methods for diffusion models as measured by common data attribution metrics like the Linear Data-modelling Score (Park et al., ) or retraining without top influences.
Finally, we also discuss interesting empirical observations that challenge our current understanding of influence functions in the context of diffusion models.

Figure: Figure 1: Most influential training data points as identified by K-FAC Influence Functions for samples generated by a denoising diffusion probabilistic model trained on CIFAR-10. The top influences are those whose omission from the training set is predicted to most increase the loss of the generated sample. Negative influences are those predicted to most decrease the loss, and the most neutral are those that should change the loss the least.
Refer to caption: x1.png

## 2 Background

This section introduces the general concepts of diffusion models, influence functions, and the GGN.

### 2.1 Diffusion Models

Diffusion models are a class of probabilistic generative models that fit a model $p_{\theta}(x)$ parameterised by parameters $\theta\in\mathbb{R}^{d_{\tt{param}}}$ to approximate a training data distribution $q(x)$, with the primary aim being to sample new data $x\sim p_{\theta}(\cdot)$ (Sohl-Dickstein et al., ; Ho et al., ; Turner et al., ).
This is usually done by augmenting the original data $x$ with $T$ fidelity levels as $x^{(0:T)}=[x^{(0)},\dots,x^{(T)}]$ with an augmentation distribution $q(x^{(0:T)})$ that satisfies the following criteria: 1) the highest fidelity $x^{(0)}$ equals the original training data $q(x^{(0)})=q(x)$, 2) the lowest fidelity $x^{(T)}$ has a distribution that is easy to sample from, and 3) predicting a lower fidelity level from the level directly above it is simple to model and learn.
To achieve the above goals, $q$ is typically taken to be a first-order Gaussian auto-regressive (diffusion) process: $q(x^{(t)}|x^{(0:t-1)})={\mathcal{N}}(x^{(t)}|\lambda_{t}x^{(t-1)},(1-\lambda_{
t})^{2}I)$, with hyperparameters $\lambda_{t}$ set so that the law of $x^{(T)}$ approximately matches a standard Gaussian distribution ${\mathcal{N}}(0,I)$.
In that case, the reverse conditionals $q(x^{(t-1)}|x^{(t:T)})=q(x^{(t-1)}|x^{(t)})$ are first-order Markov, and if the number of fidelity levels $T$ is high enough, they can be well approximated by a diagonal Gaussian, allowing them to be modelled with a parametric model with a simple likelihood function, hence satisfying (3) (Turner et al., ).
The marginals $q(x^{(t)}|x^{(0)})=\mathcal{N}\left(x^{(t)}|\left(\prod_{t^{\prime}=1}^{t}
\lambda_{t^{\prime}}\right)x^{(0)},\left(1-\prod_{t^{\prime}=1}^{t}\lambda_{t^
{\prime}}^{2}\right)I\right)$ also have a simple Gaussian form, allowing for the augmented samples to be sampled as:

$$ $\displaystyle x^{(t)}=\prod\nolimits_{t^{\prime}=1}^{t}\lambda_{t}x^{(0)}+\Big {(}1-\prod\nolimits_{t^{\prime}=1}^{t}\lambda_{t^{\prime}}^{2}\Big{)}^{ \nicefrac{{1}}{{2}}}\epsilon^{(t)},$ $\displaystyle\text{with }\epsilon^{(t)}\sim{\mathcal{N}}(0,I).$ (1) $$

Diffusion models are trained to approximate the reverse conditionals $p_{\theta}(x^{(t-1)}|x^{(t)})\approx q(x^{(t-1)}|x^{(t)})$ by maximising log-probabilities of samples $x^{(t-1)}$ conditioned on $x^{(t)}$, for all timesteps $t=1,\dots,T$.
We can note that $q(x^{(t-1)}|x^{(t)},x^{(0)})$ has a Gaussian distribution with mean given by:

$$ $\mu_{t-1|t,0}(x^{(t)},{\epsilon}^{(t)})=\frac{1}{\lambda_{t}}\Bigg{(}x^{(t)}- \frac{1-\lambda_{t}^{2}}{\big{(}1-\prod_{t^{\prime}=1}^{t}\lambda_{t^{\prime}} ^{2}\big{)}^{\nicefrac{{1}}{{2}}}}\epsilon^{(t)}\Bigg{)},\qquad\text{with } \epsilon^{(t)}\!\mathrel{\stackrel{{\scriptstyle\textnormal{def}}}{{=}}}\! \frac{\big{(}x^{(t)}-\prod_{t^{\prime}=1}^{t}\lambda_{t^{\prime}}x^{(0)}\big{) }}{(1-\prod_{t^{\prime}=1}^{t}\lambda_{t^{\prime}}^{2})^{\nicefrac{{1}}{{2}}}}$ $$

as in Equation 1.
In other words, the mean is a mixture of the sample $x^{(t)}$ and the noise $\epsilon^{(t)}$ that was applied to $x^{(0)}$ to produce it.
Hence, we can choose to analogously parameterise $p_{\theta}(x^{(t-1)}|x^{(t)})$ as $\mathcal{N}\big{(}x^{(t-1)}|\mu_{t-1|t,0}\big{(}x^{(t)},{\epsilon}^{t}_{\theta
}(x^{(t)})\big{)},\sigma_{t}^{2}I\big{)}$.
That way, the model $\epsilon^{(t)}_{\theta}(x^{(t)})$ simply predicts the noise $\epsilon^{(t)}$ that was added to the data to produce $x^{(t)}$. The variances $\sigma_{t}^{2}$ are usually chosen as hyperparameters (Ho et al., ).
With that parameterisation, the negative expected log-likelihood $\mathbb{E}_{q(x^{t-1},x^{(t)}|x^{(0)})}\left[-\log p(x^{(t-1)}|x^{(t)})\right]$, up to scale and shift independent of $\theta$ or $x^{(0)}$, can be written as (Ho et al., ; Turner et al., ):(^1^11Note that the two random variables $x^{(t)},\epsilon^{(t)}$ are deterministic functions of one-another.)

$$ $\ell_{t}(\theta,x^{(0)})=\mathbb{E}_{{\epsilon}^{(t)},x^{(t)}}\left[\left \lVert{\epsilon}^{(t)}-{\epsilon}_{\theta}^{t}\left(x^{(t)}\right)\right\rVert ^{2}\right]\qquad\begin{aligned} &{\epsilon}^{(t)}\sim{\mathcal{N}}(0,I)\\ &x^{(t)}=\prod\nolimits_{t^{\prime}=1}^{t}\lambda_{t}x^{(0)}+\Big{(}1-\prod \nolimits_{t^{\prime}=1}^{t}\lambda_{t^{\prime}}^{2}\Big{)}^{\nicefrac{{1}}{{2 }}}\epsilon^{(t)}\end{aligned}$ (2) $$

This leads to a training loss $\ell$ for the diffusion model ${\epsilon}_{\theta}^{t}(x^{(t)})$ that is a sum of per-diffusion timestep training losses:(^2^22Equivalently, a weighted sum of per-timestep negative log-likelihoods $-\log p_{\theta}(x^{(t-1)}|x^{(t)})$.)

$$ $\ell(\theta,x)=\mathbb{E}_{\tilde{t}}\left[\ell_{\tilde{t}}(\theta,x)\right] \qquad\tilde{t}\sim\operatorname{Uniform}([T]).$ $$

The parameters are then optimised to minimise the loss averaged over a training dataset ${\mathcal{D}}\leavevmode\resizebox{3.33333pt}{0.0pt}{=}\{x_{n}\}_{n=1}^{N}$:

$$ $\displaystyle\theta^{\star}({\mathcal{D}})=\operatorname*{arg\,min}_{\theta}{ \mathcal{L}}_{\mathcal{D}}(\theta)$ $\displaystyle{\mathcal{L}}_{\mathcal{D}}(\theta)\mathrel{\stackrel{{ \scriptstyle\textnormal{def}}}{{=}}}\frac{1}{N}\sum_{n=1}^{N}\ell(\theta,x_{n}).$ (3) $$

Other interpretations of the above procedure exist in the literature (Song & Ermon, ; Song et al., ; ; Kingma et al., ).

### 2.2 Influence Functions

The aim of influence functions is to answer questions of the sort “how would my model behave were it trained on the training dataset with some datapoints removed”.
To do so, they approximate the change in the optimal model parameters in Equation 3 when some training examples $(x_{j})_{j\in{\mathcal{I}}}$, ${\mathcal{I}}=\{i_{1},\dots,i_{M}\}\subseteq[N]$, are removed from the dataset ${\mathcal{D}}$. To arrive at a tractable approximation, it is useful to consider a continuous relaxation of this question: how would the optimum change were the training examples $(x_{j})_{j\in{\mathcal{I}}}$ down-weighted by $\varepsilon\in\mathbb{R}$ in the training loss:

$$ $r_{-{\mathcal{I}}}(\varepsilon)=\operatorname*{arg\,min}_{\theta}\frac{1}{N} \sum_{n=1}^{N}\ell(\theta,x_{n})-\varepsilon\sum_{j\in{\mathcal{I}}}\ell( \theta,x_{j})$ (4) $$

The function $r_{-{\mathcal{I}}}:\mathbb{R}\to\mathbb{R}^{d_{\tt{param}}}$ (well-defined if the optimum is unique) is the *response function*. Setting $\varepsilon$ to $\nicefrac{{1}}{{N}}$ recovers the minimum of the original objective in Equation 3 with examples $(x_{i_{1}},\dots,x_{i_{M}})$ removed.

Under suitable assumptions (see Appendix A), by the Implicit Function Theorem (Krantz & Parks, ), the response function is continuous and differentiable at $\varepsilon=0$. *Influence functions* can be defined as a linear approximation to the response function $r_{-{\mathcal{I}}}$ by a first-order Taylor expansion around $\varepsilon=0$:

$$ $\displaystyle r_{-{\mathcal{I}}}(\varepsilon)$ $\displaystyle=r_{-{\mathcal{I}}}(0)$ $\displaystyle+\frac{dr_{-{\mathcal{I}}}(\varepsilon^{\prime})}{d\varepsilon^{ \prime}}\Bigr{|}_{\varepsilon^{\prime}=0}\varepsilon$ $\displaystyle+o(\varepsilon)$ (5) $\displaystyle=\theta^{\star}({\mathcal{D}})$ $\displaystyle+\sum_{j\in{\mathcal{I}}}\left(\nabla_{\theta^{\star}}^{2}{ \mathcal{L}}_{\mathcal{D}}(\theta^{\star})\right)^{-1}\nabla_{\theta^{\star}} \ell(\theta^{\star}\!\!,x_{j})\varepsilon$ $\displaystyle+{o}(\varepsilon),$ $$

as $\varepsilon\to 0$. See Appendix A for a formal derivation and conditions. The optimal parameters with examples $(x_{i})_{i\in{\mathcal{I}}}$ removed can be approximated by setting $\varepsilon$ to $\nicefrac{{1}}{{N}}$ and dropping the $o(\varepsilon)$ terms.

Usually, we are not directly interested in the change in parameters in response to removing some data, but rather the change in some *measurement* function $m(\theta^{\star}(\mathcal{D}),x^{\prime})$ at a particular test input $x^{\prime}$ (e.g. per-example test loss). We can further make a first-order Taylor approximation to $m(\cdot,x^{\prime})$ at $\theta^{\star}(\mathcal{D})$ — $m(\theta,x^{\prime})=m(\theta^{\star},x^{\prime})+\nabla^{\mathchoice{\raise 1
.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle
\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt
\hbox{$\scriptscriptstyle\intercal$}}}_{\theta^{\star}}m(\theta^{\star},x^{
\prime})(\theta-\theta^{\star})+o\left(\|\theta-\theta^{\star}\|_{2}\right)$ — and combine it with Equation 5 to get a simple linear estimate of the change in the measurement function:

$$ $\displaystyle m(r_{-{\mathcal{I}}}(\varepsilon),x^{\prime})$ $\displaystyle=m(\theta^{\star}\!\!,x^{\prime})+\sum_{j\in{\mathcal{I}}}\nabla^ {\mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt \hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{ \raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}}}_{\theta^{\star}\!}m( \theta^{\star}\!,x^{\prime})\left(\nabla_{\theta^{\star}}^{2}{\mathcal{L}}_{ \mathcal{D}}(\theta^{\star})\right)^{-1}\nabla_{\theta^{\star}}\ell(\theta^{ \star}\!\!,x_{j})\varepsilon+o(\varepsilon).$ (6) $$

#### 2.2.1 Generalised Gauss-Newton matrix

Computing the influence function approximation in Equation 5 requires inverting the Hessian $\nabla_{\theta}^{2}{\mathcal{L}}_{\mathcal{D}}(\theta)\in\mathbb{R}^{d_{\tt{
param}}\times d_{\tt{param}}}$.
In the context of neural networks, the Hessian itself is generally computationally intractable and approximations are necessary.
A common Hessian approximation is the generalised Gauss-Newton matrix (GGN).
We will first introduce the GGN in an abstract setting of approximating the Hessian for a general training loss $\mathcal{L}(\theta)=\mathbb{E}_{z}\left[\rho(\theta,z)\right]$, to make it clear how different variants can be arrived at for diffusion models in the next section.

In general, if we have a function $\rho(\theta,z)$ of the form $h_{z}\circ f_{z}(\theta)$, with $h_{z}$ a convex function, the GGN for an expectation $\mathbb{E}_{z}[\rho(\theta,z)]$ is defined as

$$ $\displaystyle\mathrm{GGN}(\theta)={\color[rgb]{0.80078125,0.22265625,0.1640625 }\definecolor[named]{pgfstrokecolor}{rgb}{0.80078125,0.22265625,0.1640625} \mathbb{E}_{z}\left[{\color[rgb]{0.30859375,0.609375,0.5625}\definecolor[named ]{pgfstrokecolor}{rgb}{0.30859375,0.609375,0.5625}\nabla_{\theta}^{\intercal}f _{z}(\theta)}{\color[rgb]{0.171875,0.37890625,0.76171875}\definecolor[named]{ pgfstrokecolor}{rgb}{0.171875,0.37890625,0.76171875}\left(\nabla^{2}_{f_{z}( \theta)}h_{z}(f_{z}(\theta))\right)}{\color[rgb]{0.30859375,0.609375,0.5625} \definecolor[named]{pgfstrokecolor}{rgb}{0.30859375,0.609375,0.5625}\nabla_{ \theta}f_{z}(\theta)}\right]},$ $$

where $\nabla_{\theta}f_{z}(\theta)$ is the Jacobian of $f_{z}$.
Whenever $f_{z}$ is (locally) linear, the $\mathrm{GGN}$ is equal to the Hessian $\mathbb{E}_{z}[\nabla_{\theta}^{2}\rho(\theta,z)]$.
Therefore, we can consider the GGN as an approximation to the Hessian in which we “linearise” the function $f_{z}$.
Note that any decomposition of $\rho(\theta,z)$ results in a valid GGN as long as $h_{z}$ is convex (Martens, ).(^3^33 $h_{z}$ is typically required to be convex to guarantee the resulting GGN is a positive semi-definite (PSD) matrix. A valid non-PSD approximation to the Hessian can be formed with a non-convex $h_{z}$ as well; all the arguments about the exactness of the GGN approximation for a linear $f_{z}$ would still apply. However, the PSD property helps with numerical stability of the matrix inversion, and guarantees that the GGN will be invertible if a small damping term is added to the diagonal.)We give two examples below.
(^bordercolor=orange!50^bordercolor=orange!50todo: bordercolor=orange!50BM: This is a monstrous footnote. Maybe offload some to the appendix/delete unnecessary detail.)

Option 1.
A typical choice would be for $f_{z}$ to be the neural network function on a training datapoint $z$, and for $h_{z}$ to be the loss function (e.g. $\ell_{2}$-loss), with the expectation $\mathbb{E}_{z}$ being taken over the empirical (training) data distribution;
we call the GGN for this split $\mathrm{GGN}^{\texttt{model}}$.
The GGN with this split is exact for linear neural networks (or when the model has zero residuals on the training data) (Martens, ).

$$ $\displaystyle\begin{aligned} &f_{z}\coloneq\text{mapping from parameters to model output}\\ &h_{z}\coloneq\text{loss function (e.g. $\ell_{2}$-loss)}\end{aligned}$ $\displaystyle\rightarrow\mathrm{GGN}^{\texttt{model}}(\theta)$ (7) $$

Option 2. Alternatively, a different GGN can be defined by using a trivial split of the loss $\rho(\theta,z)$ into the identity map $h_{z}\coloneq\operatorname{id}$ and the loss $f_{z}\coloneq\rho(\cdot,z)$, and again taking the expectation over the empirical data distribution.
With this split, the resulting GGN is

$$ $\displaystyle\begin{aligned} &f_{z}\coloneq\rho(\cdot,z)\\ &h_{z}\coloneq\operatorname{id}\end{aligned}$ $\displaystyle\rightarrow\mathrm{GGN}^{\texttt{loss}}(\theta)={\color[rgb]{ 0.80078125,0.22265625,0.1640625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.80078125,0.22265625,0.1640625}\mathbb{E}_{z}\left[{\color[rgb]{ 0.30859375,0.609375,0.5625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.30859375,0.609375,0.5625}\nabla_{\theta}\rho(\theta,z)}{\color[rgb]{ 0.30859375,0.609375,0.5625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.30859375,0.609375,0.5625}\nabla^{\mathchoice{\raise 1.3611pt\hbox{$ \displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0 .95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$ \scriptscriptstyle\intercal$}}}_{\theta}\rho(\theta,z)}\right]}.$ (8) $$

This is also called the empirical Fisher (Kunstner et al., ).
Note that $\mathrm{GGN}^{\texttt{loss}}$ is only equal to the Hessian under the arguably more stringent condition that $\rho(\cdot,z)$ — the composition of the model and the loss function — is linear.
This is in contrast to $\mathrm{GGN}^{\texttt{model}}$, for which only the mapping from the parameters to the model output needs to be (locally) linear.
Hence, we might prefer to use $\mathrm{GGN}^{\texttt{model}}$ for Hessian approximation whenever we have a nonlinear loss, which is the case for diffusion models.
(^bordercolor=green!50^bordercolor=green!50todo: bordercolor=green!50RT: I wonder here whether you could emphasise that this might be a reason to prefer the previous GGN over this one to foreshadow what you are going to say later.)
(^bordercolor=orange!50^bordercolor=orange!50todo: bordercolor=orange!50BM: I tried emphasising it a bit more. Do you think it could be emphasised even better?)

## 3 Scalable influence functions for diffusion models

In this section, we discuss how we adapt influence functions to the diffusion modelling setting in a scalable manner.
We also recast data attribution methods for diffusion models proposed in prior work (Georgiev et al., ; Zheng et al., ) as the result of particular design decisions in our framework, and argue for our own choices that distinguish our method from the previous ones.

### 3.1 Approximating the Hessian

In diffusion models, we want to compute the Hessian of the loss of the form

$$ ${\mathcal{L}}_{\mathcal{D}}(\theta)=\mathbb{E}_{x_{n}}\left[\ell(\theta,x_{n}) \right]=\mathbb{E}_{x_{n}}\left[\mathbb{E}_{\tilde{t}}\left[\mathbb{E}_{x^{( \tilde{t})},{\epsilon}^{(\tilde{t})}}\left[\|{\epsilon}^{(\tilde{t})}-{ \epsilon}_{\theta}^{\tilde{t}}(x^{(\tilde{t})})\|^{2}\right]\right]\right],$ $$

where $\mathbb{E}_{x_{n}}\left[\cdot\right]=\big{(}\frac{1}{N}\sum_{n=1}^{N}\cdot\big
{)}$ is the expectation over the empirical data distribution. (^4^44Generally, $\mathbb{E}_{x_{n}}$ might also subsume the expectation over data augmentations applied to the training data points (see Section G.8 for details on how this is handled).)We will describe how to formulate different GGN approximations for this setting.

#### 3.1.1 GGN for diffusion models

Option 1. To arrive at a GGN approximation, as discussed in Section 2.2.1, we can partition the function $\theta\mapsto\|{\epsilon}^{(t)}-{\epsilon}_{\theta}^{t}(x^{(t)})\|^{2}$ into the model output $\theta\mapsto{\epsilon}_{\theta}^{t}(x^{(t)})$ and the $\ell_{2}$-loss function $\|{\epsilon}^{(t)}-\cdot\|^{2}$.
This results in the GGN:

$$ $\displaystyle\begin{aligned} &f_{z}\coloneq{\epsilon}_{\theta}^{\tilde{t}}(x^{ (\tilde{t})})\\ &h_{z}\coloneq\lVert{\epsilon}^{(\tilde{t})}-\cdot\rVert^{2}\end{aligned}$ $\displaystyle\rightarrow\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}(\theta)$ $\displaystyle={\color[rgb]{0.80078125,0.22265625,0.1640625}\definecolor[named] {pgfstrokecolor}{rgb}{0.80078125,0.22265625,0.1640625}\mathbb{E}_{x_{n}}\left[ \mathbb{E}_{\tilde{t}}\left[\mathbb{E}_{x^{(\tilde{t})},{\epsilon}^{(\tilde{t} )}}\left[{\color[rgb]{0.30859375,0.609375,0.5625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.30859375,0.609375,0.5625}\nabla_{\theta}^{\mathchoice{ \raise 1.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$ \textstyle\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0 .68054pt\hbox{$\scriptscriptstyle\intercal$}}}{\epsilon}_{\theta}^{\tilde{t}} \left(x^{(\tilde{t})}\right)}{\color[rgb]{0.171875,0.37890625,0.76171875} \definecolor[named]{pgfstrokecolor}{rgb}{0.171875,0.37890625,0.76171875}\left( 2I\right)}{\color[rgb]{0.30859375,0.609375,0.5625}\definecolor[named]{ pgfstrokecolor}{rgb}{0.30859375,0.609375,0.5625}\nabla_{\theta}{\epsilon}_{ \theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\right]\right]\right]},$ (9) $$

where $I$ is the identity matrix. This correspond to “linearising” the neural network ${\epsilon}^{t}_{\theta}$.
For diffusion models, the dimensionality of the output of ${\epsilon}_{\theta}^{\tilde{t}}$ is typically very large (e.g. $32\!\times\!32\!\times\!3$ for CIFAR), so computing the Jacobians $\nabla_{\theta}{\epsilon}^{t}_{\theta}$ explicitly is still intractable.
However, we can express $\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}$ as

$$ $\displaystyle\mathrm{F}_{\mathcal{D}}(\theta)=\mathbb{E}_{x_{n}}\left[\mathbb{ E}_{\tilde{t}}\left[\mathbb{E}_{x_{n}^{(\tilde{t})}}\left[\mathbb{E}_{{ \epsilon}_{\texttt{mod}}}\left[g_{n}(\theta)g_{n}(\theta)^{\mathchoice{\raise 1 .3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle \intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt \hbox{$\scriptscriptstyle\intercal$}}}\right]\right]\right]\right],$ $\displaystyle{\epsilon}_{\texttt{mod}}\sim{\mathcal{N}}\left({\epsilon}_{ \theta}^{\tilde{t}}\left(x_{n}^{(\tilde{t})}\right),I\right)$ (10) $$

where $g_{n}(\theta)=\nabla_{\theta}{\lVert{\epsilon}_{\texttt{mod}}-{\epsilon}_{
\theta}^{\tilde{t}}(x_{n}^{(\tilde{t})})\rVert^{2}}\in\mathbb{R}^{d_{\texttt{
param}}}$;
see Appendix B for the derivation.
This formulation lends itself to a Monte Carlo approximation, since we can now compute gradients using auxiliary targets ${\epsilon}_{\texttt{mod}}$ sampled from the model’s output distribution, as shown in Equation 10.
$\mathrm{F}_{\mathcal{D}}$ can be interpreted as a kind of Fisher information matrix (Amari, ; Martens, ), but it is not the Fisher for the marginal model distribution $p_{\theta}(x)$.

Option 2. Analogously to Equation 8, we can also consider the trivial decomposition of $\ell(\cdot,x)$ into the identity map and the loss,
effectively “linearising” $\ell(\cdot,x)$. The resulting GGN is:

$$ $\displaystyle\begin{aligned} &f_{z}\coloneq\ell(\cdot,x_{n})\\ &h_{z}\coloneq\operatorname{id}\end{aligned}$ $\displaystyle\rightarrow\mathrm{GGN}_{\mathcal{D}}^{\texttt{loss}}(\theta)={ \color[rgb]{0.80078125,0.22265625,0.1640625}\definecolor[named]{pgfstrokecolor }{rgb}{0.80078125,0.22265625,0.1640625}\mathbb{E}_{x_{n}}\left[{\color[rgb]{ 0.30859375,0.609375,0.5625}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.30859375,0.609375,0.5625}\nabla_{\theta}\ell(\theta,x_{n})\nabla^{\intercal} _{\theta}\ell(\theta,x_{n})}\right]},$ (11) $$

where $\ell(\theta,x)$ is the diffusion training loss defined in Equation 2.
This Hessian approximation $\mathrm{GGN}_{\mathcal{D}}^{\texttt{loss}}$ turns out to be equivalent to the ones considered in the previous works on data attribution for diffusion models (Georgiev et al., ; Zheng et al., ; Kwon et al., ).
In contrast, in this work, we opt for $\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}$ in Equation 9, or equivalently $\mathrm{F}_{\mathcal{D}}$, since it is arguably a better-motivated approximation of the Hessian than $\mathrm{GGN}_{\mathcal{D}}^{\texttt{loss}}$ (c.f. Section 2.2.1).
(^bordercolor=orange!50^bordercolor=orange!50todo: bordercolor=orange!50BM: after deadline: see if we need to introduce the other “empirical” Fisher variant)

In Zheng et al. (), the authors explored substituting different (theoretically incorrect) training loss functions into the influence function approximation.
In particular, they found that replacing the loss $\lVert{\epsilon}^{(t)}-{\epsilon}^{t}_{\theta}(x^{(t)})\rVert^{2}$ with the square norm loss $\lVert{\epsilon}^{t}_{\theta}(x^{(t)})\rVert^{2}$ (effectively replacing the “targets” ${\epsilon}^{(t)}$ with 0 0) gave the best results.
Note that the targets $\epsilon^{(t)}$ do not appear in the expression for $\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}$ in Equation 9.(^5^55This is because the Hessian of an $\ell_{2}$-loss w.r.t. the model output is a multiple of the identity matrix.)
Hence, in our method substituting different targets would not affect the Hessian approximation.
In Zheng et al. (), replacing the targets only makes a difference to the Hessian approximation because they use $\mathrm{GGN}_{\mathcal{D}}^{\texttt{loss}}$ (an empirical Fisher) to approximate the Hessian.

#### 3.1.2 K-FAC for diffusion models

While $\mathrm{F}_{\mathcal{D}}(\theta)$ and $\mathrm{GGN}_{\mathcal{D}}^{\texttt{loss}}$ do not require computing full Jacobians or the Hessian of the neural network model, they involve taking outer products of gradients of size $\mathbb{R}^{d_{\texttt{param}}}$, which is still intractable.
Kronecker-Factored Approximate Curvature (Heskes, ; Martens & Grosse, , K-FAC) is a common scalable approximation of the GGN to overcome this problem.
It approximates the GGN with a block-diagonal matrix, where each block corresponds to one neural network layer and consists of a Kronecker product of two matrices.
Due to convenient properties of the Kronecker product, this makes the inversion and multiplication with vectors needed in Equation 6 efficient enough to scale to large networks.
K-FAC is defined for linear layers, including linear layers with weight sharing like convolutions (Grosse & Martens, ).
This covers most layer types in the architectures typically used for diffusion models.
When weight sharing is used, there are two variants – K-FAC-expand and K-FAC-reduce (Eschenhagen et al., ).
For our recommended method, we choose to approximate the Hessian with a K-FAC approximation of $\mathrm{F}_{\mathcal{D}}$, akin to Grosse et al. ().

For the parameters $\theta_{l}$ of layer $l$, the GGN $\mathrm{F}_{\mathcal{D}}$ in Equation 10 is approximated by

$$ $\displaystyle\mathrm{F}_{\mathcal{D}}(\theta_{l})\approx\frac{1}{N^{2}}\sum_{n =1}^{N}\mathbb{E}_{\tilde{t}}\left[\mathbb{E}_{x_{n}^{(\tilde{t})},{\epsilon}^ {(\tilde{t})}}\left[a_{n}^{(l)}a_{n}^{(l){\mathchoice{\raise 1.3611pt\hbox{$ \displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0 .95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$ \scriptscriptstyle\intercal$}}}}\right]\right]\otimes\sum_{n=1}^{N}\mathbb{E}_ {\tilde{t}}\left[\mathbb{E}_{x_{n}^{(\tilde{t})},{\epsilon}^{(\tilde{t})},{ \epsilon}_{\texttt{mod}}^{(\tilde{t})}}\left[b_{n}^{(l)}b_{n}^{(l){\mathchoice {\raise 1.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$ \textstyle\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0 .68054pt\hbox{$\scriptscriptstyle\intercal$}}}}\right]\right],$ (12) $$

with $a_{n}^{(l)}\in\mathbb{R}^{d_{\texttt{in}}^{l}}$ being the inputs to the $l$th layer for data point $x_{n}^{(\tilde{t})}$ and $b_{n}^{(l)}\in\mathbb{R}^{d_{\texttt{out}}^{l}}$ being the gradient of the $\ell_{2}$-loss w.r.t. the output of the $l$th layer, and $\otimes$ denoting the Kronecker product.(^6^66For the sake of a simpler presentation this does not take potential weight sharing into account.)
The approximation trivially becomes an equality for a single data point and also for deep linear networks with $\ell_{2}$-loss (Bernacchia et al., ; Eschenhagen et al., ).
We approximate the expectations in Equation 12 with Monte Carlo samples and use K-FAC-expand whenever weight sharing is used since the problem formulation of diffusion models corresponds to the expand setting in Eschenhagen et al. (); in the case of convolutional layers this corresponds to Grosse & Martens ().
Lastly, to ensure the Hessian approximation is well-conditioned and invertible, we follow standard practice and add a damping term consisting of a small scalar damping factor times the identity matrix.
We ablate these design choices in Section 4 (Figures 4, 7 and 9).

### 3.2 Gradient compression and query batching

In practice, we recommend computing influence function estimates in Equation 6 by first computing and storing the approximate Hessian inverse, and then iteratively computing the preconditioned inner products $\nabla^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3
611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle
\intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}}}_{\theta^{
\star}}m(\theta^{\star},x)\left(\nabla_{\theta^{\star}}^{2}{\mathcal{L}}_{
\mathcal{D}}(\theta^{\star})\right)^{-1}\nabla_{\theta^{\star}}\ell(\theta^{
\star},x_{j})$ for different training datapoints $x_{j}$.
Following Grosse et al. (), we use query batching to avoid recomputing the gradients $\nabla_{\theta^{\star}}\ell(\theta^{\star},x_{j})$ when attributing multiple samples $x$.
We also use gradient compression; we found that compression by quantisation works much better for diffusion models compared to the SVD-based compression used by Grosse et al. () (see Appendix C), likely due to the fact that gradients $\nabla_{\theta}\ell(\theta,x_{n})$ are not low-rank in this setting.

### 3.3 What to measure

For diffusion models, arguably the most natural question to ask might be, for a given sample $x$ generated from the model, how did the training samples influence the probability of generating a sample $x$?
For example, in the context of copyright infringement, we might want to ask if removing certain copyrighted works would substantially reduce the probability of generating $x$.
With influence functions, these questions could be interpreted as setting the measurement function $m(\theta,x)$ to be the (marginal) log-probability of generating $x$ from the diffusion model: $\log p_{\theta}(x)$.

Computing the marginal log-probability introduces some challenges.
Diffusion models have originally been designed with the goal of tractable sampling, and not log-likelihood evaluation.
Ho et al. (); Sohl-Dickstein et al. () only introduce a lower-bound on the marginal log-probability.
Song et al. () show that exact log-likelihood evaluation is possible, but it only makes sense in settings where the training data distribution has a density (e.g. uniformly dequantised data), and it only corresponds to the marginal log-likelihood of the model when sampling deterministically (Song et al., ).(^7^77Unless the trained model satisfies very specific “consistency” constraints (Song et al., , Theorem 2).)Also, taking gradients of that measurement, as required for influence functions, is non-trivial.
Hence, in most cases, we might need a proxy measurement for the marginal probability.
We consider a couple of proxies in this work:

- 1.
Loss. Approximate $\log p_{\theta}(x)$ with the diffusion loss $\ell(\theta,x)$ in Equation 2 on that particular example. This corresponds to the ELBO with reweighted per-timestep loss terms (see Figure 18).
- 2.
Probability of sampling trajectory. If the entire sampling trajectory $x^{(0:T)}$ that generated sample $x$ is available, consider the probability of that trajectory $p_{\theta}(x^{(0:T)})=p(x^{T})\prod_{t=1}^{T}p_{\theta}(x^{(t-1)}|x^{(t)})$.
- 3.
ELBO. Approximate $\log p_{\theta}(x)$ with an Evidence Lower-Bound (Ho et al., , eq. (5)).

Initially, we might expect ELBO to be the best motivated proxy, as it is the only one with a clear link to the marginal log-probability.
Probability of sampling trajectory might also appear sensible, but it doesn’t take into account the fact that there are multiple trajectories $x^{(0:T)}$ that all lead to the same final sample $x^{(0)}$, and it has the disadvantage of not being reparameterisation invariant.(^8^88We can rescale the latent variables $x^{(t)}$ without affecting the marginal distribution $p_{\theta}(x)$, but changing the probability density of any particular trajectory.)We empirically investigate these different proxies in Section 4.

## 4 Experiments

Evaluating Data Attribution. To evaluate the proposed data attribution methods, we primarily focus on two metrics: Linear Data Modelling Score (LDS) and retraining without top influences.
LDS measures how well a given attribution method can predict the relative magnitude in the change in a measurement as the model is retrained on (random) subsets of the training data.
For an attribution method $a({\mathcal{D}},{\mathcal{D}}^{\prime},x)$ that approximates how a measurement $m(\theta^{\star}({\mathcal{D}}),x)$ would change if a model was trained on an altered dataset ${\mathcal{D}}^{\prime}$, LDS measures the Spearman rank correlation between the predicted change in output and actual change in output after retraining on different subsampled datasets:

$$ $\displaystyle\operatorname{spearman}\left[\left(a({\mathcal{D}},\tilde{{ \mathcal{D}}}_{i},x)\right)_{i=1}^{M};\left(m(\theta^{\star}(\tilde{{\mathcal{ D}}}_{i}),x)\right)_{i=1}^{M}\right],$ $$

where $\tilde{{\mathcal{D}}}_{i}$ are independently subsampled versions of the original dataset ${\mathcal{D}}$, each containing $50\%$ of the points sampled without replacement. However, a reality of deep learning is that, depending on the random seed used for initialisation and setting the order in which the data is presented in training, training on a fixed dataset can produce different models with functionally different behaviour.
Hence, for any given dataset ${\mathcal{D}}^{\prime}$, different measurements could be obtained depending on the random seed used.
To mitigate the issue, Park et al. () propose to use an ensemble average measurement after retraining as the “oracle” target:

$$ $\displaystyle\operatorname{LDS}=\operatorname{spearman}\left[\Big{(}a({ \mathcal{D}},\tilde{{\mathcal{D}}}_{i},x)\Big{)}_{i=1}^{M};\Big{(}\frac{1}{K} \sum\nolimits_{k=1}^{K}m(\tilde{\theta}_{k}^{\star}(\tilde{{\mathcal{D}}}_{i}) ,x)\Big{)}_{i=1}^{M}\right],$ (13) $$

where $\tilde{\theta}_{k}^{\star}({\mathcal{D}}^{\prime})\in\mathbb{R}^{d_{\texttt{
param}}}$ are the parameters resulting from training on ${\mathcal{D}}^{\prime}$ with a particular seed $k$.

Retraining without top influences, on the other hand, evaluates the ability of the data attribution method to surface the most influential data points – namely, those that would most negatively affect the measurement $m(\theta^{\star}\!({\mathcal{D}}^{\prime}),x)$ under retraining from scratch on a dataset ${\mathcal{D}}^{\prime}$ with these data points removed.
For each method, we remove a fixed percentage of the most influential datapoints from ${\mathcal{D}}$ to create the new dataset ${\mathcal{D}}^{\prime}$, and report the change in the measurement $m(\theta^{\star}\!({\mathcal{D}}^{\prime}),x)$ relative to $m(\theta^{\star}\!({\mathcal{D}}),x)$ (measurement by the model trained on the full dataset ${\mathcal{D}}$).

In all experiments, we look at measurements on samples generated by the model trained on ${\mathcal{D}}$.
We primarily focus on Denoising Diffusion Probabilistic Models (DDPM) (Ho et al., ) throughout.
(^8^8footnotetext: Better LDS results can sometimes be obtained when looking at validation examples (Zheng et al., ), but diffusion models are used primarily for sampling, so attributing generated samples is of primary practical interest.)

Baselines We compare influence functions with K-FAC and $\mathrm{GGN}^{\texttt{model}}_{\mathcal{D}}$ (MC-Fisher; Equation 10) as the Hessian approximation (K-FAC Influence) to TRAK as formulated for diffusion models in Georgiev et al. (); Zheng et al. ().
In our framework, their method can be tersely described as using $\mathrm{GGN}^{\texttt{loss}}_{\mathcal{D}}$ (Empirical Fisher) in Equation 11 as a Hessian approximation instead of $\mathrm{GGN}^{\texttt{model}}_{\mathcal{D}}$ (MC-Fisher) in Equation 10, and computing the Hessian-preconditioned inner products using random projections (Dasgupta & Gupta, ) rather than K-FAC.
We also compare to the ad-hoc changes to the measurement/training loss in the influence function approximation (D-TRAK) that were shown by Zheng et al. () to give improved performance on LDS benchmarks.
Note that, the changes in D-TRAK were directly optimised for improvements in LDS scores in the diffusion modelling setting, and lack any theoretical motivation.
Hence, a direct comparison for the changes proposed in this work (K-FAC Influence) is TRAK;
the insights from D-TRAK are orthogonal to our work.
(^bordercolor=orange!50^bordercolor=orange!50todo: bordercolor=orange!50BM: given the workshop review I really want to pre-empt confusion here)
These are the only prior works motivated by predicting the change in a model’s measurements after retraining that have been applied to the general diffusion modelling setting that we are aware of.
We also compare to naïvely using cosine similarity between the CLIP (Radford et al., ) embeddings of the training datapoints and the generated sample as a proxy for influence on the generated samples. Lastly, we report LDS results for the oracle method of “Exact Retraining”, where we actually retraining a single model to predict the changes in measurements.

Figure: (a) LDS results on the loss measurement.
Refer to caption: x2.png

LDS. The LDS results attributing the loss and ELBO measurements are shown in Figures 2(a) and 2(b).
K-FAC Influence outperforms TRAK in all settings.
K-FAC Influence using the loss measurement also outperforms the benchmark-tuned changes in D-TRAK in all settings as well.
In Figures 2(a) and 2(b), we report the results for both the best damping values from a sweep (see Appendix D), as well as for “default” values following recommendations in previous work (see Section G.4).
TRAK and D-TRAK appear to be significantly more sensitive to tuning the damping factor than K-FAC Influence. They often don’t perform at all if the damping factor is too small, and take a noticeable performance hit if the damping factor is not tuned to the problem or method (see Figures 8 and 10 in Appendix D).
However, in most applications, tuning the damping factor would be infeasible, as it requires retraining the model many times over to construct an LDS benchmark, so this is a significant limitation.
In contrast, for K-FAC Influence, we find that generally any sufficiently small value works reasonably well if enough samples are taken for estimating the loss and measurement gradients (see Figures 7 and 9).

Retraining without top influences. The counterfactual retraining results are shown in Figure 3 for CIFAR-2, CIFAR-10, with $2\%$ and $10\%$ of the data removed. In this evaluation, influence functions with K-FAC consistently pick more influential training examples (i.e. those which lead to a higher loss reduction) than the baselines.

Figure: Figure 3: Changes in measurements under counterfactual retraining without top influences for the loss measurement. The standard error in the estimate of the mean is indicated with error bars and reported after ‘$\pm$’, where the average is over different generated samples for which top influences are being identified.
Refer to caption: x4.png

Hessian Approximation Ablation. In Figure 4, we explore the impact of the Hessian approximation design choices discussed in Section 3.1.
We use K-FAC to approximate the GGN in all cases, with either the “expand” or the “reduce” variant (Section 3.1.2).
We find that the better-motivated “MC-Fisher” estimator $\operatorname{GGN}^{\texttt{model}}$ in Equation 9 does indeed perform better than the “empirical Fisher” in Equation 11 used in TRAK and D-TRAK.
Secondly, we find that K-FAC expand significantly outperforms K-FAC reduce, which stands in contrast to the results in the second-order optimisation setting where the two are on par with one another (Eschenhagen et al., ).
There are multiple differences from our setting to the one from the previous optimisation results: we use a square loss instead of a cross entropy loss, a full dataset estimate, a different architecture, and evaluate the approximation in a different application.
Notably, the expand variant is the better justified one since the diffusion modelling problem corresponds to the expand setting in Eschenhagen et al. ().
Hence, our results all seem to imply that a better Hessian approximation directly results in better downstream data attribution performance.
However, we do not directly evaluate the approximation quality of the estimates and also do not sweep over the damping value for all variants.

Figure: Figure 4: Ablation over the different Hessian approximation variants introduced in Section 3.1. We ablate two versions of the GGN: the “MC” Fisher in Equation 9 and the “Empirical” Fisher in Equation 11, as well as two settings for the K-FAC approximation: “expand” and “reduce” (Eschenhagen et al., ).
Refer to caption: x5.png

### 4.1 Potential challenges to use of influence functions for diffusion models

One peculiarity in the LDS results, similar to the findings in Zheng et al. (), is that substituting the loss measurement for the ELBO measurement when predicting changes in ELBO actually works better than using the correct measurement (see Figure 2(b) “K-FAC Influence (measurement loss)”).(^9^99Note that, unlike Zheng et al. (), we only change the measurement function for a proxy in the influence function approximation, keeping the Hessian approximation and training loss gradient in Equation 6 the same.)
To try and better understand the properties of influence functions for predicting the numerical change in behaviour after retraining, in this section we perform multiple ablations and report different interesting phenomena resulting from these that give some insight into the challenges of using influence functions in this setting.(^bordercolor=orange!50^bordercolor=orange!50todo: bordercolor=orange!50BM: maybe a slightly clunky sentence.)

As illustrated in Figure 18, gradients of the ELBO and training loss measurements, up to a constant scaling, consist of the same per-diffusion-timestep loss term gradients $\nabla_{\theta}\ell_{t}(\theta,x)$, but with a different weighting.
To try and break-down why approximating the change in ELBO with the training loss measurement gives higher LDS scores, we first look at predicting the change in the per-diffusion-timestep losses $\ell_{t}$ while substituting different per-diffusion-timestep losses into the K-FAC influence approximation.
The results are shown in Figure 11, leading to the following observation:

###### Observation 1

More specifically, changes in losses $\ell_{t}$ can in general be well approximated by substituting measurements $\ell_{t^{\prime}}$ into the influence approximation with $t^{\prime}>t$.
In some cases, using the incorrect timestep $t^{\prime}>t$ even results in significantly better LDS scores than the correct timestep $t^{\prime}=t$.

Based on 1, it is clear that influence function-based approximations have limitations when being applied to predict the numerical change in loss measurements.
We observe another pattern in how they can fail:

###### Observation 2

We show in Figures 15 and 16 that influence functions tend to overestimate how often removal of a group data points will lead to improvements in loss on a generated sample (both for aggregate diffusion training loss in Section 2.1, and the per-diffusion-timestep loss in Equation 2).

Lastly, although ELBO is perhaps the measurement with the most direct link to the marginal probability of sampling a particular example, we find it has some peculiar properties.
The below observation particularly puts its usefulness as a measurement function into question:

###### Observation 3

As illustrated in Figure 17, ELBO measurement is close to constant for any given sample generated from the model, no matter which $50\%$ subset of the training data is removed.
In particular, it is extremely rare that if one sample is more likely to be generated than another sample by one model (as measured by $\operatorname{ELBO}$), it will be less likely to be generated than another by a model trained on a different random subset of the data.
This implies that the standard DDPM $\operatorname{ELBO}$ might be a poor proxy for the real question we intend to ask: “how likely is a given sample to be generated by a model?”.
It would appear improbable that the true probabilities of generating any given sample are close-to unaffected by resampling what data subset the model was trained on.

###### Observation 1

###### Observation 2

###### Observation 3

## 5 Discussion

In this work, we extended the influence functions approach to the diffusion modelling setting, and showed different ways in which the GGN Hessian approximation can be formulated in this setting.
Our proposed method with recommended design choices improves performance copmared to existing techniques across various data attribution evaluation metrics.
Nonetheless, experimentally, we are met with two contrasting findings:
on the one hand, influence functions in the diffusion modelling setting appear to be able to identify important influences.
The surfaced influential examples do significantly impact the training loss when retraining the model without them (Figure 3), and they appear perceptually very relevant to the generated samples.
On the other hand, they fall short of accurately predicting the numerical changes in measurements after retraining.
Thes appears to be especially the case for measurement functions we would argue are most relevant in the image generative modelling setting – proxies for marginal probability of sampling a particular example.
This appears to be both due to the limitations of the influence functions approximation, but also due to the shortcomings of the considered proxy measurements (Section 4.1).

Despite these shortcomings, influence functions can still offer valuable insights: they can serve as a useful exploratory tool for understanding model behaviour in a diffusion modelling context, and can help guide data curation, identifying examples most responsible for certain behaviours.
To make them useful in settings where numerical accuracy in the predicted behaviour after retraining is required, such as copyright infringement, we believe more work is required into 1) finding better proxies for marginal probability than ELBO and probability of sampling trajectory , and 2) even further improving the influence function approximation.

#### Acknowledgments

We thank Jihao Andreas Lin for useful discussions on compression, and help with implementation of quantisation and SVD compression.
We would also like to thank Kristian Georgiev for sharing with us the diffusion model weights used for analysis in Georgiev et al. (), and Felix Dangel for help and feedback on implementing extensions to the [curvlinops](https://github.com/f-dangel/curvlinops) package used for the experiments in this paper.
Richard E. Turner is supported by Google, Amazon, ARM, Improbable, EPSRC grant EP/T005386/1, and the EPSRC Probabilistic AI Hub (ProbAI, EP/Y028783/1).
(^bordercolor=orange!50^bordercolor=orange!50todo: bordercolor=orange!50BM: Rich: would you like to keep this for the arxiv version?)

## Appendix A Derivation of Influence Functions

In this section, we state the implicit function theorem (Section A.1). Then, in Section A.2, we introduce the details of how it can be applied in the context of a loss function $\mathcal{L}(\bm{\varepsilon},\bm{\theta})$ parameterised by a continuous hyperparameter $\bm{\varepsilon}$ (which is, e.g., controlling how down-weighted the loss terms on some examples are, as in Section 2.2).

Figure: Figure 5: Illustration of the influence function approximation for a 1-dimensional parameter space $\theta\in\mathbb{R}$. Influence funcitons consider the extended loss landscape $\mathcal{L}(\varepsilon,\theta)\mathrel{\stackrel{{\scriptstyle\textnormal{def }}}{{=}}}\frac{1}{N}\sum_{n=1}^{N}\ell(x_{n},\theta)-\varepsilon\ell(x_{j},\theta)$, where the loss $\ell(x_{j},\theta)$ for some datapoint $x_{j}$ (alternatively, group of datapoints) is down-weighted by $\varepsilon$. By linearly extrapolating how the optimal set of parameters $\theta$ would change around $\varepsilon=0$ ( ), we can predicted how the optimal parameters would change when the term $\ell(x_{j},\theta)$ is fully removed from the loss ( ).
Refer to caption: x6.png

### A.1 Implicit Function Theorem

###### Theorem 1 (Implicit Function Theorem ( Krantz & Parks , 2003 ) )

Let $F:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{m}$ be a continuously differentiable function, and let $\mathbb{R}^{n}\times\mathbb{R}^{m}$ have coordinates $(\mathbf{x},\mathbf{y})$. Fix a point $(\mathbf{a},\mathbf{b})=(a_{1},\ldots,a_{n},b_{1},\ldots,b_{m})$ with $F(\mathbf{a},\mathbf{b})=\bm{0}$, where $\bm{0}\in\mathbb{R}^{m}$ is the zero vector. If the Jacobian matrix $\nabla_{\mathbf{y}}F(\mathbf{a},\mathbf{b})\in\mathbb{R}^{m\times m}$ of $\mathbf{y}\mapsto F(\mathbf{a},\mathbf{y})$

$$ $\left[\nabla_{\mathbf{y}}F(\mathbf{a},\mathbf{b})\right]_{ij}=\frac{\partial F _{i}}{\partial y_{j}}(\mathbf{a},\mathbf{b}),$ $$

is invertible, then there exists an open set $U\subset\mathbb{R}^{n}$ containing $\mathbf{a}$ such that there exists a unique function $g:U\to\mathbb{R}^{m}$ such that $g(\mathbf{a})=\mathbf{b}$, and $F(\mathbf{x},g(\mathbf{x}))=\bm{0}$ for all $\mathbf{x}\in U$. Moreover, $g$ is continuously differentiable.

###### Remark 1 (Derivative of the implicit function)

Denoting the Jacobian matrix of $\mathbf{x}\mapsto F(\mathbf{x},\mathbf{y})$ as:

$$ $\left[\nabla_{\mathbf{x}}F(\mathbf{x},\mathbf{y})\right]_{ij}=\frac{\partial F _{i}}{\partial x_{j}}(\mathbf{x},\mathbf{y}),$ $$

the derivative $\frac{\partial\mathbf{g}}{\partial\mathbf{x}}:U\to\mathbb{R}^{m\times n}$ of $g:U\to\mathbb{R}^{m}$ in Theorem 1 can be written as:

$$ $\frac{\partial g(\mathbf{x})}{\partial\mathbf{x}}=-\left[\nabla_{\mathbf{y}}F( \mathbf{x},g(\mathbf{x}))\right]^{-1}\nabla_{\mathbf{x}}F(\mathbf{x},g(\mathbf {x})).$ (14) $$

This can readily be seen by noting that, for $\mathbf{x}\in U$:

$$ $F(\mathbf{x^{\prime}},g(\mathbf{x^{\prime}}))=\bm{0}\quad\forall\mathbf{x}^{ \prime}\in U\qquad\Rightarrow\qquad\frac{dF(\mathbf{x},g(\mathbf{x}))}{d \mathbf{x}}=\bm{0}.$ $$

Hence, since $g$ is differentiable, we can apply the chain rule of differentiation to get:

$$ $\bm{0}=\frac{dF(\mathbf{x},g(\mathbf{x}))}{d\mathbf{x}}=\nabla_{\mathbf{x}}F( \mathbf{x},g(\mathbf{x}))+\nabla_{\mathbf{y}}F(\mathbf{x},g(\mathbf{x}))\frac{ \partial g(\mathbf{x})}{\partial\mathbf{x}}.$ $$

Rearranging gives equation Equation 14.

###### Theorem 1 (Implicit Function Theorem ( Krantz & Parks , 2003 ) )

###### Remark 1 (Derivative of the implicit function)

### A.2 Applying the implicit function theorem to quantify the change in the optimum of a loss

Consider a loss function $\mathcal{L}:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ that depends on some hyperparameter $\bm{\varepsilon}\in\mathbb{R}^{n}$ (in Section 2.2, this was the scalar by which certain loss terms were down-weighted) and some parameters $\bm{\theta}\in\mathbb{R}^{m}$.
At the minimum of the loss function $\mathcal{L}(\bm{\varepsilon},\bm{\theta})$, the derivative with respect to the parameters $\bm{\theta}$ will be zero.
Hence, assuming that the loss function is twice continuously differentiable (hence $\frac{\partial{L}}{\partial\bm{\varepsilon}}$ is continuously differentiable), and assuming that for some $\bm{\varepsilon}^{\prime}\in\mathbb{R}^{n}$ we have a set of parameters $\bm{\theta}^{\star}$ such that $\frac{\partial\mathcal{L}}{\partial\bm{\varepsilon}}(\bm{\varepsilon}^{\prime}
,\bm{\theta}^{\star})=\bm{0}$ and the Hessian $\frac{\partial^{2}\mathcal{L}}{\partial\bm{\theta}^{2}}(\bm{\varepsilon}^{
\prime},\bm{\theta}^{\star})$ is invertible, we can apply the implicit function theorem to the derivative of the loss function $\frac{\partial\mathcal{L}}{\partial\bm{\varepsilon}}:\mathbb{R}^{n}\times
\mathbb{R}^{m}\to\mathbb{R}^{m}$, to get the existence of a continuously differentiable function $g$ such that $\frac{\partial\mathcal{L}}{\partial\bm{\varepsilon}}(\bm{\varepsilon},g(\bm{
\varepsilon}))=\bm{0}$ for $\bm{\varepsilon}$ in some neighbourhood of $\bm{\varepsilon}^{\prime}$.

Now $g(\bm{\varepsilon})$ might not necessarily be a minimum of $\bm{\theta}\mapsto\mathcal{L}(\bm{\varepsilon},\bm{\theta})$. However, by making the further assumption that $\mathcal{L}$ is strictly convex we can ensure that whenever $\frac{\partial\mathcal{L}}{\partial\bm{\theta}}(\bm{\varepsilon},\bm{\theta})=
\bm{0}$, $\bm{\theta}$ is a unique minimum, and so $g(\bm{\varepsilon})$ represents the change in the minimum as we vary $\mathbf{\varepsilon}$.
This is summarised in the lemma below:

###### Lemma 1

Let $\mathcal{L}:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ be a twice continuously differentiable function, with coordinates denoted by $(\bm{\varepsilon},\bm{\theta})\in\mathbb{R}^{n}\times\mathbb{R}^{m}$, such that $\mathbf{\theta}\mapsto\mathcal{L}(\bm{\varepsilon},\bm{\theta})$ is strictly convex $\forall\bm{\varepsilon}\in\mathbb{R}^{n}$. Fix a point $(\bm{\varepsilon}^{\prime},\bm{\theta}^{\star})$ such that $\frac{\partial\mathcal{L}}{\partial\bm{\theta}}(\bm{\varepsilon}^{\prime},\bm{
\theta}^{\star})=\bm{0}$. Then, by the Implicit Function Theorem applied to $\frac{\partial\mathcal{L}}{\partial\bm{\theta}}$, there exists an open set $U\subset\mathbb{R}^{n}$ containing $\bm{\theta}^{\star}$ such that there exists a unique function $g:U\to\mathbb{R}^{m}$ such that $g(\bm{\varepsilon}^{\prime})=\bm{\theta}^{\star}$, and $g(\bm{\varepsilon})$ is the unique minimum of $\mathbf{\theta}\mapsto\mathcal{L}(\bm{\varepsilon},\bm{\theta})$ for all $\bm{\varepsilon}\in U$. Moreover, $g$ is continuously differentiable with derivative:

$$ $\frac{\partial g(\bm{\varepsilon})}{\partial\bm{\varepsilon}}=-\left[\frac{ \partial^{2}\mathcal{L}}{\partial\bm{\theta}^{2}}(\bm{\varepsilon},g(\bm{ \varepsilon}))\right]^{-1}\frac{\partial^{2}\mathcal{L}}{\partial\bm{ \varepsilon}\partial\bm{\theta}}(\bm{\varepsilon},g(\bm{\varepsilon}))$ (15) $$

###### Remark 2

For a loss function $\mathcal{L}:\mathbb{R}\times\mathbb{R}^{m}$ of the form $\mathcal{L}(\varepsilon,\bm{\theta})=\mathcal{L}_{1}(\bm{\theta})+\varepsilon
\mathcal{L}_{2}(\bm{\theta})$ (such as that in Equation 4), $\frac{\partial^{2}\mathcal{L}}{\partial\bm{\varepsilon}\partial\bm{\theta}}(
\bm{\varepsilon},g(\bm{\varepsilon}))$ in the equation above simplifies to:

$$ $\frac{\partial^{2}\mathcal{L}}{\partial\varepsilon\partial\bm{\theta}}( \varepsilon,g(\varepsilon))=\frac{\partial\mathcal{L}_{2}}{\partial\bm{\theta} }(g({\varepsilon}))$ (16) $$

The above lemma and remark give the result in Equation 5.
Namely, in section 2.2:

$$ $\displaystyle\mathcal{L}(\varepsilon,\bm{\theta})=\underbrace{\frac{1}{N}\sum_ {i=1}^{N}\ell(\bm{\theta},x_{i})}_{\mathcal{L}_{1}}\overbrace{-\frac{1}{M}\sum _{j=1}^{M}\ell(\bm{\theta},x_{i_{j}})}^{\mathcal{L}_{2}}\varepsilon\qquad \overset{\text{\lx@cref{creftype~refnum}{eq:weighted-loss-as-sum-gradient}}}{ \Longrightarrow}\qquad\frac{\partial^{2}\mathcal{L}}{\partial\varepsilon \partial\bm{\theta}}=-\frac{1}{M}\sum_{j=1}^{M}\frac{\partial}{\partial\bm{ \theta}}\ell(\bm{\theta},x_{i_{j}})$ $\displaystyle\overset{\text{\lx@cref{creftype~refnum}{eq:lemma-loss-implicit-f unction}}}{\Longrightarrow}\qquad\frac{\partial g(\bm{\varepsilon})}{\partial \bm{\varepsilon}}=\left[\frac{\partial^{2}\mathcal{L}}{\partial\bm{\theta}^{2} }(\bm{\varepsilon},g(\bm{\varepsilon}))\right]^{-1}\frac{1}{M}\sum_{j=1}^{M} \frac{\partial}{\partial\bm{\theta}}\ell(\bm{\theta},x_{i_{j}})$ $$

###### Lemma 1

###### Remark 2

## Appendix B Derivation of the Fisher “GGN” formulation for Diffusion Models

As discussed in Section 2.2.1 partitioning the function $\theta\mapsto\|{\epsilon}^{(t)}-{\epsilon}_{\theta}^{t}(x^{(t)})\|^{2}$ into the model output $\theta\mapsto{\epsilon}_{\theta}^{t}(x^{(t)})$ and the $\ell_{2}$ loss function is a natural choice and results in

$$ $\displaystyle\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}(\theta)$ $\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\mathbb{E}_{\tilde{t}}\left[\mathbb{E}_ {x^{{(\tilde{t})}},{\epsilon}^{(\tilde{t})}}\left[\nabla^{\mathchoice{\raise 1 .3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle \intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt \hbox{$\scriptscriptstyle\intercal$}}}_{\theta}{\epsilon_{\theta}^{\tilde{t}} \left(x^{(\tilde{t})}\right)}\nabla_{{\epsilon_{\theta}^{\tilde{t}}\left(x^{( \tilde{t})}\right)}}^{2}{\left\lVert{\epsilon}^{(\tilde{t})}-{\epsilon_{\theta }^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\right\rVert^{2}}\nabla_{\theta}{ \epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\right]\right]$ $\displaystyle=\frac{2}{N}\sum_{n=1}^{N}\mathbb{E}_{\tilde{t}}\left[\mathbb{E}_ {x^{(\tilde{t})},{\epsilon}^{(\tilde{t})}}\left[\nabla^{\mathchoice{\raise 1.3 611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle \intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt \hbox{$\scriptscriptstyle\intercal$}}}_{\theta}{\epsilon}_{\theta}^{\tilde{t}} \left(x^{(\tilde{t})}\right)I\nabla_{\theta}{\epsilon}_{\theta}^{\tilde{t}} \left(x^{(\tilde{t})}\right)\right]\right].$ (17) $$

Note that we used

$$ $\displaystyle\frac{1}{2}\nabla_{{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde {t})}\right)}}^{2}{\left\lVert{\epsilon}^{(\tilde{t})}-{\epsilon_{\theta}^{ \tilde{t}}\left(x^{(\tilde{t})}\right)}\right\rVert^{2}}=I.$ $$

We can substitute $I$ with

$$ $\displaystyle I=\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[-\frac{1}{2}\nabla_{ {\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}^{2}\log p\left({ \epsilon_{\texttt{mod}}}|{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})} \right)}\right)\right],$ $\displaystyle p\left({\epsilon_{\texttt{mod}}}|{\epsilon_{\theta}^{\tilde{t}} \left(x^{(\tilde{t})}\right)}\right)={\mathcal{N}}\left({\epsilon_{\texttt{mod }}}|{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)},I\right),$ $$

where the mean of the Gaussian is chosen to be the model output ${\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}$.
Furthermore, by using the “score” trick:

$$ $\displaystyle\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[\nabla_{{\epsilon_{ \theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}^{2}\log p\left({\epsilon_{ \texttt{mod}}}|{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)} \right)\right]$ $\displaystyle=-\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[\nabla_{{\epsilon_{ \theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}\log p\left({\epsilon_{ \texttt{mod}}}|{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)} \right)\nabla^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$}}{ \raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$ \scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}} }_{{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}\log p\left({ \epsilon_{\texttt{mod}}}|{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})} \right)}\right)\right]$ $\displaystyle=-\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[\frac{1}{2}\nabla_{{ \epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}\left\lVert{ \epsilon_{\texttt{mod}}}-{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})} \right)}\right\rVert^{2}\frac{1}{2}\nabla^{\mathchoice{\raise 1.3611pt\hbox{$ \displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0 .95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$ \scriptscriptstyle\intercal$}}}_{{\epsilon_{\theta}^{\tilde{t}}\left(x^{( \tilde{t})}\right)}}\left\lVert{\epsilon_{\texttt{mod}}}-{\epsilon_{\theta}^{ \tilde{t}}\left(x^{(\tilde{t})}\right)}\right\rVert^{2}\right],$ $$

we can rewrite:

$$ $\displaystyle\nabla^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$ }}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$ \scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}} }_{\theta}{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\nabla_{ \theta}{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}$ $\displaystyle=-2\nabla^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle \intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt \hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle \intercal$}}}_{\theta}{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})} \right)}\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[\left(\nabla_{{\epsilon_{ \theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}^{2}\log p\left({\epsilon_{ \texttt{mod}}}|{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)} \right)\right)\right]\nabla_{\theta}{\epsilon_{\theta}^{\tilde{t}}\left(x^{( \tilde{t})}\right)}$ $\displaystyle=\frac{1}{2}\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[\nabla^{ \mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt \hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{ \raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}}}_{\theta}{\epsilon_{ \theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\nabla_{{\epsilon_{\theta}^{ \tilde{t}}\left(x^{(\tilde{t})}\right)}}\left\lVert{\epsilon_{\texttt{mod}}}-{ \epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\right\rVert^{2} \nabla^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3 611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle \intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}}}_{{\epsilon _{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}}\left\lVert{\epsilon_{ \texttt{mod}}}-{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)} \right\rVert^{2}\nabla_{\theta}{\epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{ t})}\right)}\right]$ $\displaystyle=\frac{1}{2}\mathbb{E}_{\epsilon_{\texttt{mod}}}\left[\nabla_{ \theta}\left\lVert{\epsilon_{\texttt{mod}}}-{\epsilon_{\theta}^{\tilde{t}} \left(x^{(\tilde{t})}\right)}\right\rVert^{2}\nabla^{\mathchoice{\raise 1.3611 pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$} }{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$ \scriptscriptstyle\intercal$}}}_{\theta}\left\lVert{\epsilon_{\texttt{mod}}}-{ \epsilon_{\theta}^{\tilde{t}}\left(x^{(\tilde{t})}\right)}\right\rVert^{2} \right],$ $$

where the last equality follows by the chain rule of differentiation.
We can thus rewrite the expression for the $\mathrm{GGN}$ in Equation 17 as

$$ $\displaystyle\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}(\theta)$ $\displaystyle=\frac{1}{N}\sum_{n=1}^{N}\mathbb{E}_{\tilde{t}}\left[\mathbb{E}_ {x^{{(\tilde{t})}},{\epsilon}^{(\tilde{t})},{\epsilon_{\texttt{mod}}}}\left[ \nabla_{\theta}g_{n}(\theta)\nabla_{\theta}g_{n}(\theta)^{\mathchoice{\raise 1 .3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle \intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt \hbox{$\scriptscriptstyle\intercal$}}}\right]\right]$ $\displaystyle g(\theta)\mathrel{\stackrel{{\scriptstyle\textnormal{def}}}{{=}} }\left\lVert{\epsilon_{\texttt{mod}}}-{\epsilon_{\theta}^{\tilde{t}}\left(x^{( \tilde{t})}\right)}\right\rVert^{2}.$ $$

## Appendix C Gradient compression ablation

Figure: Figure 6: Comparison of gradient compression methods for the influence function approximation.
Refer to caption: x7.png

In Figure 6, we ablate different compression methods by computing the per training datapoint influence scores with compressed query (measurement) gradients, and looking at the Pearson correlation and the rank correlation to the scores compute with the uncompressed gradients.
We hope to see a correlation of close to $100\%$, in which case the results for our method would be unaffected by compression.
We find that using quantisation for compression results in almost no change to the ordering over training datapoints, even when quantising down to $8$ bits.
This is in contrast to the SVD compression scheme used in Grosse et al. ().
This is likely because the per-example gradients naturally have a low-rank (Kronecker) structure in the classification, regression, or autoregressive language modelling settings, such as that in Grosse et al. ().
On the other hand, the diffusion training loss and other measurement functions considered in this work do not have this low-rank structure.
This is because computing them requires multiple forward passes; for example, for the diffusion training loss we need to average the mean-squared error loss in Equation 2 over multiple noise samples $\epsilon^{(t)}$ and multiple diffusion timesteps.
We use $8$ bit quantisation with query gradient batching (Grosse et al., ) for all KFAC experiments throughout this work.

## Appendix D Damping LDS ablations

We report an ablation over the LDS scores with GGN approximated with different damping factors for TRAK/D-TRAK and K-FAC influence in Figures 7, 8, 9 and 10. The reported damping factors for TRAK are normalised by the dataset size so that they correspond to the equivalent damping factors for our method when viewing TRAK as an altenrative approximation to the GGN (see Section 3.1).

Figure: Figure 7: Effect of damping on the LDS scores for K-FAC influence on CIFAR-2. In this plot, K-FAC GGN approximation was always computed with $1000$ samples, and the number of samples used for computing a Monte Carlo estimate of the training loss/measurement gradient is indicated on the legend.
Refer to caption: x8.png

Figure: Figure 8: Effect of damping on the LDS scores for TRAK (random projection) based influence on CIFAR-2. $250$ samples were used for Monte Carlo estiamtion of all quantities (GGN and the training loss/measurement gradients). In the legend: Target indicates what measurement we’re trying to predict the change in after retraining, Measure indicates what measurement function was substituted into the influence function approximation, and Train.Loss indicates what function was substituted for the training loss in the computation of the GGN and gradient of the training loss in the influence function approximation.
Refer to caption: x9.png

Figure: Figure 9: Effect of damping on the LDS scores for K-FAC based influence on CIFAR-10. $100$ samples were used for computing the K-FAC GGN approximation, and $250$ for computing a Monte Carlo estimate of the training loss/measurement gradients. $\bm{\times}$ indicates a NaN result (the computation was not sufficiently numerically stable with that damping factor).
Refer to caption: x10.png

Figure: Figure 10: Effect of damping on the LDS scores for TRAK (random projection) based influence on CIFAR-10. $250$ samples were used for Monte Carlo estiamtion of all quantities (GGN and the training loss/measurement gradients). In the legend: Target indicates what measurement we’re trying to predict the change in after retraining, Measure indicates what measurement function was substituted into the influence function approximation, and Train.Loss indicates what function was substituted for the training loss in the computation of the GGN and gradient of the training loss in the influence function approximation.
Refer to caption: x11.png

## Appendix E Empirical ablations for challenges to use of influence functions for diffusion models

In this section, we describe the results for the observations discussed in Section 4.1.

1 is based on Figures 11 and 12.
Figure 11 shows the LDS scores on CIFAR-2 when attributing per-timestep diffusion losses $\ell_{t}$ (see Equation 2) using influence functions, whilst varying what (possibly wrong) per-timestep diffusion loss $\ell_{t^{\prime}}$ is used as a measurement function in the influence function approximation (Equation 6).
Figure 12 is a counter-equivalent to Figure 16 where instead of using influence functions to approximate the change in measurement, we actually retrain a model on the randomly subsampled subset of data and compute the measurement.

Figure: Figure 11: Rank correlation (LDS scores) between influence function estimates with different measurement functions and different true measurements CIFAR-2. The plot shows how well different per-timestep diffusion losses $\ell_{t}$ work as measurement functions in the influence function approximation, when trying to approximate changes in the actual measurements when retraining a model.
Refer to caption: x12.png

A natural question to ask with regards to 1 is: does this effect go away in settings where the influence function approximation should more exact?
Note that, bar the non-convexity of the training loss function ${\mathcal{L}}_{\mathcal{D}}$, the influence function approximation in Equation 6 is a linearisation of the actual change in the measurement for the optimum of the training loss functions with some examples down-weighted by $\varepsilon$ around $\varepsilon=0$.
Hence, we might expect the approximation to be more exact when instead of fully removing some data points from the dataset (setting $\varepsilon=\nicefrac{{1}}{{N}}$), we instead down-weight their contribution to the training loss by a smaller non-zero factor.
To investigate whether this is the case, we repeat the LDS analysis in Figures 11 and 12, but with $\varepsilon=\nicefrac{{1}}{{2N}}$; in other words, the training loss terms corresponding to the “removed” examples are simply down-weighted by a factor of $\nicefrac{{1}}{{2}}$ in the retrained models.
The results are shown in Figures 13 and 14.
Perhaps somewhat surprisingly, a contrasting effect can be observed, where using per-timestep diffusion losses for larger times yields a higher absolute rank correlation, but with the opposing sign.
The negative correlation between measurement $\ell_{t},\ell_{t^{\prime}}$ for $t\neq t^{\prime}$ can also be observed for the true measurements in the retrained models in Figure 14.
We also observe that in this setting, influence functions fail completely to predict changes in $\ell_{t}$ with the correct measurement function for $t\leq 200$.

Figure: Figure 12: Rank correlation between true measurements for losses at different diffusion timesteps on CIFAR-2.
Refer to caption: x13.png

Figure: Figure 13: Rank correlation (LDS scores) between influence function estimates with different measurement functions and different true measurements CIFAR-2, but with the retrained models trained on the full dataset with a random subset of examples having a down-weighted contribution to a training loss by a factor of $\times 0.5$.
Refer to caption: x14.png

Figure: Figure 14: Rank correlation between true measurements for losses at different diffusion timesteps on CIFAR-2, but with the retrained models trained on the full dataset with a random subset of examples having a down-weighted contribution to a training loss by a factor of $\times 0.5$.
Refer to caption: x15.png

1 Figure 15 shows the changes in losses after retraining the model on half the data removed against the predicted changes in losses using K-FAC Influence for two datasets: CIFAR-2 and CIFAR-10.
In both cases, for a vast majority of retrained models, the loss measurement on a sample increases after retraining.
On the other hand, the influence functions predict roughly evenly that the loss will increase and decrease.
This trend is amplified if we instead look at influence predicted for per-timestep diffusion losses $\ell_{t}$ (Equation 2) for earlier timesteps $t$, which can be seen in Figure 16.
On CIFAR-2, actual changes in $\ell_{1},\ell_{5}0,\ell_{1}00$ measurements are actually always positive, which the influence functions approximation completely misses.
For all plots, K-FAC Influence was ran with a damping factor of $10^{-8}$ and $250$ samples for all gradient computations.

Figure: Figure 15: Change in diffusion loss $\ell$ in Section 2.1 when retraining with random subsets of $50\%$ of the training data removed, as predicted by K-FAC influence ($x$-axis), against the actual change in the measurement ($y$-axis). Results are plotted for measurements $\ell(x,\theta)$ for $50$ samples $x$ generated from the diffusion model trained on all of the data. The scatter color indicates the sample $x$ for which the change in measurement is plotted. The figure shows that influence functions tend to overestimate how often the loss will decrease when some training samples are removed; in reality, it happens quite rarely.
Refer to caption: x16.png

Figure: Figure 16: Change in per-diffusion-timestep losses $\ell_{t}$ when retraining with random subsets of $50\%$ of the training data removed, as predicted by K-FAC influence ($x$-axis), against the actual change in the measurement ($y$-axis). Results are plotted for the CIFAR-2 dataset, for measurements $\ell_{t}(x,\theta)$ for $50$ samples $x$ generated from the diffusion model trained on all of the data. The scatter color indicates the sample $x$ for which the change in measurement is plotted. The figure shows that: 1) influence functions predict that the losses $\ell_{t}$ will increase or decrease roughly equally frequently when some samples are removed, but, in reality, the losses almost always increase; 2) for sufficiently large time-steps ($\ell_{5}00$), this pattern seems to subside. Losses $\ell_{t}$ in the $\leavevmode\nobreak\ 200-500$ range seem to work well for predicting changes in other losses Figure 11.
Refer to caption: x17.png

3 Lastly, the observations that the ELBO measurements remain essentially constant for models trained on different subsets of data is based on Figure 17.
There, we plot the values of the ELBO measurement for different pairs of models trained on different subsets of data, where we find near perfect correlation.
The only pairs of models that exhibit an ELBO measurement correlation of less that $0.99$ are the CIFAR-2 model trained on the full dataset compared to any model trained on a $50\%$ subset, which is likely due to the fact that the $50\%$ subset models are trained for half as many gradient iterations, and so may have not fully converged yet.
For CIFAR-10, where we train for $5\times$ as many training steps due to a larger dataset size, we observe near-perfect correlation in the ELBO measurements across all models.
Each ELBO measurement was computed with a Monte-Carlo estimate using $5000$ samples.

Figure: Figure 17: Correlation of the $\operatorname{ELBO}(x,\theta)$ measurements on different data points $x$ (samples generated from the model trained on full data), for models trained on different subsets of data. Each subplot plots $\operatorname{ELBO}(x,\theta)$ measurements for $200$ generated samples $x$, as measured by two models trained from scratch on different subsets of data, with the $x$-label and the $y$-label identifying the respective split of data used for training (either full dataset, or randomly subsampled $50\%$-subset). Each subplot shows the Pearson correlation coefficient ($r$) and the Spearman rank correlation ($\rho$) for the $\operatorname{ELBO}(x,\theta)$ measurements as measured by the two models trained on different subsets of data. The two parts of the figure show results for two different datasets: CIFAR-2 on the left, and CIFAR-10 on the right.
Refer to caption: x18.png

Figure: Figure 18: The diffusion loss and diffusion ELBO as formulated in (Ho et al., ) (ignoring the reconstruction term that accounts for the quantisation of images back to pixel space) are equal up to the weighting of the individual per-diffusion-timestep loss terms and a constant independent of the parameters. This plot illustrates the relatives difference in the weighting for per-diffusion-timestep losses applied in the ELBO vs. in the training loss.
Refer to caption: x19.png

## Appendix F LDS results for probability of sampling trajectory

The results for the “log probability of sampling trajectory” measurements are shown in Figure 19.
The probability of sampling trajectory appears to be a measurement with a particularly low correlation across different models trained with the same data, but different random seeds.
This is perhaps unsurprising, since the measurement comprises the log-densities of particular values of $1000$ latent variables.

Figure: Figure 19: Linear Data-modelling Score (LDS) for the probability of sampling trajectory. The plot follows the same format as that of Figures 2(a) and 2(b). Overall, probability of the sampling trajectory appears to be a difficult proxy for the marginal probability of sampling a given example, given that it suffers from the same issues as the ELBO on CIFAR-2 (it’s better approximated by the wrong measurement function), and there is extremely little correlation in the measurement across the retrained models on larger datasets (CIFAR-10).
Refer to caption: x20.png

## Appendix G Experimental details

In this section, we describe the implementation details for the methods and baselines, as well as the evaluations reported in Section 4.

### G.1 Datasets

We focus on the following dataset in this paper:

CIFAR-10 CIFAR-10 is a dataset of small RGB images of size $32\times 32$ Krizhevsky (). We use $50000$ images (the train split) for training.

CIFAR-2 For CIFAR-2, we follow Zheng et al. () and create a subset of CIFAR-10 with $5000$ examples of images only corresponding to classes car and horse. $2500$ examples of class car and $2500$ examples of class horse are randomly subsampled without replacement from among all CIFAR-10 images of that class.

### G.2 Models

For all CIFAR datasets, we train a regular Denoising Diffusion Probabilistic Model using a standard U-Net architecture as described for CIFAR-10 in Ho et al. ().
This U-Net architecture contains both convolutional and attention layers.
We use the same noise schedule as described for the CIFAR dataset in Ho et al. ().

Sampling We follow the standard DDPM sampling procedure with a full $1000$ timesteps to create the generated samples as described by Ho et al. ().
DDPM sampling usually gives better samples (in terms of visual fidelity) than Denoising Diffusion Implicit Models (DDIM) sampling Song et al. () when a large number of sampling steps is used.
As described in Section 2.1, when parameterising the conditionals $p_{\theta}(x^{(t-1)}|x^{(t)})$ with neural networks as $\mathcal{N}\big{(}x^{(t-1)}|\mu_{t-1|t,0}\big{(}x^{(t)},{\epsilon}^{t}_{\theta
}(x^{(t)})\big{)},\sigma_{t}^{2}I\big{)}$ we have a choice in how to set the variance hyperparameters $\{\sigma_{t}^{2}\}_{t=1}^{T}$.
The $\sigma_{t}^{2}$ hyperparameters do not appear in the training loss; however, they do make a difference when sampling.
We use the “small” variance variant from Ho et al. (, §3.2), i.e. we set:

$$ $\sigma_{t}^{2}=\frac{1-\prod_{t^{\prime}=1}^{t-1}\lambda_{t^{\prime}}}{1-\prod _{t^{\prime}=1}^{t}\lambda_{t^{\prime}}}\left(1-\lambda_{t}\right)$ $$

### G.3 Details on data attribution methods

TRAK For TRAK baselines, we adapt the implementation of Park et al. (); Georgiev et al. () to the diffusion modelling setting.
When running TRAK, there are several settings the authors recommend to consider: 1) the projection dimension $d_{\texttt{proj}}$ for the random projections, 2) the damping factor $\lambda$, and 3) the numerical precision used for storing the projected gradients.
For (1), we use a relatively large projection dimension of $32768$ as done in most experiments in Zheng et al. ().
We found that the projection dimension affected the best obtainable results significantly, and so we couldn’t get away with a smaller one.
We also found that using the default float16 precision in the TRAK codebase for (3) results in significantly degraded results (see Figure 20, and so we recommend using float32 precision for these methods for diffusion models.
In all experiments, we use float32 throughout.
For the damping factor, we report the sweeps over LDS scores in Figures 8 and 10, and use the best result in each benchmark, as these methods fail drastically if the damping factor is too small.
The damping factor reported in the plots is normalised by the dataset size $N$, to match the definition of the GGN, and to make it comparable with the damping reported for other influence functions methods introduced in this paper.
For non-LDS experiments, we use the best damping value from the corresponding LDS benchmark.

CLIP cosine similarity One of the data attribution baselines used for the LDS experiments is CLIP cosine similarity (Radford et al., ).
For this baseline, we compute the CLIP embeddings (Radford et al., ) of the generated sample and training datapoints, and consider the cosine similarity between the two as the “influence” of that training datapoint on that particular target sample.
See (Park et al., ) for details of how this influence is aggregated for the LDS benchmark.
Of course, this computation does not in any way depend on the diffusion model or the measurement function used, so it is a pretty naïve method for estimating influence.

K-FAC We build on the [https://github.com/f-dangel/curvlinops](https://github.com/f-dangel/curvlinops) package for our implementation of K-FAC for diffusion models.
Except for the ablation in Figure 4, we use the K-FAC expand variant throughout.
We compute K-FAC for PyTorch nn.Conv2d and nn.Linear modules (including in attention), ignoring the parameters in the normalisation layers.

Compression for all K-FAC influence functions results, we use int8 quantisation for the query gradients.

Monte Carlo computation of gradients and the GGN for influence functions Computing the per-example training loss $\ell(\theta,x_{n})$ in Section 2.1, the gradients of which are necessary for computing the influence function approximation (Equation 6), includes multiple nested expectations over diffusion timestep $\tilde{t}$ and noise added to the data $\epsilon^{(t)}$.
This is also the case for the $\mathrm{GGN}^{\texttt{model}}_{\mathcal{D}}$ in Equation 9 and for the gradients $\nabla_{\theta}\ell(\theta,x_{n})$ in the computation of $\mathrm{GGN}^{\texttt{loss}}_{\mathcal{D}}$ in Equation 11, as well as for the computation of the measurement functions.
Unless specified otherwise, we use the same number of samples for a Monte Carlo estimation of the expectations for all quantities considered.
For example, if we use $K$ samples, that means that for the computation of the gradient of the per-example-loss $\nabla_{\theta}\ell(\theta,x_{n})$ we’ll sample tuples of $(\tilde{t},\epsilon^{(\tilde{t})},x^{(\tilde{t})})$ independently $K$ times to form a Monte Carlo estimate.
For $\mathrm{GGN}^{\texttt{model}}_{\mathcal{D}}$, we explicitly iterate over all training data points, and draw $K$ samples of $\left(\tilde{t},\epsilon^{(\tilde{t})},x^{(\tilde{t})}_{n}\right)$ for each datapoint.
For $\mathrm{GGN}^{\texttt{loss}}_{\mathcal{D}}$, we explicitly iterate over all training data points, and draw $K$ samples of $\left(\tilde{t},\epsilon^{(\tilde{t})},x^{(\tilde{t})}_{n}\right)$ to compute the gradients $\nabla_{\theta}\ell(\theta,x_{n})$ before taking an outer product.
Note that, for $\mathrm{GGN}^{\texttt{loss}}_{\mathcal{D}}$, because we’re averaging over the samples before taking the outer product of the gradients, the estimator of the GGN is no longer unbiased.
Similarly, $K$ samples are also used for computing the gradients of the measurement function.

For all CIFAR experiments, we use $250$ samples throughout for all methods (including all gradient and GGN computations for K-FAC Influence, TRAK, D-TRAK), unless explicitly indicated in the caption otherwise.

Figure: Figure 20: LDS scores on for TRAK (random projection) based influence on CIFAR-2 when using half-precision (float16) for influence computations. Compare with Figure 8. NaN results are indicated with $\bm{\times}$.
Refer to caption: x21.png

### G.4 Damping

For all influence function-like methods (including TRAK and D-TRAK), we use damping to improve the numerical stability of the Hessian inversion. Namely, for any method that computes the inverse of the approximation to the Hessian $H\approx\nabla_{\theta}^{2}{\mathcal{L}}_{\mathcal{D}}=\nabla_{\theta}^{2}
\nicefrac{{1}}{{N}}\sum\ell(theta,x_{n})$, we add a damping factor $\lambda$ to the diagonal before inversion:

$$ $(H+\lambda I)^{-1},$ $$

where $I$ is a ${d_{\texttt{param}}\times d_{\texttt{param}}}$ identity matrix.
This is particularly important for methods where the Hessian approximation is at a high risk of being low-rank (for example, when using the empirical GGN in Equation 11, which is the default setting for TRAK and D-TRAK).
For TRAK/D-TRAK, the approximate Hessian inverse is computed in a smaller projected space, and so we add $\lambda$ to the diagonal directly in that projected space, as done in Zheng et al. ()).
In other words, if $P\in\mathbb{R}^{d_{\texttt{proj}}\times d_{\texttt{param}}}$ is the projection matrix (see (Park et al., ) for details), then damped Hessian-inverse preconditioned vector inner products between two vectors $v_{1},v_{2}\in\mathbb{R}^{d_{\texttt{param}}}$ (e.g. the gradients in Equation 6) would be computed as:

$$ $\left(Pv_{1}\right)^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle\intercal$ }}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt\hbox{$ \scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle\intercal$}} }\left(H+\lambda I\right)^{-1}Pv_{,}.$ $$

where $H\approx P\nabla_{\theta}^{2}{\mathcal{L}}_{\mathcal{D}}P^{\mathchoice{\raise 1
.3611pt\hbox{$\displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle
\intercal$}}{\raise 0.95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt
\hbox{$\scriptscriptstyle\intercal$}}}\in\mathbb{R}^{d_{\texttt{proj}}\times{d
_{\texttt{proj}}}}$ is an approximation to the Hessian in the projected space.

For the “default” values used for damping for TRAK, D-TRAK and K-FAC Influence, we primarily follow recommendations from prior work.
For K-FAC Influence, the default is a small damping value $10^{-8}$ throughout added for numerical stability of inversion, as done in prior work (Bae et al., ).
For TRAK-based methods, Park et al. () recommend using no damping. Hence, we use the lowest numerically stable value of $\leavevmode\nobreak\ 10^{-9}$ as the default value throughout.

Note that all damping values reported in this paper are reported as if being added to the GGN for the Hessian of the loss normalised by dataset size . This differs from the damping factor in the TRAK implementation ([https://github.com/MadryLab/trak](https://github.com/MadryLab/trak)), which is added to the GGN for the Hessian of an unnormalised loss ($\sum_{n}\ell(\theta,x_{n})$). Hence, the damping values reported in (Zheng et al., ) are larger by a factor of $N$ (the dataset size) than the equivalent damping values reported in this paper.

### G.5 LDS Benchmarks

For all LDS benchmarks Park et al. (), we sample $100$ sub-sampled datasets ($M\coloneq 100$ in Equation 13), and we train $5$ models with different random seeds ($K\coloneq 5$ in Equation 13), each with $50\%$ of the examples in the full dataset, for a total of $500$ retrained models for each benchmark.
We compute the LDS scores for $200$ samples generated by the model trained on the full dataset.

Monte Carlo sampling of measurements For all computations of the “true” measurement functions for the retrained models in the LDS benchmarks we use $5000$ samples to estimate the measurement.

### G.6 Retraining without top influences

For the retraining without top influences experiments (Figure 3), we pick $5$ samples generated by the model trained on the full dataset, and, for each, train a model with a fixed percentage of most influential examples for that sample removed from the training dataset, using the same procedure as training on the full dataset (with the same number of training steps).
We then report the change in the measurement on the sample for which top influences were removed.

Monte Carlo sampling of measurements Again, for all computations of the “true” measurement functions for the original and the retrained models used for calculating the difference in loss after retraining we use $5000$ samples to estimate the measurement.

### G.7 Training details

For CIFAR-10 and CIFAR-2 we again follow the training procedure outlined in Ho et al. (), with the only difference being a shortened number of training iterations.
For CIFAR-10, we train for $160000$ steps (compared to $800000$ in Ho et al. ()) for the full model, and $80000$ steps for the subsampled datasets (410 epochs in each case).
On CIFAR-2, we train for $32000$ steps for the model trained on the full dataset, and $16000$ steps for the subsampled datasets ($\leavevmode\nobreak\ 800$ epochs).
We train for significantly longer than Zheng et al. (), as we noticed the models trained using their procedure were somewhat significantly undertrained (some per-diffusion-timestep training losses $\ell_{t}(\theta,x)$ have not converged).
We also use a cosine learning-rate schedule for the CIFAR-2 models.

### G.8 Handling of data augmentations

In the presentation in Section 2, we ignore for the sake of clear presentation the reality that in most diffusion modelling applications we also apply data augmentations to the data.
For example, the training loss $\mathcal{L}_{\mathcal{D}}$ in Equation 3 in practice often takes the form:

$$ $\mathcal{L}_{\mathcal{D}}=\frac{1}{N}\sum_{n=1}^{N}\mathbb{E}_{\tilde{x}_{n}} \left[\ell(\theta,\tilde{x}_{n})\right],$ $$

where $\tilde{x}_{n}$ is the data point $x_{n}$ after applying a (random) data augmentation to it.
This needs to be taken into account 1) when defining the GGN, as the expectation over the data augmentations $\mathbb{E}_{\tilde{x}_{n}}$ can either be considered as part of the outer expectation $\mathbb{E}_{z}$, or as part of the loss $\rho$ (see Section 2.2.1), 2) when computing the per-example train loss gradients for influence functions, 3) when computing the loss measurement function.

When computing $\mathrm{GGN}^{\texttt{model}}_{\mathcal{D}}$ in Equation 9, we treat data augmentations as being part of the out “empirical data distribution”. In other words, we would simply replace the expectation $\mathbb{E}_{x_{n}}$ in the definition of the GGN with a nested expectation $\mathbb{E}_{x_{n}}\mathbb{E}_{\tilde{x}_{n}}$:

$$ $\mathrm{GGN}_{\mathcal{D}}^{\texttt{model}}(\theta)={\mathbb{E}_{x_{n}}\left[{ \color[rgb]{0.89453125,0.7734375,0.46484375}\definecolor[named]{pgfstrokecolor }{rgb}{0.89453125,0.7734375,0.46484375}\mathbb{E}_{\tilde{x}_{n}}}\left[ \mathbb{E}_{\tilde{t}}\left[\mathbb{E}_{{\color[rgb]{ 0.89453125,0.7734375,0.46484375}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.89453125,0.7734375,0.46484375}x^{(\tilde{t})}},{\epsilon}^{(\tilde{t})}} \left[{\nabla_{\theta}^{\mathchoice{\raise 1.3611pt\hbox{$\displaystyle \intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0.95276pt \hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$\scriptscriptstyle \intercal$}}}{\epsilon}_{\theta}^{\tilde{t}}\left({\color[rgb]{ 0.89453125,0.7734375,0.46484375}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.89453125,0.7734375,0.46484375}{x}^{(\tilde{t})}}\right)}{\left(2I\right)}{ \nabla_{\theta}{\epsilon}_{\theta}^{\tilde{t}}\left({\color[rgb]{ 0.89453125,0.7734375,0.46484375}\definecolor[named]{pgfstrokecolor}{rgb}{ 0.89453125,0.7734375,0.46484375}{x}^{(\tilde{t})}}\right)}\right]\right]\right ]\right]}.$ $$

with ${\color[rgb]{0.89453125,0.7734375,0.46484375}\definecolor[named]{
pgfstrokecolor}{rgb}{0.89453125,0.7734375,0.46484375}x^{(\tilde{t})}}$ now being sampled from the diffusion process $q(x^{(\tilde{t})}|\tilde{x}_{n})$ conditioned on the augmented sample $\tilde{x}_{n}$.
The terms changing from the original equation are indicated in yellow.
The “Fisher” expression amenable to MC sampling takes the form:

$$ $\displaystyle\mathrm{F}_{\mathcal{D}}(\theta)=\mathbb{E}_{x_{n}}\left[{\mathbb {E}_{\tilde{x}_{n}}}\left[\mathbb{E}_{\tilde{t}}\left[\mathbb{E}_{x_{n}^{( \tilde{t})},{\epsilon}^{(\tilde{t})}}\mathbb{E}_{{\epsilon}_{\texttt{mod}}} \left[g_{n}(\theta)g_{n}(\theta)^{\mathchoice{\raise 1.3611pt\hbox{$ \displaystyle\intercal$}}{\raise 1.3611pt\hbox{$\textstyle\intercal$}}{\raise 0 .95276pt\hbox{$\scriptstyle\intercal$}}{\raise 0.68054pt\hbox{$ \scriptscriptstyle\intercal$}}}\right]\right]\right]\right],$ $\displaystyle{\epsilon}_{\texttt{mod}}\sim{\mathcal{N}}\left({\epsilon}_{ \theta}^{\tilde{t}}\left({{x}_{n}^{(\tilde{t})}}\right),I\right),$ $$

where, again, $g_{n}(\theta)=\nabla_{\theta}{\lVert{\epsilon}_{\texttt{mod}}-{\epsilon}_{
\theta}^{\tilde{t}}(x_{n}^{(\tilde{t})})\rVert^{2}}$.

When computing $\mathrm{GGN}^{\texttt{loss}}_{\mathcal{D}}$ in Equation 11, however, we treat the expectation over daea augmentations as being part of the loss $\rho$, in order to be more compatible with the implementations of TRAK (Park et al., ) in prior works that rely on an empirical GGN (Zheng et al., ; Georgiev et al., ).(^10^1010The implementations of these methods store the (randomly projected) per-example training loss gradients for each example before computing the Hessian approximation. Hence, unless data augmentation is considered to be part of the per-example training loss, the number of gradients to be stored would be increased by the number of data augmentation samples taken.)Hence, the GGN in Equation 11 takes the form:

$$ $\displaystyle\mathrm{GGN}_{\mathcal{D}}^{\texttt{loss}}(\theta)$ $\displaystyle=\mathbb{E}_{x_{n}}\left[\nabla_{\theta}\left(\mathbb{E}_{\tilde{ x}_{n}}\left[\ell(\theta,\tilde{x}_{n})\right]\right)\nabla^{\intercal}_{ \theta}\underbrace{\left(\mathbb{E}_{\tilde{x}_{n}}\left[\ell(\theta,\tilde{x} _{n})\right]\right)}_{\tilde{\ell}(\theta,x_{n})}\right]$ $\displaystyle=\mathbb{E}_{x_{n}}\left[\nabla_{\theta}\tilde{\ell}(\theta, \tilde{x}_{n})\nabla^{\intercal}_{\theta}\tilde{\ell}(\theta,\tilde{x}_{n}) \right],$ $$

where $\tilde{\ell}$ is the per-example loss in expectation over data-augmentations.
This is how the Hessian approximation is computed both when we’re using K-FAC with $\mathrm{GGN}^{\texttt{model}}_{\mathcal{D}}$ in presence of data augmentations, or when we’re using random projections (TRAK and D-TRAK).

When computing the training loss gradient in influence function approximation in equation Equation 5, we again simply replace the per-example training loss $\ell(\theta^{\star},x_{j})$ with the per-example training loss averaged over data augmentations $\tilde{\ell}(\theta^{\star},x_{j})$, so that the training loss ${\mathcal{L}}_{\mathcal{D}}$ can still be written as a finite sum of per-example losses as required for the derivation of influence functions.

For the measurement function $m$ in Equation 6, we assume we are interested in the log probability of (or loss on) a particular query example in the particular variation in which it has appeared, so we do not take data augmentations into account in the measurement function.

Lastly, since computing the training loss gradients for the influence function approximation for diffusion models usually requires drawing MC samples anyways (e.g. averaging per-diffusion timestep losses over the diffusion times $\tilde{t}$ and noise samples $\epsilon^{(t)}$), we simply report the total number of MC samples per data point, where data augmentations, diffusion time $\tilde{t}$, etc. are all drawn independently for each sample.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 2
questions: - In practice, how would we tune the damping parameter for K-FAC Influence? How would the tuning procedure differ from that for D-TRAK?
- Observation 2 shows that F-FAC Influence overestimates the impact of data removal, whereas Koh et al. (2019) show that influence functions underestimates the impact of data removal. What are possible explanations for these differing observations?
- Observation 3 seems to have been shown in Kadkhodaie et al. (2024; Figure 2 in particular). What additional insights are provided? For example, does Observation 3 only hold for ELBO, or does it also hold for the simple diffusion loss?


References
- Koh et al. (2019) - On the Accuracy of Influence Functions for Measuring Group Effects.
- Kadkhodaie et al. (2024) - Generalization in Diffusion Models Arises from Geometry-Adaptive Harmonic Representations.
soundness: 3
strengths: - A theoretical framework is shown to unify influence functions, Journey-TRAK, and D-TRAK for diffusion models. This provides some clarity to the field of data attribution for diffusion models, which currently seems rather empirical.
- The theoretical framework motivates the design choices for approximating influence functions, which actually lead to better performance (as shown in Figure 4).
- Empirical observations (Section 4.1) are made to better understand how to apply influence functions to diffusion models.
summary: This paper extends influence functions for data attribution of diffusion models. With an exposition focusing on approximating a Hessian with a generalized Gauss-Newton matrix (GCN), this paper justifies specific design choices (linearizing the neural networks combined with K-FAC-expand approximation) that lead to better performing influence functions. The theoretical framework also subsumes previous work such as Journey-TRAK and D-TRAK. Empirically, influence functions perform similarly to D-TRAK on CIFAR-2 and CIFAR-10, with lower sensitivity to the choice of the damping parameter. Finally, empirical observations are made to show the challenges when applying influence functions to diffusion models.
weaknesses: - The LDS and counterfactual evaluations are limited to DDPM on CIFAR. The paper would be strengthened if K-FAC Influence is also evaluated on LDM and LoRA-fine-tuned Stable Diffusion (as in the Journey-TRAK and D-TRAK papers). Furthermore, since the observations in Section 4.1 are empirical, not theoretical, they need to be validated in other model architectures and datasets to make the claims more general.
- The proposed method K-FAC Influence is also sensitive to the choice of the damping parameter (LDS can range from ~15% to ~5% on CIFAR-10). This range of variation seems not much different from D-TRAK.
- Observation 3 might be sensitive to the amount of training data removed. For example, if 90% of the data are removed, there could be more variation in the retrained diffusion models, as suggested by Kadkhodaie et al. (2024; Figure 2 in particular). In this sense, the claim in Observation 3 might need to be adjusted to be less general.
- This is minor: It would be useful to include a brief summary of Grosse and Martens (e.g., the reduce vs. expand settings), perhaps in the Appendix.

References
 - Kadkhodaie et al. (2024) - Generalization in Diffusion Models Arises from Geometry-Adaptive Harmonic Representations.

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: You explicitly mention that K-FAC is defined for linear and convolutional layers (page 6), but this doesn’t seem to include attention layers, which all diffusion architectures that I’m aware of use. Is K-FAC also defined for attention in prior work, or could you define it?

Could you provide a runtime analysis, or some idea of what the runtime is (e.g. how much time would it take to compute a all influence values for a single sample on an A100)? How does this runtime (and possibly space requirement) scale with the complexity of the model? How much computational resources would I need if I wanted to run this on my favorite existing diffusion model?
soundness: 3
strengths: This paper combines interesting ideas to tractably approximate influence functions for diffusion models, and achieve state-of-the-art performance on training data attribution on diffusion models, which is a challenging problem.
summary: In this work the authors adapt influence functions to work with diffusion models. Computing influence functions typically requires one to compute the inverse of the Hessian, which is intractable for deep learning models. Therefore the authors make use of a series of approximations. Firstly, they approximate the Hessian with the generalized Gauss-Newton matrix (GGN). Then they approximate the GGN with a block diagonal matrix, where the diagonal blocks are defined by layers. Finally, they approximate each block as a kronecker product.

The authors then evaluate their method on “linear data modelling score” and “retraining without top influencers”, and find their methods achieve state-of-the-art on both metrics.
weaknesses: While it is clear that the proposed method is much more scalable than computing the raw influence function, it is difficult to assess how scalable this method actually is based on the main body of this paper (e.g. things like what is the runtime, how does it scale with the model size, what are the space requirements). A reader ought to have a good idea of how feasible this would be to implement for themselves.

It is not clear whether the proposed method outperforms the existing methods given the overlaps in standard error. Perhaps statistical significance could be reported here to give a better idea?

The experimental evaluation is only performed on CIFARs 2 and 10, so it is unclear how well this method scales to larger image generators, both in the quality of attributions and in terms of compute.

---

confidence: 4
contribution: 4
initial_rating: 10
presentation: 3
questions: In the last paragraph before the discussion (starting from line (507), the authors argue that the ELBO could be a poor proxy of the marginal probability due to the fact that it is rather insensitive to the resampling of up 50% of the dataset, meaning that such a resampling do not meaningfully change the ELBO of the generated images.

I think that this conclusion is incorrect. There is in fact evidence that diffusion training enters a 'generalization phase' after a finite dataset size, after which the generative probabilities remain approximately constant even under a non-overlapping split of the training set. This for example can be seen in (Zahra, 2023).

References:
Kadkhodaie, Zahra, et al. "Generalization in diffusion models arises from geometry-adaptive harmonic representation." arXiv preprint arXiv:2310.02557 (2023).
soundness: 3
strengths: This work is very valuable for anyone interested in IF estimation in generative diffusion models. This is likely going to be a very important tool due to privacy and copyright concerns associated with image and video generation. 

Both the theoretical discussion and the experimental analysis is systematic and well executed. The paper is very easy to follow and it manages to both provide the intuition and the theoretical justification for its design choices.

I particularly appreciated the discussion of the results, which highlight both the strengths and weaknesses of each approach and also the limitation of the approximate IF estimation when compared with exact retraining.
summary: The paper provides a systematic treatment of tractable influence function (IF) estimation for generative diffusion models. Several forms of IFs are considered: IF of the ELBO, IF of the loss, IF of the path probabilities. These approaches are compared both conceptually and experimentally. The IF of the marginal probabilities of the images is discussed but not implemented.
A large fraction of the paper deals with the intractability of the Jacobian estimation involved in the calculation of the IF and offers and compares several approximate solutions.

The paper provides a thorough experimental analysis of the different IF approaches on several datasets and compares the results with the exact retraining method.

[edit] I increased my score due to the new analysis and results provided during the rebuttal phase, which answered some of the opened questions and provided links with the recent literature.
weaknesses: Arguably, the marginal likelihood logp(x) is the ideal target for a IFa analysis. While it is somewhat challenging to compute using the deterministic integration formula with stochastic or deterministic  trace estimation, it is not intractable and in my opinion it should have been analyzed. For example, having the analysis of the IF of logp(x) would help to make sense of the counterintuitive behavior of the IF of the ELBO.

It would have also appreciated to see an analysis on more tractable low-dimensional models with small networks, where the Hessian can be computed exactly. This analysis can be used to evaluate the relative performance of the Hessian approximation techniques, at least in this simple regime.

Meta-review: This paper extends influence functions to diffusion models, addressing the computational challenges of Hessian inversions with scalable approximations, including Kronecker-factored Approximate Curvature (K-FAC) and its eigenvalue-corrected variant (EK-FAC). The proposed method demonstrates state-of-the-art performance on benchmarks like Linear Data Modeling Score (LDS). The authors provide theoretical insights that unify influence functions with prior work (e.g., TRAK and D-TRAK).

Decision: Accept (Oral)