# A Bayesian Perspective On Training Data Attribution

Elisa Nguyen
Tübingen AI Center
University of Tübingen
&Minjoon Seo
KAIST AI
&Seong Joon Oh
Tübingen AI Center
University of Tübingen

###### Abstract

Training data attribution (TDA) techniques find influential training data for
the model’s prediction on the test data of interest. They approximate the
impact of down- or up-weighting a particular training sample. While
conceptually useful, they are hardly applicable in practice, particularly
because of their sensitivity to different model initialisation. In this paper,
we introduce a Bayesian perspective on the TDA task, where the learned model
is treated as a Bayesian posterior and the TDA estimates as random variables.
From this novel viewpoint, we observe that the influence of an individual
training sample is often overshadowed by the noise stemming from model
initialisation and SGD batch composition. Based on this observation, we argue
that TDA can only be reliably used for explaining model predictions that are
consistently influenced by certain training data, independent of other noise
factors. Our experiments demonstrate the rarity of such noise-independent
training-test data pairs but confirm their existence. We recommend that future
researchers and practitioners trust TDA estimates only in such cases. Further,
we find a disagreement between ground truth and estimated TDA distributions
and encourage future work to study this gap. Code is provided at
https://github.com/ElisaNguyen/bayesian-tda.

##  1 Introduction

Understanding how machine learning models arrive at decisions is desirable for
social, legal and ethical reasons, particularly for opaque deep learning
models [1]. One approach to explanations is the data-centric approach of
training data attribution (TDA). As the name suggests, TDA finds attributing
training samples for a model decision, uncovering which part of the training
data is relevant. The attribution $\tau$ of a training sample $z_{j}$ on
another sample $z$ is usually defined as the change of model loss
$\mathcal{L}$ on $z$ when the model is retrained without $z_{j}$ [2]:

$\displaystyle\tau(z_{j},z):=\mathcal{L}(z;\theta_{\setminus
j})-\mathcal{L}(z;\theta)$ (1)

where $\theta$ is a model trained on the entire dataset $\mathcal{D}$ and
$\theta_{\setminus j}$ is a model trained on the same set without $z_{j}$.
Since the direct computation of Equation 1 is expensive, various TDA
techniques for approximating the quantity have been proposed, such as
influence functions [3] and TracIn [4]. Their approximations are often based
on some form of inner product between the parameter gradients
$\nabla_{\theta}\mathcal{L}(z;\theta)$ and
$\nabla_{\theta}\mathcal{L}(z_{j};\theta)$.

Knowing how training samples attribute to a model decision provides an
actionable understanding of the training data distribution, especially in
cases of model error. TDA methods can identify the training samples that are
most relevant to an error and therefore enable users to understand why the
error occurred (e.g. due to domain mismatch of test and training data or
wrongly labelled training data) [3]. Additionally, TDA gives them the tool to
address the errors by e.g. changing the model directly through the training
data. Even in non-erroneous cases, understanding the attributing training data
may enable users affected by model decisions to contest the decisions if the
attributing training data is noisy or of low quality [5].

At the same time, TDA methods, especially influence functions [3], have been
criticised for their fragility when applied to deep models [6, 7, 8]. The main
reasons are model complexity and the stochasticity of deep model training.
While the former poses a challenge specifically for influence functions as
they rely on strong convexity assumptions, the latter is a more general
challenge [6, 9]. The randomness inherent to the training process does not
only lead to variation in the learned model parameters but also in TDA scores,
which makes them untrustworthy. Hence, K & Søgaard [9] recommend using
expected TDA scores for increased stability.

We argue that solely considering the expectation is not sufficient to ensure
the reliability of TDA but requires inspecting the variance, too. We introduce
a Bayesian perspective on the TDA task, noticing that there is no
deterministic mapping from a dataset $\mathcal{D}$ to the corresponding model
$\theta$ for deep neural networks. The learned model depends on the
initialisation and batch composition in the stochastic gradient descent (SGD)
optimiser. We capture the resulting randomness via Bayesian model posterior
$p(\theta|\mathcal{D})$ over the parameter space [10, 11, 12]. In turn, the
TDA estimate (Equation 1) is a random variable that depends on two posteriors,
$p(\theta|\mathcal{D})$ and $p(\theta_{\setminus j}|\mathcal{D}_{\setminus
j})$.

This viewpoint leads to a few insights into the practical usage and evaluation
of TDA techniques. We confirm quantitatively that the ground-truth influence
$\tau(z_{j},z)$ is often dominated by the noise:
$\sqrt{\text{Var}(\tau)}>\mathbb{E}|\tau|$. We argue that it is practically
difficult to apply any TDA technique on pairs $(z_{j},z)$ whose ground-truth
attributions $\tau(z_{j},z)$ are noisy in the first place. Likewise, any
evaluation of TDA methods on such high-variance pairs would not be reliable.

Nonetheless, we are optimistic that TDA techniques are useful in practice,
particularly for train-test pairs with high signal-to-noise ratios:
$\sqrt{\text{Var}(\tau)}\ll\mathbb{E}|\tau|$. We observe that such pairs are
rare but consistently present in multiple experiments. We recommend that
researchers and practitioners confine their usage to scenarios where the
signal-to-noise ratios are expected to be large enough.

Our contributions are as follows: (1) Bayesian formulation of the training
data attribution (TDA) task. (2) Observation that the ground-truth TDA values
are often unreliable and highly variable. (3) Recommendation for the community
to use the TDA tools only when the expected noise level is low. (4)
Experimental analysis of the contributing factors to the variance of ground-
truth TDA values. (5) Observation that the estimated TDA is closer to a local
perturbation in the model, while ground-truth TDA is more global.

##  2 Background

We cover the background materials for the paper, including the concept,
method, and evaluation of training data attribution (TDA) methods and Bayesian
deep learning.

###  2.1 Training data attribution (TDA)

We introduce the TDA task, a few representative TDA methods, and existing
evaluation strategies.

#### TDA task.

Given a model $f_{\theta}$ parametrised by $\theta$, a training set
$\mathcal{D}:=\\{z_{1},\cdots,z_{N}\\}$, and a test sample $z$, one is
interested in the impact of a training sample $z_{j}$ on the model’s behaviour
on the test sample $z$. In the TDA context, one is often interested in the
counterfactual change in the loss value for $z$ after leave-one-out (LOO)
training, when $z_{j}$ is excluded from the training set (Equation 1). TDA has
been considered in different use cases, such as understanding the bias in word
embeddings [13], fact tracing in language model outputs [14] and measuring the
robustness of model predictions [5].

#### TDA methods.

The conceptually most straightforward way to compute the difference due to LOO
training (Equation 1) is to compute it directly. However, this is
computationally expensive, as it involves the learning algorithm for obtaining
$\theta_{\setminus j}$ for every $j$. This gives rise to various TDA
techniques that find _approximate_ estimates $\tau^{\prime}(z_{j},z)$ of LOO.

A prominent example of such approximation is the influence function (IF)
method [3] by Koh & Liang. Under strong smoothness assumptions, they have
approximated Equation 1 by:

$\tau^{\prime}(z_{j},z):=-\nabla_{\theta}\mathcal{L}(z;\theta)^{\top}H_{\theta}^{-1}\nabla_{\theta}\mathcal{L}(z_{j};\theta)$
(2)

where $\nabla_{\theta}\mathcal{L}(z;\theta)$ and
$\nabla_{\theta}\mathcal{L}(z_{j};\theta)$ refer to the parameter gradients of
$f_{\theta}$ for $z$ and $z_{j}$ respectively. Recognising the difficulty of
scaling up the inverse Hessian computation $H_{\theta}^{-1}$ and the high
dimensionality of operations in Equation 2, subsequent papers have proposed
further approximations to speed up the computation [15, 16]. Charpiat et al.
[17] have analysed the influence of $z_{j}$ on $z$ by dropping the need to
compute the Hessian and formulating influence as the loss change when an
additional training step (ATS) on $z_{j}$ is taken:

$\tau(z_{j},z):=\mathcal{L}(z;\theta_{+j})-\mathcal{L}(z;\theta)$ (3)

where $\theta_{+j}$ is a learned model parameter with $\mathcal{D}$ and an
additional step on $z_{j}$. They propose two approximations:

Grad-Dot (GD):
$\displaystyle\quad\tau^{\prime}(z_{j},z):=\nabla_{\theta}\mathcal{L}(z_{j};\theta)^{\top}\nabla_{\theta}\mathcal{L}(z;\theta)$
(4) Grad-Cos (GC):
$\displaystyle\quad\tau^{\prime}(z_{j},z):=\frac{\nabla_{\theta}\mathcal{L}(z_{j};\theta)}{\|\nabla_{\theta}\mathcal{L}(z_{j};\theta)\|}^{\top}\frac{\nabla_{\theta}\mathcal{L}(z;\theta)}{\|\nabla_{\theta}\mathcal{L}(z;\theta)\|}$
(5)

This method is closely linked to TracIn [4] which computes the Grad-Dot not
just at the end of the training, but averages the regular Grad-Dot
similarities throughout the model training iterations. We note later in our
analysis that within our Bayesian treatment of TDA, the TracIn method
coincides conceptually with the Grad-Dot method. In our analysis, we study the
sensitivity of LOO and the above TDA methods against noise.

#### TDA evaluation.

The primal aim of TDA methods is to measure how well they approximate the
ground-truth LOO values. This is often done by measuring the correlation
between the estimates from each TDA method and the ground-truth LOO values
(Equation 1) [3, 6, 7, 9]. They use either a linear (Pearson) correlation or a
rank (Spearman) correlation over a small number of train-test sample pairs
$(z_{j},z)$ due to the computational burden of computing the actual LOO
values, especially for larger models. Usually, a few samples $z$ are chosen
for a comparison against LOO, e.g. Koh & Liang [3] report results for one $z$
and Guo et al. [15] for 10 samples $z$. Others have adopted indirect
evaluation metrics such as the retrieval performance of mislabelled or
poisoned training data based on the TDA estimates [3, 18, 15, 9, 16]. In this
work, we adopt the Pearson and Spearman correlation metrics and discuss ways
to extend them when the target (LOO) and estimates (TDA) are both random
variables.

###  2.2 Bayesian deep learning.

Bayesian machine learning treats the learned model as a posterior distribution
over the parameter space, rather than a single point:

$\displaystyle
p(\theta|\mathcal{D})=p(\mathcal{D}|\theta)p(\theta)/p(\mathcal{D}).$ (6)

Bayesian ML nicely captures the intuition that the mapping from a training set
$\mathcal{D}$ to the learned model $p(\theta|\mathcal{D})$ is not a
deterministic mapping, especially for non-convex models like deep neural
networks (DNNs). Depending on the initialisation, among other factors, DNN
training almost always ends up learning vastly different parameters.

The estimation of the true posterior is indeed difficult for complex models
like DNNs. The field of Bayesian deep learning is dedicated to the
interpretation of certain random elements in DNN training as sources of
randomness for the approximated Bayesian posteriors. For example, if Dropout
[19] is used for training a model, it may be used at test time to let users
sample $\theta$ from the posterior distribution $p(\theta|\mathcal{D})$ [20].
More generally-used components like stochastic gradient descent (SGD) have
also been interpreted as sources of randomness. The random walk induced by SGD
iterations in the parameter space can be viewed as a Markov Chain Monte-Carlo
sampler from the posterior distribution, after a slight modification of the
optimisation algorithm (Stochastic Gradient Langevin Dynamics [10]).
Similarly, the last few iterations of the vanilla SGD iterations may also be
treated as samples from the posterior, resulting in more widely applicable
Bayesian methods like Stochastic Weight Averaging (SWA) [12, 21]. Finally, the
random initialisation of DNNs has also been exploited for modelling posterior
randomness; training multiple versions of the same model with different
initial parameters may be interpreted as samples from the posterior [11]. We
show in the next section how the Bayesian viewpoint will help us model the
sources of stochasticity for TDA estimates.

Figure 1: A Bayesian interpretation of training data attribution (TDA).

##  3 A Bayesian perspective on training data attribution

Training data attribution (TDA) $\tau(z_{j},z)$ is defined as the attribution
of one training sample $z_{j}$ to another sample $z$ in terms of how a target
metric like the loss of a sample $\mathcal{L}(z,\theta)$ changes when the
model is trained without $z_{j}$ (Equation 1). We note here that according to
the definition, we are interested in the impact of the change in the dataset
from $\mathcal{D}$ to $\mathcal{D}_{\setminus j}$, rather than the change in
the model parameter. From a Bayesian perspective, a change in the training
dataset leads to a shift in the posterior distribution,
$p(\theta|\mathcal{D})\rightarrow p(\theta_{\setminus
j}|\mathcal{D}_{\setminus j})$, leading to the definition of TDA as a random
variable:

$\displaystyle\tau(z_{j},z)$
$\displaystyle:=\mathcal{L}(z;\mathcal{D}_{\setminus
j})-\mathcal{L}(z;\mathcal{D})=\mathcal{L}(z;\theta_{\setminus
j})-\mathcal{L}(z;\theta)$ (7)

where $\theta\sim p(\theta|\mathcal{D})$ and $\theta_{\setminus j}\sim
p(\theta_{\setminus j}|\mathcal{D}_{\setminus j})$. This interpretation is
more natural, given the non-uniqueness of the mapping from a training dataset
$\mathcal{D}$ to the optimal model parameter $\theta$ for general, non-convex
models like DNNs.

#### Sampling TDA values.

One could plug in various Bayesian DL techniques (§2.2) to compute samples of
$p(\theta|\mathcal{D})$, which can be used to get the samples of
$\tau(z_{j},z)$. In our work, we use the Stochastic Weight Averaging (SWA)
[12, 21] and Deep Ensemble (DE) [11] which are applicable to a wide class of
deep models. More specifically, we obtain $T$ samples
$\theta^{(1)},\cdots,\theta^{(T)}\sim p(\theta|\mathcal{D})$ either by taking
the last $T$ model checkpoints of the SGD iterations (SWA) or by taking the
last model checkpoints from $T$ different model initialisations (DE). The same
is done for the counterfactual posterior $\theta^{(1)}_{\setminus
j},\cdots,\theta^{(T)}_{\setminus j}\sim p(\theta_{\setminus
j}|\mathcal{D}_{\setminus j})$.

#### Statistical analysis on TDA.

The simplest statistics for the TDA $\tau(z_{j},z)$ are the mean and variance:

$\displaystyle\mathbb{E}[\tau(z_{j},z)]$
$\displaystyle=\frac{1}{T}\sum_{t}\mathcal{L}(z;\theta^{(t)}_{\setminus
j})-\mathcal{L}(z;\theta^{(t)})$ (8) $\displaystyle\text{Var}[\tau(z_{j},z)]$
$\displaystyle=\frac{1}{T^{2}}\sum_{t,t^{\prime}}\left(\mathcal{L}(z;\theta^{(t)}_{\setminus
j})-\mathcal{L}(z;\theta^{(t^{\prime})})-\mathbb{E}[\tau(z_{j},z)]\right)^{2}$
(9)

Our main interest lies in whether the influence of the training data $z_{j}$
on the test data $z$ is statistically significant and not dominated by the
inherent noise. For this purpose, we design a Student t-test for quantifying
the statistical significance. Our null and alternative hypotheses are:

$\displaystyle H_{0}:\tau(z_{j},z)=0\quad\quad H_{1}:\tau(z_{j},z)\neq 0.$
(10)

We consider the test statistic based on sample mean and variance:

$\displaystyle
t=\frac{{\tau}(z_{j},z)-\mathbb{E}[\tau(z_{j},z)]}{\sqrt{\text{Vars}[\tau(z_{j},z)]/T^{2}}}.$
(11)

Vars refers to the sample variance where the denominator in Equation 9 is
$T^{2}-1$ instead. We report the significance of the absolute TDA
$|\tau(z_{j},z)|$ for every train-test pair $(z_{j},z)$ by computing the
p-value corresponding to the t-test statistic. The greater the p-value is, the
greater the dominance of noise is for the TDA estimate.

#### TDA methods likewise estimate random quantities.

Approximate TDA methods like influence functions (IF), Grad-Dot, and Grad-Cos
(§2.1) also predict random quantities $\tau^{\prime}(z_{j},z)$. For example,
IF predicts
$\tau(z_{j},z)\approx\tau^{\prime}(z_{j},z):=-\nabla_{\theta}\mathcal{L}(z_{j};\theta)^{\top}H_{\theta}^{-1}\nabla_{\theta}\mathcal{L}(z;\theta)$,
where one may sample $\theta$ from the posterior $p(\theta|\mathcal{D})$.
Similar statistical analyses on the TDA estimations can be performed as above,
including the mean and variance computations and statistical testing for the
significance of influence.

#### Evaluating TDA as a random variable.

Previously, the LOO-based TDA values $\tau(z_{j},z)$ and the estimates from
various approximate TDA methods $\tau^{\prime}(z_{j},z)$ are compared via
correlation measures like Pearson or Spearman. Our treatment of those
quantities as 1-D random variables poses a novel challenge for evaluation
because there exists no inborn notion of ordering among 1-D random variables.
We address the challenge by examining the approximation ability of TDA methods
for both the first and second moments of the true TDA values $\tau(z_{j},z)$.
More specifically, we compute the Pearson and Spearman correlation for both
the mean (Equation 8) and variance (Equation 9) between the ground-truth
$\tau(z_{j},z)$ and estimated TDA $\tau^{\prime}(z_{j},z)$ values across
multiple train-test pairs $(z_{j},z)$.

##  4 Experiments

We introduce our experimental settings, present analyses on factors
contributing to the reliability of TDA values, compare TDA methods, and draw
suggestions on the evaluation practice of TDA.

Figure 2: Sources of randomness for Bayesian posteriors. In each case, the
training starts from initialisation $\theta_{0}$. Depending on whether $z_{j}$
is included in the training data, one has either samples from the original
posterior $p(\theta|\mathcal{D})$ or from the counterfactual posterior
$p(\theta_{\setminus j}|\mathcal{D}_{\setminus j})$. For deep ensemble [11],
the randomness stems either from random initialisation (DE-Init) or from SGD
batch composition (DE-Batch). For stochastic weight averaging (SWA) [12, 21],
last few checkpoints of the training are treated as posterior samples.

###  4.1 Implementation details

We illustrate the specific details of our implementation. See the Appendix for
further information.

#### TDA methods.

We study different TDA methods from a Bayesian perspective. We test the
methods introduced in §2.1 for estimating TDA: influence functions (IF) [3],
Grad-Dot (GD) and Grad-Cos (GC) [17]. We use the PyTorch implementation of IF
from Guo et al. [15] and modify it for our models. As the ground-truth target,
we consider Leave-one-out training (LOO) [3]. For LOO, we remove of $z_{j}$
from the training set $\mathcal{D}$ by zeroing out the weight for sample
$z_{j}$ towards the loss. Additionally, we include Charpiat et al.’s [17]
notion of TDA that a training data point $z_{j}$ attributes more if an
additional training step (ATS) on it changes the test loss more significantly.

#### Inducing randomness in posterior $p(\theta|\mathcal{D})$.

In §2.2, we have introduced the interpretation of various elements around
model training as sources of randomness for Bayesian posterior. We summarise
our methods for inducing randomness in Figure 2. We use the notion of the Deep
Ensemble (DE) [11] to sample from the posterior. In a variant of DE with the
initialisation as the source of randomness (DE-Init), we train each of
$T_{\text{DE}}$ randomly initialised parameters $\theta^{(t)}_{0}$ on either
$\mathcal{D}$ or $\mathcal{D}_{\setminus j}$. The resulting parameter sets,
$\theta^{(t)}$ and $\theta^{(t)}_{\setminus j}$, are treated as samples from
respective posteriors. We also consider the batch composition in stochastic
gradient descent (SGD) as the source of randomness (DE-Batch). In this case,
we train from one initial parameter $\theta_{0}$ with $T_{\text{DE}}$
different random shuffles $\pi^{(t)}$ of the training sets $\mathcal{D}$ and
$\mathcal{D}_{\setminus j}$. This results in two sets of samples from the
original and counterfactual posteriors. We increase the number of samples by
taking the last $T_{\text{SWA}}$ checkpoints as the Stochastic Weight
Averaging (SWA) samples [12, 21]. For Grad-Dot, this coincides with the
definition of TracIn [4] as we average the dot products across checkpoints. In
total, we take $T=T_{\text{DE}}\times T_{\text{SWA}}=10\times 5$ samples from
$p(\theta|\mathcal{D})$ and $p(\theta_{\setminus j}|\mathcal{D}_{\setminus
j})$ to estimate $\tau(z_{j},z)$.

#### Datasets $\mathcal{D}$.

To enable an exhaustive analysis of every train-test pair $(z_{j},z)$, we
define smaller datasets. We use variants of MNIST [22] limited to three
classes (MNIST3), and CIFAR10 [23]. For MNIST3, we sample a training set of
size 150 and a test set of size 900, i.e. 135,000 train-test pairs. For
CIFAR10, we define the training and test set at size 500, i.e. 250,000 train-
test pairs.

#### Models.

We consider two types of image classifiers, convolutional neural networks
(CNN, [24]) and visual transformers (ViT, [25]). For ViT variants, instead of
full finetuning, we use LoRA adapter layers [26] to minimise the number of
parameters being tuned. The number of trainable parameters of ViT+LoRA
(597,514) is comparable to adding another layer to our CNN (225,034).

###  4.2 Reliability of TDA evaluation

We assess the reliability of TDA evaluation by measuring the degrees of noise
in both the ground-truth TDA (LOO) $\tau(z_{j},z)$ and the estimated TDA
$\tau^{\prime}(z_{j},z)$. The noise level is measured with the p-value of the
Student-t hypothesis testing to determine if the absolute TDA values are
significantly greater than the sample noise (§3).

Table 1: Stability of TDA estimates. We report p-values for the ground-truth TDA $\tau(z_{j},z)$ (LOO) and the estimated TDA values $\tau^{\prime}(z_{j},z)$ (rest 4 columns). The p-values are averaged across all train-test pairs $(z_{j},z)$. We use the CNN model throughout. Data | Randomness | LOO | ATS | IF | GD | GC
---|---|---|---|---|---|---
MNIST3 |  SWA$+$DE-Init | 0.331 | 0.254 | 0.352 | 0.363 | 0.003
SWA$+$DE-Batch | 0.025 | 0.039 | 0.000 | 0.000 | 0.000
CIFAR10 |  SWA$+$DE-Init | 0.692 | 0.437 | 0.575 | 0.587 | 0.356
SWA$+$DE-Batch | 0.487 | 0.296 | 0.484 | 0.517 | 0.236
MNIST3, SWA+DE-I | MNIST3, SWA+DE-B | CIFAR10, SWA+DE-I | CIFAR10, SWA+DE-B
---|---|---|---
Figure 3: Stability of TDA estimates per train-test pair. Distribution of
p-values for ground-truth TDA (LOO) for different experiments.

We report the results in Table 1. Generally, we observe many TDA measurements,
ground-truth and estimations likewise, are unstable with non-significant
p-values ($>0.05$). In particular, even the ground-truth LOO shows p-values of
0.331 on MNIST3 and 0.692 for CIFAR10 (SWA+DE-Init). In these cases, the noise
effectively dominates the signal and any evaluation that does not consider the
variance in the posterior $p(\theta|\mathcal{D})$ is likely to be misleading.
This confirms the reports in [9] that TDA values are sensitive to model
initialisation.

TDA methods often show similar levels of instability. For example, the IF
attains p-values 0.352 and 0.575 on MNIST3 and CIFAR10, respectively, roughly
matching the LOO case. Grad-Cos is an exception: it attains lower p-values
than the other TDA methods (0.003 and 0.356 for MNIST3 and CIFAR10,
respectively). We interpret this as an overconfident TDA estimation.
Practitioners shall be wary of using TDA methods that are unreasonably stable
when the ground-truth TDA itself is not.

###  4.3 Factors influencing the variability of TDA

(a) CNN, MNIST3

(b) CNN, CIFAR10

Figure 4: Impact of training data size. Mean p-values of TDA methods with
randomness induced by SWA+DE-Init.

Based on the observation in §4.2 that TDA values are often dominated by noise,
we delve into the factors that lead to the instability of data attributions.
We inspect the contribution of model initialisation, training set size and
model complexity.

#### Source of randomness.

From a Bayesian ML perspective, the stochasticity of TDA stems from the
inherent uncertainty of the learned model posterior $p(\theta|\mathcal{D})$.
We consider two sources of randomness, model initialisation (DE-Init) and SGD
batch composition (DE-Batch). Results are reported in Table 1. For MNIST3 and
CIFAR10, we observe that DE-Batch introduces lower levels of noise in the TDA
estimates (lower p-values). Particularly on MNIST3, both LOO and other TDA
methods result in statistically significant p-values ($<0.005$). This implies
that almost every training data $z_{j}$ is influencing every test data $z$
consistently across various batch compositions. We conclude that the greater
source of variations for the TDA estimates is the model initialisation.

#### Training set size.

We study how training set size is a source of noise (cf. Figure 4). We train
the CNN with different-size datasets of MNIST3 and CIFAR10, where we vary the
number of samples per class. To limit the effect of noise introduced by new
data, we subsample the smaller datasets from the larger ones. Batches are
composed differently depending on the dataset size, meaning that parameter
updates are made after processing different data. The results show a tendency
for high variation in TDA scores with larger datasets. This makes sense as the
number of combinations for batch composition increases with dataset size. As
the batching is initialised randomly during training, batches are likely to be
composed of different data for larger datasets. This leads to variation in the
learned model parameters, in turn affecting the reliability of TDA.

Table 2: Impact of model complexity. Mean p-values of the ground-truth TDA $\tau(z_{j},z)$ (LOO) and estimated TDA values $\tau^{\prime}(z_{j},z)$ (the other 4 columns) with randomness induced by SWA+DE-Batch for MNIST3 and DE-Batch for CIFAR10. Model | Data | LOO | ATS | IF | GD | GC
---|---|---|---|---|---|---
CNN | MNIST3 | 0.025 | 0.039 | 0.000 | 0.000 | 0.000
ViT+LoRA | MNIST3 | 0.786 | 0.573 | 0.369 | 0.365 | 0.093
CNN | CIFAR10 | 0.623 | 0.374 | 0.535 | 0.534 | 0.314
ViT+LoRA | CIFAR10 | 0.777 | 0.766 | 0.686 | 0.686 | 0.522

#### Model complexity.

We study how model complexity is linked to the reliability of TDA estimates.
See Table 2. We observe that, compared to a small CNN model, a large ViT model
trained with LoRA results in dramatically greater p-values. For example, for
LOO on MNIST3, the p-value increases from 0.025 to 0.786. A similar trend is
observed for other TDA methods. This implies that the reliability of TDA
estimates decreases with increasing model complexity. While we limit the
number of trainable parameters in our ViT by using LoRA to be comparable to
e.g. adding another layer in the CNN, the p-values computed from TDA estimates
are significantly larger. Larger models exhibit a larger parameter space so
that noise stemming from model initialisation or batch composition is
amplified. While we fix the model initialisation and dataset size, the batch
composition still varies across the model parameters $\theta$ sampled from the
posterior $p(\theta|\mathcal{D})$ per model. As both the CNN and ViT are
trained with the same sampled batch compositions, we attribute the strong
increase of p-value to the model complexity.

###  4.4 (Dis)agreement between TDA methods

Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 5: Correlation of TDA methods. Pearson and Spearman correlation
coefficients among ground-truth TDA and approximate TDA methods. We show
correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$,
and TDA p-values. All results are based on the setting: CNN, MNIST3, SWA+DE-
Init.

We test the reliability of different TDA methods. Ideally, all methods
approximate the ground-truth TDA (LOO). Yet the results suggest that there are
substantial differences among the methods. For example, Grad-Cos is much more
stable than all others. Hence, we study TDA methods with respect to both their
correlation with LOO and among each other using Pearson and Spearman
correlation of mean and variance of the TDA distributions, as proposed in §3.

Figure 5 shows the correlation matrices for one experiment (all experimental
results in the Appendix). The results show that viewing TDA scores as
distributions gives insights into the reliability of TDA methods: None of the
tested TDA method’s expected values $\hat{\mu}$ correlates with LOO. This
implies that none of the TDA methods is a good approximation when the random
factors are considered. The poor correlation of p-values of LOO indicates a
disagreement in the train-test pairs considered low noise. We conclude that
none of the tested methods reliably capture ground-truth TDA distributions.

Interestingly, we notice a stronger correlation between all other methods
particularly when looking at the correlations of p-values. We identify two
groups based on positive correlation, i.e. ATS with IF and GD with GC. Among
the two groups, there is a negative correlation which indicates that methods
interpret the sign of the attribution differently. Between IF, GD and GC this
makes sense, as there is a negative sign in the definition of IF (Equation 2)
which is not present in GD and GC. Considering absolute correlation, IF and GD
are strongly correlated which shows that the dot product is a valid
alternative for IF as they produce similar score distributions. The
correlation between GD and GC indicates that the normalisation of gradients
does not have a strong impact on the estimated TDA.

IF, GD and GC correlate considerably with ATS, which measures how the loss of
a model on $z$ changes after doing one additional training step on $z_{j}$.
Practically, ATS represents the gradient update after $z_{j}$, which is the
same as the gradient $z_{j}$ itself. Therefore, it makes sense that the
gradient-based approximation methods are close to ATS. We recognise a
difference in the scope LOO and ATS address. LOO looks at TDA globally and
encapsulates the whole training, whereas ATS considers a local scope with a
small model change. As IF, GD and GC correlate with ATS, we observe that they
also correspond to a local change in the model, which underlines and extends
the argument of [7]: There is a gap between LOO and IF, and more generally
between the global and local view on TDA.

The TDA variance $\hat{\sigma}$ is noticeably well-correlated for TDA
estimators and LOO, except for GC. This implies the existence of a consistent
ranking of train-test pairs with stable attribution relationships. In
particular, stable train-test pairs predicted by LOO are also likely to be
stable pairs for TDA methods like ATS, IF, and GD. This motivates our final
analysis and recommendation for evaluation in §4.5.

###  4.5 Considerations on TDA evaluation from a Bayesian perspective

Our analysis shows that both TDA estimates and ground-truth TDA values are
affected by the noise stemming from the stochastic nature of model training.
Hence, the practice of comparing against such a ground truth is destined to
result in fragile estimates. We propose to treat TDA estimates as random
variables which allows us to look at the evaluation from a Bayesian
perspective: The comparison of TDA estimates against target TDA values is a
comparison of two random variables. Since it is impossible to get rid of
noise, it is better to compare distributions rather than point estimates. This
provides an understanding of how well methods approximate the ground-truth
distribution.

We observe that p-values vary between individual train-test sample pairs
$(z_{j},z)$; not all TDA estimates are equally affected by stochasticity.
Interestingly, the presence of low-noise pairs is consistent across the
majority of our experiments (cf. Figure 3), with varying sizes of the low-
noise fraction. We find that fixing model initialisation and a small dataset
size gives rise to a larger number of low-noise pairs. We propose to focus on
such low-noise pairs in TDA evaluation as their estimates are low in variance,
leading to a more reliable evaluation. Identifying such pairs requires an
analysis similar to this work: treating TDA values as distributions and
sampling multiple times from the posterior to get an estimate of the noise. It
is crucial to find low-noise pairs to base evaluations on and understand when
TDA is applicable. If no low-variance pairs exist, TDA cannot be used.

##  5 Related work

We study the reliability of TDA methods and add to the existing body of work
on the fragility of TDA methods. Previous studies focused primarily on IF [27,
6, 8, 9, 7]. We extend the analysis by additionally studying other TDA
methods. While IFs are theoretically grounded in robust statistics, they are
based on two assumptions which are not always fulfilled in the context of deep
learning: Twice-differentiability and strict convexity of the loss [3]. Zhang
& Zhang [27] and Basu et al. [6] point to the fragility of the influence
scores due to the non-convexity of deep learning. Particularly increasing
model size is connected to increased model curvature, which means that
influence estimates are more fragile with larger models. They find that strong
regularisation is needed to improve estimation quality. Our experiments verify
the observation that fragility increases with model size, which we observe
across methods. We add that sources of randomness in the training process
attribute to the fragility of TDA methods with increasing model size.
Furthermore, related work found that the size of the training set contributes
to the fragility of influence estimates. The attribution of one sample in a
large training set is marginal so both influence estimates and ground-truth
influence scores (i.e., from retraining the model) are noisy [6, 8, 9].
Through a Bayesian lens, we connect the increased fragility with increasing
dataset size to batch composition as well. Not only is the attribution of a
single sample in a large dataset marginal [6] but batches have vastly
different compositions in larger datasets, introducing noise. A recent work
[7] states that influence functions in deep learning do not correspond to LOO
and quantify gaps in the estimation stemming from model non-linearity. K &
Søgaard [9] recommend reporting expected TDA scores to increase estimation
stability. This approach is closest to our work but misses the consideration
of variance in TDA estimates which we include by taking a Bayesian viewpoint.

In contrast to related work, we treat TDA values as distributions, which
enables a novel perspective on the TDA task. We highlight the importance of
considering the variance when studying reliability.

##  6 Conclusion

We adopt a Bayesian perspective on the training data attribution (TDA) methods
to study their reliability, given the stochastic nature of deep model
training. By modelling TDA scores as distributions, we find that randomness in
the training process, particularly due to parameter initialisation and batch
composition, translates to variation in ground-truth TDA. We empirically
observe that current estimation methods, such as influence functions, model a
local change in the model whereas the ground truth attribution considers a
global model change. Therefore, TDA is subject to inherent variance, leading
us to suggest to the community: (1) When proposing a novel TDA method, one
should view TDA from a Bayesian perspective and study the TDA estimates as
distributions. (2) When using TDA, one should consider the variance to
understand when the estimate can be trusted.

#### Limitations.

We perform an exhaustive analysis of TDA values $\tau$ and the estimates
$\tau^{\prime}$ for all train-test pairs $(z_{j},z)$. Because of considerable
computational costs, we have subsampled the datasets. In practice, datasets
are considerably larger. Moreover, we choose simple tasks to eliminate the
need for an additional hyperparameter search for model training, as the
principal focus is on studying TDA methods. We choose gradient-based TDA
methods but acknowledge that there exist many more, that we do not address.
Hence, we encourage further study of TDA methods to fill these gaps and
recommend investigating TDA from a Bayesian perspective, particularly in the
low-data regime.

#### Broader impact.

This paper contributes to the field of data-driven XAI which aims at helping
humans understand the inner workings of opaque models through data-centric
approaches. Our work contributes to understanding the reliability of TDA
methods and rethinking their evaluation against a noisy ground truth, which
could help assess when TDA is appropriate and reliable.

## Acknowledgments and Disclosure of Funding

Kay Choi has helped in designing Figures 1 and 2. The authors thank the
International Max Planck Research School for Intelligent Systems (IMPRS-IS)
for supporting Elisa Nguyen.

## References

  * Guidotti et al. [2018] Riccardo Guidotti, Anna Monreale, Salvatore Ruggieri, Franco Turini, Fosca Giannotti, and Dino Pedreschi.  A survey of methods for explaining black box models.  _ACM Comput. Surv._ , 51(5), 2018.  ISSN 0360-0300.  doi: 10.1145/3236009.  URL https://doi.org/10.1145/3236009.
  * Hammoudeh and Lowd [2022] Zayd Hammoudeh and Daniel Lowd.  Training data influence analysis and estimation: A survey.  _ArXiv_ , abs/2212.04612, 2022.
  * Koh and Liang [2017] Pang Wei Koh and Percy Liang.  Understanding black-box predictions via influence functions.  In Doina Precup and Yee Whye Teh, editors, _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 1885–1894. PMLR, 06–11 Aug 2017.  URL https://proceedings.mlr.press/v70/koh17a.html.
  * Pruthi et al. [2020] Garima Pruthi, Frederick Liu, Satyen Kale, and Mukund Sundararajan.  Estimating training data influence by tracing gradient descent.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 19920–19930. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper_files/paper/2020/file/e6385d39ec9394f2f3a354d9d2b88eec-Paper.pdf.
  * Yang et al. [2023] Jinghan Yang, Sarthak Jain, and Byron Wallace.  How many and which training points would need to be removed to flip this prediction?  In _Proceedings of the 17th Conference of the European Chapter of the Association for Computational Linguistics_ , pages 2563–2576, Dubrovnik, Croatia, May 2023. Association for Computational Linguistics.  URL https://aclanthology.org/2023.eacl-main.188.
  * Basu et al. [2021] Samyadeep Basu, Phil Pope, and Soheil Feizi.  Influence functions in deep learning are fragile.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=xHKVVHGDOEk.
  * Bae et al. [2022] Juhan Bae, Nathan Ng, Alston Lo, Marzyeh Ghassemi, and Roger B Grosse.  If influence functions are the answer, then what is the question?  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 17953–17967. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/7234e0c36fdbcb23e7bd56b68838999b-Paper-Conference.pdf.
  * Epifano et al. [2023] Jacob Epifano, Ravichandran Ramachandran, Aaron J. Masino, and Ghulam Rasool.  Revisiting the fragility of influence functions.  _Neural networks : the official journal of the International Neural Network Society_ , 162:581–588, 2023.
  * K and Søgaard [2021] Karthikeyan K and Anders Søgaard.  Revisiting methods for finding influential examples.  _CoRR_ , abs/2111.04683, 2021.  URL https://arxiv.org/abs/2111.04683.
  * Welling and Teh [2011] Max Welling and Yee W Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th international conference on machine learning (ICML-11)_ , pages 681–688, 2011.
  * Lakshminarayanan et al. [2017] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  _Advances in neural information processing systems_ , 30, 2017.
  * Izmailov et al. [2018] Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson.  Averaging weights leads to wider optima and better generalization.  In _34th Conference on Uncertainty in Artificial Intelligence 2018, UAI 2018_ , 34th Conference on Uncertainty in Artificial Intelligence 2018, UAI 2018, pages 876–885. Association For Uncertainty in Artificial Intelligence (AUAI), 2018.
  * Brunet et al. [2019] Marc-Etienne Brunet, Colleen Alkalay-Houlihan, Ashton Anderson, and Richard Zemel.  Understanding the origins of bias in word embeddings.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 803–811. PMLR, 09–15 Jun 2019.  URL https://proceedings.mlr.press/v97/brunet19a.html.
  * Akyurek et al. [2022] Ekin Akyurek, Tolga Bolukbasi, Frederick Liu, Binbin Xiong, Ian Tenney, Jacob Andreas, and Kelvin Guu.  Towards tracing knowledge in language models back to the training data.  In _Findings of the Association for Computational Linguistics: EMNLP 2022_ , pages 2429–2446, Abu Dhabi, United Arab Emirates, December 2022\. Association for Computational Linguistics.  URL https://aclanthology.org/2022.findings-emnlp.180.
  * Guo et al. [2021] Han Guo, Nazneen Rajani, Peter Hase, Mohit Bansal, and Caiming Xiong.  FastIF: Scalable influence functions for efficient model interpretation and debugging.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 10333–10350, Online and Punta Cana, Dominican Republic, November 2021. Association for Computational Linguistics.  doi: 10.18653/v1/2021.emnlp-main.808.  URL https://aclanthology.org/2021.emnlp-main.808.
  * Schioppa et al. [2022] Andrea Schioppa, Polina Zablotskaia, David Vilar, and Artem Sokolov.  Scaling up influence functions.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(8):8179–8186, Jun. 2022.  doi: 10.1609/aaai.v36i8.20791.  URL https://ojs.aaai.org/index.php/AAAI/article/view/20791.
  * Charpiat et al. [2019] Guillaume Charpiat, Nicolas Girard, Loris Felardos, and Yuliya Tarabalka.  Input similarity from the neural network perspective.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/c61f571dbd2fb949d3fe5ae1608dd48b-Paper.pdf.
  * Khanna et al. [2019] Rajiv Khanna, Been Kim, Joydeep Ghosh, and Sanmi Koyejo.  Interpreting black box predictions using fisher kernels.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 3382–3390. PMLR, 2019.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * Gal and Ghahramani [2016] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _international conference on machine learning_ , pages 1050–1059. PMLR, 2016.
  * Maddox et al. [2019] Wesley J Maddox, Pavel Izmailov, Timur Garipov, Dmitry P Vetrov, and Andrew Gordon Wilson.  A simple baseline for bayesian uncertainty in deep learning.  _Advances in neural information processing systems_ , 32, 2019.
  * LeCun and Cortes [2005] Yann LeCun and Corinna Cortes.  The mnist database of handwritten digits.  2005\.
  * Krizhevsky [2009] Alex Krizhevsky.  Learning multiple layers of features from tiny images.  2009\.
  * LeCun et al. [1998] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Dosovitskiy et al. [2021] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=YicbFdNTTy.
  * Hu et al. [2022] Edward J Hu, yelong shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen.  LoRA: Low-rank adaptation of large language models.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=nZeVKeeFYf9.
  * Zhang and Zhang [2022] Rui Zhang and Shihua Zhang.  Rethinking influence functions of neural networks in the over-parameterized regime.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 36(8):9082–9090, Jun. 2022.  doi: 10.1609/aaai.v36i8.20893.  URL https://ojs.aaai.org/index.php/AAAI/article/view/20893.
  * Mangrulkar et al. [2022] Sourab Mangrulkar, Sylvain Gugger, Lysandre Debut, Younes Belkada, and Sayak Paul.  Peft: State-of-the-art parameter-efficient fine-tuning methods.  https://github.com/huggingface/peft, 2022.
  * Wolf et al. [2020] Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander M. Rush.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online, October 2020. Association for Computational Linguistics.  URL https://www.aclweb.org/anthology/2020.emnlp-demos.6.
  * Wu et al. [2020] Bichen Wu, Chenfeng Xu, Xiaoliang Dai, Alvin Wan, Peizhao Zhang, Zhicheng Yan, Masayoshi Tomizuka, Joseph Gonzalez, Kurt Keutzer, and Peter Vajda.  Visual transformers: Token-based image representation and processing for computer vision, 2020.

##  Appendix A Model training details

We provide the source code at https://github.com/ElisaNguyen/bayesian-tda. All
experiments were run on a single Nvidia 2080ti GPU.

###  A.1 Data sampling

We use subsampled versions of the openly available MNIST [22] and CIFAR10 [23]
datasets. For this, we first add an index which we use for randomly sampling a
fixed number of images from each class. Table 3 includes the dataset sizes of
the different experiments.

###  A.2 Model training

The CNN has two convolutional layers followed by two fully connected linear
layers, with GeLU activation. We use the Adam optimizer with a learning rate
of 0.001 and a weight decay of 0.005. We use the cross-entropy loss and train
the model for 15 epochs on MNIST3 and for 30 epochs on CIFAR10 with a batch
size of 32. For training the ViT with LoRA, we use the peft [28] and
HuggingFace transformers library [29]. We start from the pretrained model
checkpoint of [30] and finetune the LoRA model with the same hyperparameters
as the CNN.

An overview of the predictive performance measured in accuracy on the
subsampled training and test sets is provided in Table 3.

#### Hint for reproducibility.

In particular, we use CrossEntropyLoss(reduction=‘none’) during model training
and also update this in the ViT training script modeling_vit.py. This is
important for the LOO experiments, where we exclude a sample $z_{j}$ from
contributing to the training by zeroing out the loss.

Table 3: Predictive performance at 95% CI across 10 runs (computed as 1.96*SE) Experiment |  |
---|---|---
Model | Data | Randomness | $|\mathcal{D}_{\text{train}}|$ | $|\mathcal{D}_{\text{test}}|$ | Accuracy${}_{\text{train}}$ | Accuracy${}_{\text{test}}$
CNN | MNIST3 | SWA+DE-Init | 30 | 900 | 0.987$\pm$0.010 | 0.953$\pm$0.003
CNN | MNIST3 | SWA+DE-Init | 60 | 900 | 0.985$\pm$0.007 | 0.970$\pm$0.004
CNN | MNIST3 | SWA+DE-Init | 150 | 900 | 0.998$\pm$0.004 | 0.970$\pm$0.003
CNN | CIFAR10 | SWA+DE-Init | 100 | 500 | 0.989$\pm$0.020 | 0.260$\pm$0.010
CNN | CIFAR10 | SWA+DE-Init | 200 | 500 | 1.000$\pm$0.000 | 0.328$\pm$0.007
CNN | CIFAR10 | SWA+DE-Init | 500 | 500 | 0.992$\pm$0.014 | 0.377$\pm$0.012
CNN | MNIST3 | SWA+DE-Batch | 150 | 900 | 0.993$\pm$0.002 | 0.975$\pm$0.002
CNN | CIFAR10 | SWA+DE-Batch | 500 | 500 | 0.991$\pm$0.010 | 0.364$\pm$0.003
ViT+LoRA | MNIST3 | SWA+DE-Batch | 150 | 900 | 0.945$\pm$0.008 | 0.935$\pm$0.005
ViT+LoRA | CIFAR10 | SWA+DE-Batch | 500 | 500 | 0.934$\pm$0.006 | 0.892$\pm$0.008

##  Appendix B Additional experimental results

We study the reliability of TDA estimates and values through a hypothesis
test, where we report the p-values as an indication of the statistical
significance of the TDA estimate. In this appendix, we provide a complete
overview of the analyses of p-values and correlations between different TDA
methods.

###  B.1 DE vs. DE+SWA

In our work, we sample trained models $\theta$ sampled from the posterior
$p(\theta|\mathcal{D})$. Concretely, we train the model across 10 different
random seeds and record the checkpoints after each of the last five epochs of
training. Each of these models represents a sample. We test how the stability
of the TDA values $\tau(z_{j},z)$ (LOO) behaves when we ensemble different
numbers of checkpoints by investigating the mean p-values across all train-
test pairs $(z_{j},z)$. Figure 6 shows that higher numbers of samples $\theta$
increase stability, therefore we use all available samples in our subsequent
analyses.

MNIST3 | CIFAR10
---|---
Figure 6: P-values of ground-truth TDA (LOO) with increasing number of SWA
samples (i.e., number of model checkpoints used as samples of $\theta$) for
the CNN trained on MNIST3 and CIFAR10.

###  B.2 All results: Stability of TDA values and estimates

Table 4 presents the complete table of p-values of all tested TDA methods
across all experiments of our work. Below, we display the histograms of
p-values per experiment corresponding to each line in the table captioned with
the experiment ID. The histograms show that low-noise train-test pairs
$(z_{j},z)$ are present in all experiments involving the CNN model (i.e.,
experiments 1-8), where the number of pairs varies. Generally, we observe that
there is no connection between the size of the dataset and the distribution of
p-values. Furthermore, we notice that fixing the model initialisation (i.e.,
randomness induced by SWA+DE-Batch) increases the number of stable train-test
pairs (cf. experiment 3 to 7, 4 to 8). However, in the case of the ViT
experiments, stable train-test pairs are practically non-existent which shows
that model complexity affects the stability of TDA.

Table 4: Complete list of mean p-values of TDA values for all experiments. Experiment |  |  |  |  |
---|---|---|---|---|---
ID | Model | Data | Randomness | $|\mathcal{D}_{\text{train}}|$ | LOO | ATS | IF | GD | GC
1 | CNN | MNIST3 | SWA+DE-Init | 30 | 0.058 | 0.088 | 0.111 | 0.110 | 0.020
2 | CNN | MNIST3 | SWA+DE-Init | 60 | 0.421 | 0.148 | 0.251 | 0.253 | 0.002
3 | CNN | MNIST3 | SWA+DE-Init | 150 | 0.331 | 0.254 | 0.352 | 0.363 | 0.003
4 | CNN | CIFAR10 | SWA+DE-Init | 100 | 0.665 | 0.424 | 0.607 | 0.608 | 0.352
5 | CNN | CIFAR10 | SWA+DE-Init | 200 | 0.552 | 0.397 | 0.450 | 0.452 | 0.399
6 | CNN | CIFAR10 | SWA+DE-Init | 500 | 0.692 | 0.438 | 0.575 | 0.587 | 0.356
7 | CNN | MNIST3 | SWA+DE-Batch | 150 | 0.025 | 0.039 | 0.000 | 0.000 | 0.000
8 | CNN | CIFAR10 | SWA+DE-Batch | 500 | 0.623 | 0.374 | 0.535 | 0.534 | 0.314
9 | ViT+LoRA | MNIST3 | SWA+DE-Batch | 150 | 0.786 | 0.573 | 0.369 | 0.365 | 0.093
10 | ViT+LoRA | CIFAR10 | SWA+DE-Batch | 500 | 0.777 | 0.766 | 0.686 | 0.686 | 0.522
1 | 2 | 3 | 4
---|---|---|---
5 | 6 | 7 | 8
9 | 10 |  |
Figure 7: Distribution of p-values for ground-truth TDA (LOO) for all
experiments (IDs corresponding to IDs in Table 4).

###  B.3 All results: Correlation analysis

In the main body of this paper, we report the Pearson and Spearman correlation
matrices for experiment 3 (CNN trained on MNIST3 with 50 samples per class and
randomness induced by SWA+DE-Init). This section presents the complete
overview of correlations between TDA methods across all experiments in Figures
8 \- 17. We note that the observations and analyses in the main paper hold
across experiments.

Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 8: Experiment 1 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 9: Experiment 2 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 10: Experiment 3 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 11: Experiment 4 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 12: Experiment 5 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 13: Experiment 6 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 14: Experiment 7 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 15: Experiment 8 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 16: Experiment 9 (cf. Table 4): Pearson and Spearman correlation coefficients among ground-truth TDA and approximate TDA methods. We show correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$, and TDA p-values. Pearson’s $r$ |  |  |
---|---|---|---
Spearman’s $\rho$ |  |  |
|  Mean $\hat{\mu}$ |  Standard deviation $\hat{\sigma}$ | p-value
Figure 17: Experiment 10 (cf. Table 4): Pearson and Spearman correlation
coefficients among ground-truth TDA and approximate TDA methods. We show
correlations for TDA mean $\hat{\mu}$, TDA standard deviation $\hat{\sigma}$,
and TDA p-values.


11institutetext: Department of Physics and Astronomy, Clemson University,
Kinard Lab of Physics, Clemson, SC 29634, USA 22institutetext: Harvard-
Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138,
USA 33institutetext: INAF - Osservatorio di Astrofisica e Scienza dello Spazio
di Bologna, Via Piero Gobetti, 93/3, 40129, Bologna, Italy

# A New Mid-Infrared and X-ray Machine Learning Algorithm to Discover Compton-
thick AGN

R. Silver 11 N. Torres-Albà 11 X. Zhao 22 1 1 S. Marchesi 33 1 1 A. Pizzetti
11 I.Cox 11 M. Ajello 11

We present a new method to predict the line-of-sight column density (NH)
values of active galactic nuclei (AGN) based on mid-infrared (MIR), soft, and
hard X-ray data. We developed a multiple linear regression machine learning
algorithm trained with WISE colors, Swift-BAT count rates, soft X-ray hardness
ratios, and an MIR$-$soft X-ray flux ratio. Our algorithm was trained off 451
AGN from the Swift-BAT sample with known NH and has the ability to accurately
predict NH values for AGN of all levels of obscuration, as evidenced by its
Spearman correlation coefficient value of 0.86 and its 75% classification
accuracy. This is significant as few other methods can be reliably applied to
AGN with Log(N${}_{\rm H}<$) 22.5. It was determined that the two soft X-ray
hardness ratios and the MIR$-$soft X-ray flux ratio were the largest
contributors towards accurate NH determination. This algorithm will contribute
significantly to finding Compton-thick (CT-) AGN (N${}_{\rm H}\geq$ 1024
cm-2), thus enabling us to determine the true intrinsic fraction of CT-AGN in
the local universe and their contribution to the Cosmic X-ray Background.

###### Key Words.:

Galaxies: active – Galaxies: nucleus – Infrared: galaxies – X-rays: galaxies

##  1 Introduction

Active Galactic Nuclei (AGN) are supermassive black holes (SMBHs) that reside
in the center of nearly all massive galaxies and accrete nearby material.
These are one of the most powerful sources classes in the Universe, and emit
over the entire electromagnetic spectrum. It has been shown that the masses of
the SMBHs correlate with that of the host galaxy bulge, velocity dispersion,
and luminosity (Magorrian et al., 1998; Richstone et al., 1998; Gebhardt et
al., 2000; Merritt & Ferrarese, 2001; Ferrarese & Ford, 2005; Kormendy & Ho,
2013). This trend indicates that SMBHs may determine star formation rates, due
to molecular and ionized outflows (Ferrarese & Merritt, 2000; Gebhardt et al.,
2000; Di Matteo et al., 2005; Merloni et al., 2010; Fiore, F. et al., 2017;
Martín-Navarro et al., 2018). If true, then the cosmic evolution of SMBHs and
their host galaxies are inextricably linked. Therefore, being able to study
the properties of SMBHs, including the gas and dust that surrounds them,
becomes crucial.
One of the best ways to study AGN through cosmic time is the cosmic X-ray
background (CXB), i.e., the diffuse X-ray emission from 1 to 200–300 keV
(e.g., Alexander et al., 2003; Gilli et al., 2007; Treister et al., 2009; Ueda
et al., 2014; Brandt & Yang, 2021). Models have shown that a significant
fraction (15-20%; Gilli et al., 2007; Ananna et al., 2019) of the peak of the
CXB ($\sim$30 keV, Ajello et al., 2008) is generated by a population of AGN
with large obscuring column densities, NH,los $\geq$ 1024 cm-2, labeled as
Compton-thick (CT-) AGN. Additionally, population synthesis models designed to
accurately describe the origins of the CXB estimate that between 20% (Ueda et
al., 2014) and 50% (Ananna et al., 2019) of all AGN are CT. Nonetheless, the
current fraction of observed CT-AGN is only between 5% and 10%, even at low
redshifts (i.e., z $<$ 0.01; Burlon et al., 2011; Ricci et al., 2015; Torres-
Albà et al., 2021).
CT-AGN are challenging to discover because the majority of their emission,
from the optical through the soft X-rays, is obscured by the surrounding dust
and gas (i.e., the torus; Urry & Padovani, 1995). However, the hard X-rays
($>$10 keV) and the mid-infrared (MIR, 3–30$\mu$m) are able to pierce through
the torus up to high column densities, making them the least biased bands
against the detection of heavily obscured AGN (Treister et al., 2004; Stern et
al., 2005; Alexander et al., 2008). The hard X-ray emission is created when UV
light from the accretion disk interacts with hot electrons in the corona above
the disk, thus Compton up-scattering into the hard X-ray band (Haardt &
Maraschi, 1993). Additionally, the same UV radiation is absorbed by the dust,
which in turn emits thermally in the IR (Almeida & Ricci, 2017; Hönig, 2019).
Because of this, the emission in these two bands is expected to correlate
significantly in AGN. Therefore, targeting the X-rays and the MIR is the ideal
way to discover new CT-AGN.
Observing and analyzing spectra of AGN in the X-rays and infrared to identify
strong CT candidates is a time and resource intensive endeavor (see, e.g.,
Marchesi et al., 2017; Andonie et al., 2022; Silver et al., 2022). The Burst
Alert Telescope (BAT) on board Swift (Gehrels et al., 2004) detected 1390
sources in its 150-month
catalog111https://science.clemson.edu/ctagn/bat-150-month-catalog/, almost 500
more than its predecessor, the 100-month catalog. With the addition of
hundreds of sources in every catalog release, an efficient and accurate method
to identify potential heavily obscured AGN is necessary. For this reason, our
team has developed a new multiple linear regression machine learning algorithm
to predict the line-of-sight column density of AGN. We have constructed a
large sample of AGN with known NH values and trained the algorithm using their
MIR from the Wide-field Infrared Survey Explorer (WISE, Wright et al., 2010),
soft X-ray data from the Swift-X-ray Telescope (XRT), and hard X-ray data from
Swift-BAT, to accurately predict the column density of new AGN. This work
proceeds as follows: Section 2 discusses the creation of our sample while
Section 3 describes how the NH values were determined. Section 4 details the
algorithm implemented and the input parameters included. Section 5 discusses
the results of our algorithm and compares our predictive capabilities to other
recent methods based on linear data modeling, rather than on multi-parameter
machine learning algorithms like what is presented in this paper.

##  2 Sample Selection

Our sample is taken from the 1390 sources detected in the BAT 150 Month
catalog222The online version of the catalog can be found at:
https://science.clemson.edu/ctagn/bat-150-month-catalog/ (Imam et al. in
preparation). Of these 1390 sources, 568 are AGN with reliable NH
determinations (see Section 3 for details). Asmus et al. (2015) showed that
the ratio of the MIR and X-ray flux can be a strong predictor of column
density. Therefore, our machine learning algorithm (see Section 4.2) includes
data from XRT and WISE, so we cross-matched (using the BAT counterpart
coordinates) with the 2SXPS (Evans et al., 2020) and AllWISE (Cutri et al.,
2013) with 5$\arcsec$ and 10$\arcsec$, respectively. For the 2SXPS and
AllWISE, we found an average separation of $\sim$1.7$\arcsec$ and
$\sim$1.8$\arcsec$, and a standard deviation of $\sim$1.1$\arcsec$ and
$\sim$1.5$\arcsec$, respectively. This left us with a sample of 451 sources to
train and test our machine learning algorithm (see Section 4).

##  3 Data Analysis

The majority of the sources (361) in our sample of 451 are in the BAT 70-month
catalog (Ricci et al., 2017), which provides NH values based on spectral
analysis of soft X-ray (ASCA, Chandra, Suzaku, Swift-XRT, and XMM-Newton) and
BAT spectra. For the remaining 90 sources, we modeled their soft X-ray jointly
with their Swift-BAT spectra. XMM-Newton data was available for 18 sources,
while Chandra data was available for an additional 24. For the remaining 48
sources, the soft X-ray data were provided by Swift-XRT. As the greater part
of the sources in the sample were unobscured (Log(N${}_{\rm H})<$ 22) or
mildly obscured (22 $<$ Log(N${}_{\rm H})<$ 23), they were sufficiently
modeled with an absorbed powerlaw, as shown below:

$Model1=constant_{1}*phabs*(zphabs*zpowerlw),$ (1)

However, Compton-thin (23 $<$ Log(N${}_{\rm H})<$ 24) sources required a more
complex model to account for the Fe K$\alpha$ emission and the fraction of
intrinsic emission that leaks through the torus rather than being absorbed by
the obscuring material. These sources were modeled as such:

$Model2=constant_{1}*phabs*(zphabs*zpowerlw\\\
+zgauss+constant_{2}*zpowerlw),$ (2)

where $constant_{1}$ accounts for cross-normalization differences between the
soft X-ray instrument and Swift-BAT, $phabs$ models the galactic absorption,
$zphabs*zpowerlw$ is the absorbed power-law modeling the intrinsic emission,
$zgauss$ models the Fe K$\alpha$ emission line, and $constant_{2}*zpowerlw$
represents the scattered emission that leaks through the torus.
When sources approach or surpass the Compton-thick limit, they require even
more sophisticated modeling. These sources were modeled with physically
motivated models such as MYTorus (Murphy & Yaqoob, 2009) and borus02
(Baloković et al., 2018), and have been described in detail in Zhao et al.
(2019a, b); Torres-Albà et al. (2021); Silver et al. (2022). These models are
used for heavily obscured AGN because they account for the photons that
interact with the dust and gas surrounding the SMBH and are reflected into the
observer line of sight.

##  4 Machine Learning

###  4.1 Multiple Linear Regression

Linear regression is one of the most commonly used machine learning techniques
(see, e.g., Chen et al., 2021; Mizukoshi et al., 2022). Simply, linear
regression models the linear relationship between an explanatory variable
(input parameter) and the response variable (output parameter). Since few
quantities can be accurately modeled using only one explanatory variable,
using numerous can improve the predictive capability of an algorithm. This is
referred to as multiple linear regression, or just multiple regression for
short, and is modeled as shown below:

$y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...\beta_{i}x_{i},$ (3)

where y is the response variable, $x_{i}$ are the explanatory variables,
$\beta_{0}$ is the y-intercept (if necessary), and $\beta_{i}$ are the slope
coefficients corresponding to each explanatory variable. Using a large sample
of sources with data for every explanatory variable and a known value for the
response variable, the algorithm trains itself to determine which combination
of $\beta_{i}$ values is optimal for reproducing the response variable. Out of
our sample of 451 sources, 80% were randomly selected to be used in our
training sample, thus leaving 20% (91 sources) left for our test sample. We
determined this was the optimal ratio as it used enough sources to accurately
train the algorithm while simultaneously leaving a statistically significant
sample to verify this accuracy 333We note that neural networks are another
commonly used machine learning algorithm (see, e.g., Finke et al., 2021;
Chainakun et al., 2022; Zubovas et al., 2022). We applied this technique to
our data set and after finding the optimal configuration, yielded very similar
results to those generated by our linear regression model. For this reason, we
have optioned to present the results from the comparatively simpler linear
regression model in this paper..

###  4.2 Parameters Used

The algorithm will only be as accurate as however strong the relationship is
between the chosen input parameters and the desired output parameter. We have
selected WISE colors, BAT count rates, soft X-ray hardness ratios (HRs), and
an MIR$-$soft X-ray flux ratio as all have been previously shown to correlate
with NH. These parameters are described below.

####  4.2.1 MIR Colors

Roughly half of the intrinsic emission from the AGN is absorbed by the dusty
torus (see, e.g., Almeida & Ricci, 2017; Hönig, 2019). As a consequence, the
dust present in the torus is heated to temperatures of several hundred Kelvin,
and thus radiates thermally. This emission peaks in the MIR ($\sim$3–30$\mu$m)
and is much less prone to absorption than the optical and UV, making it a
crucial tool to study obscured AGN. With the launch of WISE, we have an all-
sky instrument with superb resolution ($\sim$6”) capable of studying these
obscured sources. WISE observes the entire sky in four bands, 3.4, 4.6, 12,
and 22 $\mu$m (W1, W2, W3, and W4, respectively), and to date has detected
nearly 750 million sources and reached flux limits of 7.1$\times$10-14 erg
cm-2 s-1. The differences between these bands have been proven to be a good
predictor for different levels of obscuration. Kilerci Eser et al. (2020) used
a sample of AGN from the BAT 105-month catalog (Oh et al., 2018) to create new
CT-AGN selection criteria based on MIR colors. They find that the median
values for different colors have an increasing trend with NH (see Figure 9
from their paper). Therefore, our algorithm includes six WISE colors: W1-W2,
W1-W3, W1-W4, W2-W3, W2-W4, and W3-W4.

####  4.2.2 MIR–X-ray Flux Ratio

As the X-ray and MIR emission are both reprocessed from the same material, it
is expected that a correlation exists between them. Asmus et al. (2015) shows
the trend between the observed 12 $\mu$m flux and the 2–10 keV flux (see
Figure 1). Moreover, a shift is evident in the trend based on obscuration. As
source obscuration increases, the observed 2–10 keV flux decreases, thus
causing the source to fall to the left of the predicted trend. This is
evidenced in the figure by Seyfert 2 galaxies (red squares) falling to the
left of traditionally unobscured Seyfert 1 galaxies (blue circles). Moreover,
the confirmed CT-AGN (black stars) fall well to the left of even Seyfert 2
galaxies, suggesting an extremely suppressed observed X-ray flux. As a result
of this trend, Asmus et al. (2015) used the log ratio of the 12 $\mu$m flux
density and the 2–10 keV flux to predict the column density of an AGN. We have
included this parameter in our algorithm, using the 12 $\mu$m flux density
measurement from WISE and the 2–10 kev flux from Swift-XRT.

Figure 1: The relationship between the observed 12$\mu$m flux and the observed
2–10 keV flux, adapted from Figure 1 in Asmus et al. (2015). The closed shapes
represent unobscured AGN (LogN${}_{\rm H}<$22) while open shapes represent
obscured AGN (LogN${}_{\rm H}>$22). The blue circles are Seyfert (Sy) 1-1.5
galaxies; green triangles are Sy 1.8-1.9 galaxies; and red squares are Sy 2
galaxies. The black stars represent confirmed CT-AGN. The obscured sources
fall to the left of the trend, signifying that the ratio between these two
quantities can be a tracer of obscuration.

####  4.2.3 Soft X-ray Hardness Ratios

Soft X-rays (0.3–10 keV) are very susceptible to changes in column density, as
evidenced in Figure 2. It can be seen that the 0.3–10 keV emission is far more
suppressed in a source with Log(NH)=24 compared to a source with Log(NH)=23.
Therefore, the ratio between the counts in different energy bands, or hardness
ratios, covering this energy band are highly dependent on column density. For
this reason, we have included two hardness ratios from the latest Swift-XRT
point source catalog, the 2SXPS (Evans et al., 2020); (M-S)/(M+S) and
(H-M)/(H+M) where S, M, and H correspond to the 0.3–1, 1–2, and 2–10 keV
bands.

Figure 2: Simulated X-ray spectra of an AGN with line-of-sight column density
Log(NH)=23 (black solid line) and Log(NH)=24 (red dotted line). The vertical
regions denoted by S (blue), M (orange), and H (green) represent the different
bands used in our two hardness ratios where they correspond to the 0.3–1, 1–2,
and 2–10 keV bands, respectively. The two spectra show extreme differences in
the soft X-rays, particularly in the 2–10 keV band. Thus, hardness ratios
targeting this band are helpful in determining the column density of AGN.

####  4.2.4 Hard X-ray Count Rates

While significantly less affected than soft X-rays, hard X-rays do display an
increased curvature with higher column densities. Koss et al. (2016) analyzed
sources with Swift-BAT data and found a correlation between the spectral
curvature and column density. Using simulated data of CT-AGN, they generated
the following equation:

$\text{SC}_{\text{BAT}}=\frac{-3.42\times A-0.82\times B+1.65\times
C+3.58\times D}{\text{Total Rate}},$ (4)

where A, B, C, and D refer to the 14–20 keV, 20–24 keV, 24-35 keV, and 35–50
keV bands, while the total rate is the 14–50 keV band. As plotted in Figure 3,
an increase in this value (calculated via Equation 4) was linked to an
increase in line-of-sight column density. Two different models that calculate
NH are plotted and all agree that the spectral curvature value increases with
NH. However, we note that this method is only valid up to NH = 4$\times$1024
cm-2.
We used this principle to improve our algorithm. Whereas Koss et al. (2016)
only included data up to 50 keV, we found our algorithm performed better when
including data up to 150 keV. Because of this, we included BAT count rates for
nine different energy bands in our algorithm: 14–20 keV, 20–24 keV, 24–34 keV,
34–45 keV, 45–60 keV, 60–85 keV, 85–110 keV, 110–150 keV, and 14–150 keV.
Including each band as a parameter accounts for both the curvature in the
spectrum while also serving as a proxy for the BAT flux. For this reason, we
elected to include every band instead of just the spectral curvature value.

Figure 3: The spectral curvature value based on different column densities,
displaying some of the configurations adapted from Figure 3 in Koss et al.
(2016). Each curve represents how the SC value changes with NH based on
different input parameters used in their MYTorus model simulations. The red
circles show the curve when the SC equation is calibrated to Swift-BAT data,
while the blue triangles show the different SC equation when calibrated to
NuSTAR data. The dashed grey line indicates the cutoff for CT-AGN determined
by Koss et al. (2016). Both lines illustrate how the increased curvature of
hard X-rays of AGN is related to an increase in column density.

##  5 Results

Figure 4 shows the X-ray-confirmed NH values for the 91 sources plotted
against the predictions by our algorithm (blue circles). We used the Spearman
rank correlation coefficient to measure the strength of the correlation
between the two sets of NH values. Our algorithm yielded a Spearman
coefficient of 0.86, signifying that it performs very well in recreating the
true NH values of these sources. Moreover, of the 31 heavily obscured (Log(NH)
$\geq$ 23) sources in our test sample, our algorithm correctly predicted 25 of
them (80%) with Log(NH) $>$ 23 and 30 (97%) with Log(NH) ¿ 22.80. We note that
we are currently unable to distinguish this obscuration as being caused by the
nuclei or the host galaxy, particularly for edge-on galaxies.
In order to determine which input parameters were most impactful in training
our algorithm, we used the percent difference of the Spearman correlation
coefficient when all parameters were used (0.86) and the coefficient when only
the parameter listed was excluded. The larger this difference, the worse our
algorithm performed without including said input parameter. Since removing one
WISE color or BAT count rate had little effect, we grouped the parameters as
such: all six WISE colors; WISE colors + the MIR-X-ray flux ratio (MIR); the
MIR-X-ray flux ratio; the two XRT hardness ratios; the two hardness ratios +
the MIR-X-ray flux ratio (Soft X-ray); and the BAT count rates. Since the MIR-
X-ray flux ratio includes both information from the infrared and the X-rays,
we included separate categories without it to determine which wavelength
influenced our algorithm the most. Figure 5 shows that the three parameters
using soft X-ray data (the two XRT hardness ratios and the MIR-X-ray flux
ratio) were the largest contributors towards our algorithm producing accurate
results.

Figure 4: The x-axis shows the “true” line-of-sight Log(NH) values, as
determined by spectral fitting. The y-axis shows the Log(NH) values predicted
by our machine learning algorithm (blue circles) and those predicted by the
Asmus et al. (2015) equation (orange stars). Our algorithm shows superior
predictive capabilities, particularly for lower levels of obscuration
(Log(N${}_{\rm H})<$ 23), where our algorithm does not incorrectly classify
unobscured sources as heavily obscured as displayed by the grey dash-dotted
line. The black dotted line represents the one-to-one ratio between the “true”
and predicted NH values. The errors from our algorithm were calculated
statistically. No errors are included on the orange points for readability
purposes. Figure 5: The percent difference between the Spearman correlation
coefficient including all parameters and the coefficient when the listed
parameter is excluded. The larger the difference, the worse the fit without
that parameter (i.e. the higher the importance of that parameter). Therefore,
the soft X-ray-related parameters are the highest contributors to the
predictive capability of the algorithm. “MIR” refers to the WISE colors and
the MIR-X-ray flux ratio. “HR” represents the two X-ray hardness ratios. “Soft
XR” refers to the two X-ray hardness ratios and the MIR-X-ray flux ratio.

###  5.1 Comparison with Previous Methods

####  5.1.1 Asmus et al. (2015)

Utilizing MIR data alongside soft X-ray data of a sample of 152 AGN, Asmus et
al. (2015) determined a relation to predict line-of-sight column densities.
The relation is as follows:

$\text{log}\left(\frac{N_{\text{H}}}{\text{cm}^{-2}}\right)=(14.37\pm
0.11)+(0.67\pm 0.11)\\\
\times\text{log}\left(\frac{F^{\text{nuc}}(12\mu\text{m})}{F^{\text{obs}}(2-10\,\text{keV})}\frac{\text{erg
s${}^{-1}$ cm${}^{-2}$}}{\text{mJy}}\right).$ (5)

Using the WISE 12 $\mu$m and XRT 2–10 keV fluxes, we plotted the NH values
predicted by the Asmus relation for our test sample of 91 sources in Figure 4
. While these results show a good trend for heavily obscured sources, below
Log(N${}_{\rm H})$ = 23, our machine learning algorithm performs far better.
This is quantified by the lower Spearman correlation coefficient of 0.65 for
the Asmus predictions and the real NH values. The lack of predictive
capability below 1023 cm-2 affects the whole range of possible NH values. This
is because, a priori, one does not know the ‘true’ NH of the source, and if
choosing a source with Log(N${}_{\rm H})<$ 23, the relation will confidently
place it as being heavily obscured. Therefore, sources with Log(N${}_{\rm
H})<$ 23 can actually have any value of ‘true’ Log(NH) between 20 and 23.

####  5.1.2 Pfeifle et al. (2022)

Pfeifle et al. (2022) improved upon the work of Asmus et al. (2015) by
creating a new relationship based on the ratio of the 2–10 keV and 12 $\mu$m
luminosities. Using 456 AGN detected in the 70-month BAT catalog (Ricci et
al., 2017) that also possess infrared data, their team created the relation
listed below:

$\text{log}\left(\frac{N_{\text{H}}}{\text{cm}^{-2}}\right)=20+(1.61^{+0.33}_{-0.31})\\\
\times\text{log}\left(\Bigg{|}\frac{\text{log}(\frac{L_{X,Obs.}}{L_{12\mu\text{m}}})+(0.34^{+0.06}_{-0.06})}{(-0.003^{+0.002}_{-0.005})}\Bigg{|}\right).$
(6)

With this relation, we have predicted the NH values for our sample of 91 test
sources as seen as the magenta triangles in Figure 6. Pfeifle et al. (2022)
claims that their method is most accurate when applied to sources with
Log(N${}_{\rm H})>$ 22.5, which is confirmed here. While their method is
accurate for heavily obscured sources, it is far less predictive than our
algorithm for AGN with Log(N${}_{\rm H})<$ 22.5. Just as with the Asmus
relation, this represents a significant drawback when selecting sources as we
do not know whether or not the ‘true’ Log(N${}_{\rm H})>$ 22.5. Overall, it
has a Spearman correlation coefficient of 0.27.

Figure 6: As in Figure 4, the x-axis shows the “true” line-of-sight Log(NH)
values determined by spectral fitting while the y-axis shows the Log(NH)
values predicted by our machine learning algorithm (blue circles) and by the
Pfeifle et al. (2022) relation (magenta triangles). The errors from our
algorithm were calculated statistically. No errors are included on the magenta
points for readability purposes.

####  5.1.3 Koss et al. (2016)

Koss et al. (2016) developed a method to identify new CT-AGN using weighted
averages of different Swift-BAT bands. It was determined that an SC${}_{BAT}>$
0.40 would identify a CT-AGN candidate. We applied this formula to our 91 test
sources and found 14 that would be considered CT. These sources are plotted as
red squares in Figure 7, overlapped on our machine learning predictions. We
note that this method does show promise, as 8 of the 14 sources are heavily
obscured, with Log(N${}_{\rm H})>$ 23\. However, 6 sources (43%) predicted as
CT have true Log(N${}_{\rm H})<$ 23, including two that are unobscured
(Log(N${}_{\rm H})<$ 22). Additionally, of the 14 predicted as CT, only 2
(14%) truly are. Our machine learning algorithm does not misclassify any
unobscured sources as CT and performs more accurately throughout all column
density ranges. Moreover, both sources predicted as CT by our algorithm are
truly Compton-thick.

Table 1 displays how many of the 91 test sources are divided in each of these
four categories: Compton-thick (Log(N${}_{\rm H})>$ 24), Compton-thin (23 $<$
Log(N${}_{\rm H})<$ 24), obscured (22 $<$ Log(N${}_{\rm H})<$ 23), and
unobscured (Log(N${}_{\rm H})<$ 22). As can be seen, our machine learning
algorithm performs the best overall, correctly classifying $\sim$75% of the
sources. This is particularly true for sources with Log(N${}_{\rm H})<$ 23, in
which our algorithm has a 73% accuracy, while the other two applicable methods
are both only 18% accurate.

Figure 7: As in Figure 4, the x-axis shows the “true” line-of-sight Log(NH)
values determined by spectral fitting while the y-axis shows the Log(NH)
values predicted by our machine learning algorithm (blue circles). The red
squares represent the 14 sources that were predicted to be CT based on the
spectral curvature method introduced in Koss et al. (2016). 6 of these sources
(43%) have true NH values $<$ 1023 cm-2. Our algorithm makes no such
misclassifications. The errors from our algorithm were calculated
statistically.

##  6 Conclusions

In this work, we present a new machine learning algorithm that predicts the
line-of-sight column density of AGN, thus enabling us to discover new CT-AGN
candidates. Using MIR data from WISE, soft X-ray data from Swift-XRT and hard
X-ray data from Swift-BAT, our machine learning algorithm has proven its
ability to accurately reproduce the NH values of our 91-source test sample,
correctly classifying 75% of sources based on their obscuration. Moreover, our
algorithm has shown a superior ability to predict the column density of AGN
with Log(NH) $<$ 22.5 when compared with previously published methods. In the
future, this algorithm will be used to: 1) identify promising CT-AGN
candidates and 2) efficiently determine NH values of large samples of sources
(like the Chandra and XMM-Newton source catalogs) in an effort to determine
the obscuration distribution of the entire AGN population across cosmic time.

Table 1: We have split the 91 test sources into four classifications based on their X-ray-measured column densities: Compton-thick (Log(N${}_{\rm H})>$ 24), Compton-thin (23 $<$ Log(N${}_{\rm H})<$ 24), obscured (22 $<$ Log(N${}_{\rm H})<$ 23), and unobscured (Log(N${}_{\rm H})<$ 22). The number of sources correctly classified for each of the four methods mentioned in this paper are shown below. Classification | Real Number | This Work | Asmus et al. (2015) | Pfeifle et al. (2022) | Koss et al. (2016)
---|---|---|---|---|---
Compton-thick | 3 | 2 | 2 | 3 | 2
Compton-thin | 28 | 22 | 25 | 13 | …
Obscured | 24 | 16 | 8 | 7 | …
Unobscured | 36 | 28 | 3 | 4 | …
Total | 91 | 68 (75%) | 38 (42%) | 27 (30%) | …

## References

  * Ajello et al. (2008) Ajello, M., Greiner, J., Sato, G., et al. 2008, The Astrophysical Journal, 689, 666
  * Alexander et al. (2003) Alexander, D. M., Bauer, F. E., Brandt, W. N., et al. 2003, The Astronomical Journal, 126, 539
  * Alexander et al. (2008) Alexander, D. M., Chary, R. R., Pope, A., et al. 2008, ApJ, 687, 835
  * Almeida & Ricci (2017) Almeida, C. & Ricci, C. 2017, Nature Astronomy, 1, 679
  * Ananna et al. (2019) Ananna, T. T., Treister, E., Urry, C. M., et al. 2019, ApJ, 871, 240
  * Andonie et al. (2022) Andonie, C., Alexander, D. M., Rosario, D., et al. 2022, MNRAS, 517, 2577
  * Asmus et al. (2015) Asmus, D., Gandhi, P., Hönig, S. F., Smette, A., & Duschl, W. J. 2015, MNRAS, 454, 766
  * Baloković et al. (2018) Baloković, M., Brightman, M., Harrison, F. A., et al. 2018, ApJ, 854, 42
  * Brandt & Yang (2021) Brandt, W. N. & Yang, G. 2021, arXiv e-prints, arXiv:2111.01156
  * Burlon et al. (2011) Burlon, D., Ajello, M., Greiner, J., et al. 2011, The Astrophysical Journal, 728, 58
  * Chainakun et al. (2022) Chainakun, P., Fongkaew, I., Hancock, S., & Young, A. J. 2022, MNRAS, 513, 648
  * Chen et al. (2021) Chen, Y., Gu, Q., Fan, J., et al. 2021, ApJ, 913, 93
  * Cutri et al. (2013) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, VizieR Online Data Catalog, II/328
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604
  * Evans et al. (2020) Evans, P. A., Page, K. L., Osborne, J. P., et al. 2020, The Astrophysical Journal Supplement Series, 247, 54
  * Ferrarese & Ford (2005) Ferrarese, L. & Ford, H. 2005, Space Sci. Rev., 116, 523
  * Ferrarese & Merritt (2000) Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9
  * Finke et al. (2021) Finke, T., Krämer, M., & Manconi, S. 2021, MNRAS, 507, 4061
  * Fiore, F. et al. (2017) Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, The Astrophysical Journal, 539, L13
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005
  * Gilli et al. (2007) Gilli, Comastri, A., & Hasinger, G. 2007, A&A, 463, 79
  * Haardt & Maraschi (1993) Haardt, F. & Maraschi, L. 1993, ApJ, 413, 507
  * Hönig (2019) Hönig, S. F. 2019, ApJ, 884, 171
  * Kilerci Eser et al. (2020) Kilerci Eser, E., Goto, T., Güver, T., Tuncer, A., & Ataş, O. H. 2020, Monthly Notices of the Royal Astronomical Society, 494, 5793
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Koss et al. (2016) Koss, M. J., Assef, R., Baloković, M., et al. 2016, ApJ, 825, 85
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, The Astronomical Journal, 115, 2285
  * Marchesi et al. (2017) Marchesi, S., Ajello, M., Comastri, A., et al. 2017, ApJ, 836, 116
  * Martín-Navarro et al. (2018) Martín-Navarro, I., Brodie, J. P., Romanowsky, A. J., Ruiz-Lara, T., & van de Ven, G. 2018, Nature, 553, 307
  * Merloni et al. (2010) Merloni, A., Bongiorno, A., Bolzonella, M., et al. 2010, The Astrophysical Journal, 708, 137
  * Merritt & Ferrarese (2001) Merritt, D. & Ferrarese, L. 2001, The Astrophysical Journal, 547, 140
  * Mizukoshi et al. (2022) Mizukoshi, S., Minezaki, T., Tsunetsugu, S., et al. 2022, MNRAS, 516, 2876
  * Murphy & Yaqoob (2009) Murphy, K. D. & Yaqoob, T. 2009, MNRAS, 397, 1549
  * Oh et al. (2018) Oh, K., Koss, M., Markwardt, C. B., et al. 2018, The Astrophysical Journal Supplement Series, 235, 4
  * Pfeifle et al. (2022) Pfeifle, R. W., Ricci, C., Boorman, P. G., et al. 2022, ApJS, 261, 3
  * Ricci et al. (2017) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, ApJS, 233, 17
  * Ricci et al. (2015) Ricci, C., Ueda, Y., Koss, M. J., et al. 2015, The Astrophysical Journal, 815, L13
  * Richstone et al. (1998) Richstone, D., Ajhar, E. A., Bender, R., et al. 1998, Nature, 385, A14
  * Silver et al. (2022) Silver, R., Torres-Albà, N., Zhao, X., et al. 2022, ApJ, 932, 43
  * Stern et al. (2005) Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, The Astrophysical Journal, 631, 163
  * Torres-Albà et al. (2021) Torres-Albà, N., Marchesi, S., Zhao, X., et al. 2021, ApJ, 922, 252
  * Treister et al. (2004) Treister, E., Urry, C. M., Chatzichristou, E., et al. 2004, The Astrophysical Journal, 616, 123
  * Treister et al. (2009) Treister, E., Urry, C. M., & Virani, S. 2009, The Astrophysical Journal, 696, 110
  * Ueda et al. (2014) Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, The Astrophysical Journal, 786, 104
  * Urry & Padovani (1995) Urry, C. M. & Padovani, P. 1995, Publications of the Astronomical Society of the Pacific, 107, 803
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Zhao et al. (2019a) Zhao, X., Marchesi, S., & Ajello, M. 2019a, ApJ, 871, 182
  * Zhao et al. (2019b) Zhao, X., Marchesi, S., Ajello, M., et al. 2019b, ApJ, 870, 60
  * Zubovas et al. (2022) Zubovas, K., Bialopetravičius, J., & Kazlauskaitė, M. 2022, MNRAS, 515, 1705