# Out-of-Distribution Detection with a Single Unconditional Diffusion Model

Alvin Heng1, Alexandre H. Thiery2, Harold Soh1,3
1Department of Computer Science, National University of Singapore
2Department of Statistics and Data Science, National University of Singapore
3Smart Systems Institute, National University of Singapore
{alvinh<EMAIL_ADDRESS>
###### Abstract

Out-of-distribution (OOD) detection is a critical task in machine learning
that seeks to identify abnormal samples. Traditionally, unsupervised methods
utilize a deep generative model for OOD detection. However, such approaches
necessitate a different model when evaluating abnormality against a new
distribution. With the emergence of foundational generative models, this paper
explores whether a single generalist model can also perform OOD detection
across diverse tasks. To that end, we introduce our method, Diffusion Paths,
(DiffPath) in this work. DiffPath proposes to utilize a single diffusion model
originally trained to perform unconditional generation for OOD detection.
Specifically, we introduce a novel technique of measuring the rate-of-change
and curvature of the diffusion paths connecting samples to the standard
normal. Extensive experiments show that with a single model, DiffPath
outperforms prior work on a variety of OOD tasks involving different
distributions. Our code is publicly available at https://github.com/clear-
nus/diffpath.

##  1 Introduction

Out-of-distribution (OOD) detection, or anomaly detection, seeks to detect
abnormal samples that are drawn far from a given distribution. This is a vital
problem as deep neural networks are known to be overconfident when making
incorrect predictions on OOD samples [1, 2], leading to potential issues in
safety-critical applications such as healthcare, finance and criminal justice
[3]. Traditionally, OOD detection using only unlabeled data relies on training
a generative model on in-distribution (ID) data. Thereafter, measures such as
model likelihood or its variants are used as an OOD detection score [4, 5, 6].
An alternative approach is to utilize the excellent sampling capabilities of
diffusion models (DMs) to reconstruct corrupted samples, and use the
reconstruction loss as an OOD measure [7, 8, 9].

However, these conventional methods necessitate separate generative models
tailored to specific ID datasets and require retraining if the ID data
evolves, such as in continual learning setups. This prompts the question: can
OOD detection be performed using a single generative model? We answer in the
affirmative and present Diffusion Paths (DiffPath) in this paper. While the
use of a single model for OOD detection has been proposed in the
discriminative setting [10], to the best of our knowledge we are the first to
explore this for generative models. We believe that the generative setting is
particularly salient in light of recent trends where large foundational
generative models achieve state-of-the-art (SOTA) results across various tasks
[11, 12].

Our method utilizes a single DM trained on a diverse dataset (e.g., ImageNet).
In a departure from prior works that utilize variants of likelihoods [6, 5, 4]
or reconstruction losses [7, 8, 9], we propose to perform OOD detection by
measuring characteristics of the forward diffusion trajectory, specifically
its rate-of-change and curvature, which can be computed from the score
predicted by the diffusion model. We provide theoretical and empirical
analyses that motivate these quantities as useful OOD detectors; their
magnitudes are similar for samples from the same distribution and different
otherwise. We summarize our contributions as follows:

Figure 1: Illustration of the diffusion paths of two different datasets
(CIFAR10 and SVHN) obtained via DDIM integration. The paths have different
first and second derivatives (rate-of-change and curvature). We propose to
measure these quantities for OOD detection.

  1. 1.

We introduce a novel approach to OOD detection by examining the rate-of-change
and curvature along the diffusion path connecting different distributions to
standard normal.

  2. 2.

Through comprehensive experiments with various OOD datasets, we show that a
single generative model outperforms baselines that necessitates separate
models for each distribution.

  3. 3.

We offer a theoretical framework demonstrating that our method characterizes
properties of the optimal transport (OT) path between the data distribution
and the standard normal.

##  2 Background

#### Score-based Diffusion Models.

Let $p_{0}({\mathbf{x}})$ denote the data distribution. We define a stochastic
differential equation (SDE), also known as the forward process, to diffuse
$p_{0}({\mathbf{x}})$ to a noise distribution $p_{T}({\mathbf{x}})$:

${\mathrm{d}}{\mathbf{x}}_{t}={\bm{f}}({\mathbf{x}}_{t},t){\mathrm{d}}t+g(t){\mathrm{d}}{\mathbf{w}}_{t},\
\ \ {\mathbf{x}}_{0}\sim p_{0}({\mathbf{x}})$ (1)

where ${\bm{f}}(\cdot,t):\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ is the drift
coefficient, $g(t)\in\mathbb{R}$ is the diffusion coefficient and
${\mathbf{w}}_{t}\in\mathbb{R}^{D}$ is the standard Wiener process (Brownian
motion). We denote $p_{t}$ as the marginal distribution of Eq. 1 at time $t$.
By starting from noise samples ${\mathbf{x}}_{T}\sim p_{T}$, new samples
${\mathbf{x}}_{0}\sim p_{0}({\mathbf{x}})$ can be sampled by simulating the
reverse SDE

${\mathrm{d}}{\mathbf{x}}_{t}=[{\bm{f}}({\mathbf{x}}_{t},t)-g(t)^{2}\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}}_{t})]{\mathrm{d}}t+g(t)d\bar{{\mathbf{w}}}_{t},\ \ \
{\mathbf{x}}_{T}\sim p_{T}({\mathbf{x}})$ (2)

where $\bar{{\mathbf{w}}}_{t}$ is the standard Wiener process when time flows
backwards from $T$ to 0, and ${\mathrm{d}}t$ is an infinitesimal negative
timestep. The diffusion process described by Eq. 1 also has an equivalent ODE
formulation, termed the probability flow (PF) ODE [13], given by

${\mathrm{d}}{\mathbf{x}}_{t}=\left[{\bm{f}}({\mathbf{x}}_{t},t)-\frac{1}{2}g(t)^{2}\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}}_{t})\right]{\mathrm{d}}t.$ (3)

The ODE and SDE formulations are equivalent in the sense that trajectories
under both processes share the same marginal distribution
$p_{t}({\mathbf{x}}_{t})$. Hence, given an estimate of the score function
$s_{\theta}({\mathbf{x}}_{t},t)\approx\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}}_{t})$, which can be obtained using score-matching
approaches [14, 13], one can sample from the diffusion model by solving the
reverse SDE or integrating the PF ODE backwards in time.

In this work, we focus on the variance-preserving formulation used in DDPM
[15], which is given by an Ornstein-Uhlenbeck forward process

${\mathrm{d}}{\mathbf{x}}_{t}=-\frac{1}{2}\beta_{t}{\mathbf{x}}_{t}{\mathrm{d}}t+\sqrt{\beta_{t}}{\mathrm{d}}{\mathbf{w}}_{t},\
\ \ {\mathbf{x}}_{0}\sim p_{0}({\mathbf{x}})$ (4)

where $\beta_{t}$ are time-dependent constants. Under Eq. 4, diffused samples
${\mathbf{x}}_{t}$ can be sampled analytically via
$p_{t}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=\mathcal{N}({\mathbf{x}}_{t};\sqrt{\bar{\alpha}_{t}}{\mathbf{x}}_{0},\sigma_{t}^{2}{\mathbf{I}})$,
where $\beta_{t}=-\frac{{\mathrm{d}}}{{\mathrm{d}}t}\log\bar{\alpha}_{t}$ and
$\sigma_{t}^{2}=1-\bar{\alpha}_{t}$. The score estimator,
$\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\approx-\sigma_{t}\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}}_{t})$, can be trained via the following objective

$\min_{\theta}\mathbb{E}_{t\sim\mathcal{U}[0,1]{\mathbf{x}}_{0}\sim
p_{0}({\mathbf{x}}_{0}){\mathbf{x}}_{t}\sim
p_{t}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})}\left[\left\lVert\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)-\bm{\epsilon}\right\rVert^{2}_{2}\right],$
(5)

where $\bm{\epsilon}=-\sigma_{t}\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}}_{t}|{\mathbf{x}}_{0})=({\mathbf{x}}_{t}-\sqrt{\bar{\alpha}}_{t}{\mathbf{x}}_{0})/\sigma_{t}$.

#### Unsupervised OOD Detection.

We formalize the task of unsupervised OOD detection here. Given a distribution
of interest $p({\mathbf{x}})$, the goal of OOD detection is to construct a
scoring function which outputs a quantity
$S_{\theta}({\mathbf{x}})\in\mathbb{R}$ that identifies if a given test point
${\mathbf{x}}_{\textrm{test}}$ is from $p({\mathbf{x}})$. In this work, a
higher value of $S_{\theta}({\mathbf{x}}_{\textrm{test}})$ indicates that the
sample is more likely to be drawn from $p({\mathbf{x}})$. We will use the
notation “A vs B” to denote the task of distinguishing samples between A and
B, where A is the in-distribution (ID) and B is the out-of-distribution (OOD)
dataset. In unsupervised OOD detection, one must construct the function
$S_{\theta}$ using only knowledge of A.

##  3 Diffusion Models for OOD Detection

An overview of our method, DiffPath, is illustrated in Fig. 1. DiffPath is
based on the insight that the rate-of-change and curvature of the diffusion
path connecting samples to standard normal differ between distributions,
making them effective indicators for OOD detection. This section outlines the
methodology behind DiffPath We begin in Sec. 3.1, where we provide evidence
that likelihoods of a diffusion model are insufficient for OOD detection.
Next, we show that the score function is a measure of the rate-of-change, as
well as motivate the use of a single generative model in Sec. 3.2. Following
which, we motivate the curvature as the derivative of the score in Sec. 3.3.
We consider the curvature statistic as one variation of our method and
abbreviate it as DiffPath-1D. In Sec. 3.4, we contextualize our method in
terms of the optimal transport path between samples and standard normal, and
finally propose a higher-order, hybrid variation called DiffPath-6D in Sec.
3.5, which incorporates both the rate-of-change and curvature quantities.

###  3.1 Likelihoods Do Not Work

When leveraging a likelihood-based generative model for OOD detection, the
most natural statistic to consider is the likelihood itself. As DMs are
trained to maximize the evidence lower bound (ELBO) of the data, one would
expect that in-distribution samples have higher ELBO under the model compared
to out-of-distribution samples. However, prior works [2, 16] have shown that
the opposite behavior was observed in deep generative models, such as
normalizing flows, where the model assigned higher likelihoods to OOD samples.

In Fig. 2, we plot the distributions of the negative ELBO (denoted NLL) of
various training sets for a DM trained on CIFAR10. Our results corroborate
earlier findings that likelihoods are not good OOD detectors; the NLL of
CIFAR10 samples are higher than SVHN samples, meaning ID samples have lower
likelihoods than OOD samples. Furthermore, the NLL of CIFAR10 overlaps
significantly with that of CelebA samples. The poor AUROC score in Table 1
quantitatively demonstrates the inability of likelihoods to distinguish
between ID and OOD samples. This motivates us to search for better statistics
that we can extract from DMs for OOD detection.

###  3.2 Scores as OOD Statistic

Figure 2: Histograms of various statistics of the respective training sets.
The NLL is calculated using a diffusion model trained on CIFAR10, while the
other two statistics are calculated with a model trained on ImageNet.

Table 1: AUROC of statistics shown in Fig. 2.

| C10 vs
---|---
Method | SVHN | CelebA
NLL | 0.091 | 0.573
$\sqrt{\sum_{t}\left\lVert\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\right\rVert_{2}^{2}}$ | 0.856 | 0.638
$\sqrt{\sum_{t}\left\lVert\partial_{t}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\right\rVert_{2}^{2}}$ | 0.965 | 0.998

#### Scores as KL Divergence Proxy.

We start by rewriting the PF ODE, Eq. 3, in the following form:

$\frac{{\mathrm{d}}{\mathbf{x}}_{t}}{{\mathrm{d}}t}={\bm{f}}({\mathbf{x}}_{t},t)+\frac{g(t)^{2}}{2\sigma_{t}}\bm{\epsilon}_{p}({\mathbf{x}}_{t},t)$
(6)

where we have parameterized the score as
$\bm{\epsilon}_{p}({\mathbf{x}}_{t},t)=-\sigma_{t}\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}})$.

######  Theorem 1.

Denote $\phi_{t}$ and $\psi_{t}$ as the marginals from evolving two distinct
distributions $\phi_{0}$ and $\psi_{0}$ via their respective probability flow
ODEs (Eq. 6) forward in time. We consider the case with the same forward
process, i.e., the two PF ODEs have the same
${\bm{f}}({\mathbf{x}}_{t},t),g(t)\textrm{ and }\sigma_{t}$. Under some
regularity conditions stated in Appendix A.1,

$D_{\mathrm{KL}}(\phi_{0}\|\psi_{0})=\frac{1}{2}\int_{0}^{T}\mathbb{E}_{{\mathbf{x}}\sim\phi_{t}}\frac{g(t)^{2}}{\sigma_{t}}\left\lVert\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t)-\bm{\epsilon}_{\psi}({\mathbf{x}}_{t},t)\right\rVert^{2}_{2}{\mathrm{d}}t+D_{\mathrm{KL}}(\phi_{T}\|\psi_{T}).$

The term $D_{\mathrm{KL}}(\phi_{T}\|\psi_{T})$ vanishes as
$\phi_{T}=\psi_{T}=\mathcal{N}(\bm{0},{\mathbf{I}})$ by construction, assuming
the true scores are available. In practice, one does not have access to the
true score, and we rely on a score estimator $\bm{\epsilon}_{\theta}$ obtained
via score matching approaches. Theorem 1 suggests that the scores of the
marginal distributions along the ODE path serve as a proxy for the KL
divergence: as $D_{\mathrm{KL}}(\phi_{0}\|\psi_{0})$ increases, so should the
difference in the norms of their scores. Another interpretation is that this
difference,
$\mathbb{E}[\left\lVert\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t)-\bm{\epsilon}_{\psi}({\mathbf{x}}_{t},t)\right\rVert^{2}_{2}]$,
is a measure of the Fisher divergence between the two distributions, which
forms the foundation for score matching [17]. Thus, this motivates the use of
the norm of the scores as an OOD statistic in distinguishing two distributions
apart.

However, Theorem 1 is not immediately useful as it requires a priori knowledge
of both distributions, whereas in unsupervised OOD detection only knowledge of
the in-distribution is available. Interestingly, we empirically observe that
it is possible to approximate the forward probability flow ODE for different
distributions using a single diffusion model when it is trained on a large,
diverse dataset $p_{0}({\mathbf{x}})$. Recall that as the PF ODE has the same
marginal as the forward SDE, if the score estimate $\bm{\epsilon}_{\theta}$
has converged to the true score, then forward integration of a sample
${\mathbf{x}}_{0}$ using Eq. 6 should bring the sample to approximately
standard normal, ${\mathbf{x}}_{T}\sim\mathcal{N}(\bm{0},{\mathbf{I}})$.

Specifically, we consider the following parameterization [18] of the PF ODE

${\mathrm{d}}\bar{{\mathbf{x}}}_{t}=\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t){\mathrm{d}}\gamma_{t}$
(7)

where $\gamma_{t}=\sqrt{\frac{1-\bar{\alpha}_{t}^{2}}{\bar{\alpha}_{t}}}$ and
$\bar{{\mathbf{x}}}_{t}={\mathbf{x}}_{t}\sqrt{1+\gamma_{t}^{2}}$. Let
$\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)=-\sigma_{t}\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}})$ be a score model trained on $p_{0}({\mathbf{x}})$. It is
known that the DDIM sampler [19] is Euler’s method applied to Eq. 7. In Fig.
3, we integrate Eq. 7 forward in time using DDIM for samples from
distributions unseen by the model during training. We observe the surprising
fact that when $p_{0}({\mathbf{x}})$ is ImageNet, the model is able to bring
the samples approximately to the standard normal. However, this is not true
for the model trained on SVHN, suggesting that $p_{0}$ needs to be a large and
diverse dataset. We ablate the choice of $p_{0}$ in Sec. 5.2.

Figure 3: Illustration of the forward integration of Eq. 7 on samples from
CIFAR10, SVHN and CelebA. The ImageNet model is able bring the samples to
approximately standard normal despite never being trained explicitly on them.
The SVHN model is only able to bring the SVHN sample to standard normal. This
highlights the necessity for the diffusion model to be trained on a large,
diverse dataset so that it can estimate the scores for a variety of data
accurately.

This motivates $\bm{\epsilon}_{\theta}$ as a replacement for arbitrary
$\bm{\epsilon}_{\phi}$ when integrating Eq. 7 forward with samples from
$\phi_{0}$. In Fig. 2, we see that the distributions of
$\sqrt{\sum_{t}\left\lVert\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\right\rVert_{2}^{2}}$
, the square root of the sum of $L^{2}$ norms of scores over time, applied to
the three datasets using a single model trained on ImageNet are better
separated than the likelihoods for CIFAR10 vs SVHN. We leverage this insight
to design a simple scheme for OOD detection: fit a Kernel Density Estimator
(KDE) to
$\sqrt{\sum_{t}\left\lVert\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\right\rVert_{2}^{2}}$
of the training set for a given distribution, then use the KDE likelihoods of
a test sample as the OOD score $S_{\theta}$. We provide pseudocode in
Algorithm 1. The AUROC scores in Table 1 show a large improvement in CIFAR10
vs SVHN scores over likelihoods, although CIFAR10 vs CelebA is still lacking.

Algorithm 1 OOD detection with DiffPath

1:Trained DM $\bm{\epsilon}_{\theta}$, ID train set
${\mathbf{X}}_{\textrm{train}}$, test samples ${\mathbf{X}}_{\textrm{test}}$,
empty lists $L_{\textrm{train}}$ and $L_{\textrm{test}}$

2:OOD scores of test samples $S_{\theta}({\mathbf{X}}_{\textrm{test}})$

3:for ${\mathbf{x}}_{0}$ in ${\mathbf{X}}_{\textrm{train}}$ do

4:
$\\{\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\\}_{t=0}^{T}\leftarrow\texttt{DDIMInversion}({\mathbf{x}}_{0},\bm{\epsilon}_{\theta})$
$\triangleright$ Integrate Eq. 7 from $t=0$ to $T$

5: Calculate OOD statistic using
$\\{\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\\}_{t=0}^{T}$

6: Append statistic to $L_{\textrm{train}}$

7:end for

8:$p_{\textrm{train}}(\cdot)\leftarrow$ fit density estimate to
$L_{\textrm{train}}$ $\triangleright$ e.g., KDE, GMM

9:$L_{\textrm{test}}\leftarrow$ Repeat lines $1-5$ with
${\mathbf{X}}_{\textrm{test}}$

10:return $p_{\textrm{train}}(l)$ for every $l$ in $L_{\textrm{test}}$

#### Score as First-Order Taylor Expansion.

We provide a second interpretation of the score, and subsequently motivate a
new statistic that can be used for OOD detection. Recall that the numerical
DDIM solver is the first-order Euler’s method applied to Eq. 7. In general, we
can expand the ODE to higher-order terms using the truncated Taylor method
[20, 18]:

$\displaystyle\begin{split}\bar{{\mathbf{x}}}_{t_{n+1}}&=\bar{{\mathbf{x}}}_{t_{n}}+h_{n}\frac{{\mathrm{d}}\bar{{\mathbf{x}}}_{t_{n}}}{{\mathrm{d}}\gamma_{t}}+\frac{1}{2}h^{2}_{n}\frac{{\mathrm{d}}^{2}\bar{{\mathbf{x}}}_{t_{n}}}{{\mathrm{d}}\gamma_{t}^{2}}\big{|}_{(\bar{{\mathbf{x}}}_{t_{n}},t_{n})}+\dots\\\
&=\bar{{\mathbf{x}}}_{t_{n}}+h_{n}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t_{n}},t_{n})+\frac{1}{2}h^{2}_{n}\frac{{\mathrm{d}}\bm{\epsilon}_{\theta}}{{\mathrm{d}}\gamma_{t}}\big{|}_{(\bar{{\mathbf{x}}}_{t_{n}},t_{n})}+\dots\end{split}$
(8)

where $h_{n}=\gamma_{t_{n+1}}-\gamma_{t_{n}}$. The norm of the first-order
score, $\left\lVert\bm{\epsilon}_{\theta}\right\rVert_{2}$, therefore measures
the rate-of-change of the ODE integration path. Intuitively, the ODE
integration path necessary to bring different distributions to the standard
normal in finite time differs (c.f. the PF ODE path is also the optimal
transport path, see Sec. 3.4) , hence the rate-of-change differs as well.

###  3.3 Second-Order Taylor Expansion (DiffPath-1D)

Based on the preceding discussion, it is natural to consider if higher-order
terms in the ODE Taylor expansion can also serve as an effective OOD
statistic. We answer in the affirmative by considering the second order term,
$\frac{{\mathrm{d}}\bm{\epsilon}_{\theta}}{{\mathrm{d}}\gamma_{t}}$.
Intuitively, the second-order term measures the curvature of the ODE
integration path. We expand the derivative as follows [18]:

$\displaystyle\begin{split}\frac{{\mathrm{d}}\bm{\epsilon}_{\theta}}{{\mathrm{d}}\gamma_{t}}&=\frac{\partial\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\partial{\mathbf{x}}_{t}}\frac{{\mathrm{d}}{\mathbf{x}}_{t}}{{\mathrm{d}}\gamma_{t}}+\frac{\partial\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\partial
t}\frac{{\mathrm{d}}t}{{\mathrm{d}}\gamma_{t}}\\\
&=\frac{1}{\sqrt{\gamma_{t}^{2}+1}}\underbrace{\frac{\partial\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\partial{\mathbf{x}}_{t}}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}_{\textrm{JVP}}-\frac{\gamma_{t}}{1+\gamma_{t}^{2}}\underbrace{\frac{\partial\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\partial{\mathbf{x}}_{t}}{\mathbf{x}}_{t}}_{\textrm{JVP}}+\frac{\partial\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\partial
t}\frac{{\mathrm{d}}t}{{\mathrm{d}}\gamma_{t}}.\end{split}$ (9)

We see that the derivative contains two Jacobian-Vector Products (JVP) and a
simple time derivative term. In principle, all three terms can be computed
using automatic differentiation. However, this makes inference twice as costly
due to the need for an additional backward pass after every forward pass of
the network, and significantly more memory-intensive due to storage of the
full computation graph. Fortunately, the time derivative term can be computed
using simple finite difference:

$\frac{\partial\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\partial
t}\approx\frac{\bm{\epsilon}_{\theta}({\mathbf{x}}_{t+\Delta t},t+\Delta
t)-\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)}{\Delta t}$ (10)

where the pairs $({\mathbf{x}}_{t},{\mathbf{x}}_{t+\Delta t})$ are obtained
during standard DDIM integration. Thus, no additional costs associated with
gradient computation are incurred.

Surprisingly, we observe that for high-dimensional distributions such as
images that we consider in this work, the time derivative in Eq. 10 alone
provides an accurate enough estimate for OOD detection. Using the same
ImageNet model as in Sec. 3.2, we see from Fig. 2 that the distributions of
$\sqrt{\sum_{t}\left\lVert\partial_{t}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\right\rVert_{2}^{2}}$
for different datasets are more spread out. This reflects quantitatively in
the AUROC scores in Table 1, where the second-order statistic improves
significantly upon the first-order scores, in particular for CIFAR10 vs CelebA
with near perfect performance.

We thus consider the second-order statistic alone,
$\sqrt{\sum_{t}\left\lVert\partial_{t}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\right\rVert_{2}^{2}}$
, as one of our main statistic for OOD detection. As it is a one-dimensional
quantity, we abbreviate it as DiffPath-1D.

###  3.4 Connections to Optimal Transport

Recent works have viewed DDIM integration as an encoder that maps the data
distribution to standard normal [21, 22]. They prove that this map is the
optimal transport (OT) path if the data distribution is Gaussian, while
providing numerical results suggesting likewise for high-dimensional data like
images. As a result, we can view the OOD statistics proposed in Sec. 3.2 and
Sec. 3.3 as characterizing different derivatives of the OT path: the score
$\left\lVert\bm{\epsilon}_{\theta}\right\rVert_{2}$ represents the rate-of-
change of the path, while the time derivative
$\left\lVert\partial_{t}\bm{\epsilon}_{\theta}\right\rVert_{2}$ represents its
curvature.

To justify our proposition that derivatives of OT paths serve as meaningful
OOD statistics, we consider the following toy example [21] of distinguishing
two multivariate Gaussians (detailed derivation in Appendix A.2). Let the
distributions be
$p^{i}_{0}({\mathbf{x}})\sim\mathcal{N}({\mathbf{a}}_{i},{\mathbf{I}}),i\in\\{0,1\\}$,
where ${\mathbf{a}}_{i}\in\mathbb{R}^{d}\textrm{ and
}{\mathbf{I}}\in\mathbb{R}^{d\times d}$. The marginal densities along the
forward diffusion can be computed exactly using the transition formulas for
SDEs [23] and is given by
$p^{i}_{t}({\mathbf{x}})\sim\mathcal{N}({\mathbf{a}}_{i}e^{-t},{\mathbf{I}})$,
with PF ODE
$\frac{{\mathrm{d}}{\mathbf{x}}_{i}}{{\mathrm{d}}t}=-{\mathbf{a}}_{i}e^{-t}$.
This path corresponds exactly to the OT map between $p^{i}_{0}$ and
$\mathcal{N}(\bm{0},{\mathbf{I}})$. In this case, the corresponding first and
second-order OOD statistics are equal and given by
$\left\lVert\frac{{\mathrm{d}}{\mathbf{x}}_{i}}{{\mathrm{d}}t}\right\rVert_{2}=\left\lVert\frac{{\mathrm{d}}^{2}{\mathbf{x}}_{i}}{{\mathrm{d}}t^{2}}\right\rVert_{2}=\left\lVert{\mathbf{a}}_{i}e^{-t}\right\rVert_{2}$.
Crucially, they are proportional to
$\left\lVert{\mathbf{a}}_{i}\right\rVert_{2}$, meaning that as the two
distributions move farther apart (i.e., as
$\left\lVert{\mathbf{a}}_{0}-{\mathbf{a}}_{1}\right\rVert_{2}$ increases), so
should the $L^{2}$ norms of the OOD statistics, thereby increasing their
ability to distinguish samples between the two.

###  3.5 Higher-dimensional Statistic (DiffPath-6D)

Table 2: AUROC of DiffPath 1D vs 6D. Method |  | C10 vs
---
neg. C10
1D | 0.500
6D | 0.994

Owing to its simplicity, the one-dimensional statistic proposed in Sec. 3.3
may suffer from edge cases or perform suboptimally as information is condensed
to a single scalar. For instance, given an image ${\mathbf{x}}_{0}$ with
pixels normalized to the range $[-1,1]$, one such edge case is distinguishing
${\mathbf{x}}_{0}$ from itself with the sign of the pixels flipped,
$-{\mathbf{x}}_{0}$. The two samples will produce symmetric OT paths differing
only by a negative sign, resulting in identical statistics after taking the
$L^{2}$ norm. We can see this from Table 2 where DiffPath-1D fails to
distinguish CIFAR10 samples from itself with signs flipped, which we call
negative CIFAR10. As such, we are motivated to propose a higher-dimensional
statistic that does not utilize the standard form of the $L^{p}$ norm:
$\left\lVert{\mathbf{x}}\right\rVert_{p}=\sum_{i}|{\mathbf{x}}_{i}|^{p}$. We
define a new scalar quantity,
$\langle{\mathbf{x}}\rangle_{p}=\sum_{i}({\mathbf{x}}_{i})^{p}$, which omits
the absolute value so that sign information is retained, and propose a new
six-dimensional statistic we dub DiffPath-6D:

$\Big{[}\sum_{t}\langle\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\rangle_{1},\sum_{t}\langle\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\rangle_{2},\sum_{t}\langle\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\rangle_{3},\sum_{t}\langle\partial_{t}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\rangle_{1},\sum_{t}\langle\partial_{t}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\rangle_{2},\sum_{t}\langle\partial_{t}\bm{\epsilon}_{\theta}({\mathbf{x}}_{t},t)\rangle_{3}\Big{]}^{\top}$

which concatenates scalars based on the first, second and third powers of
$\bm{\epsilon}_{\theta}$ and $\partial_{t}\bm{\epsilon}_{\theta}$ into a
vector. From Table 2, DiffPath-6D is able to distinguish CIFAR10 from neg.
CIFAR10 near perfectly. We validate both DiffPath-1D and DiffPath-6D on a
wider suite of experiments in Sec. 5.

##  4 Related Works

Modern OOD detection can be divided roughly into three categories: feature-
based, likelihood-based, reconstruction-based. Feature-based methods extract
features from ID samples and fit a likelihood or distance function as an OOD
detector. For instance, one can obtain the latent representations of a test
sample using an autoencoder and measure its Mahalanobis distance to the
representations of ID samples [24]. Distances between features derived from
self-supervised learning models are also utilized in similar contexts [25,
26]. Similar to our work, Xiao et al. [10] showed that one can perform OOD
detection using features from a single discriminative model trained on
ImageNet.

Likelihood-based approaches leverage a generative model trained on ID samples.
These methods typically employ variants of the log-likelihood of a test sample
under the model as the OOD detection score. Nalisnick et al. [2] first pointed
out that deep generative models might erroneously assign higher likelihoods to
OOD samples. Several explanations have been proposed, such as the input
complexity [16] and typicality [27] of samples. As a result, just as we show
in Sec. 3.1, vanilla likelihoods are rarely used. Instead, variants derived
from the log-likelihood have been proposed, such as likelihood ratios [4],
ensembles of the likelihood [5], density of states [6], energy-based models
[28] and typicality tests [27]. A similar line of work to ours is MSMA [29],
which uses the score function over discrete noise levels for OOD detection.
One can view MSMA as a specific case of DiffPath, where the former only
utilizes the first-order statistic, while we generalize to higher-order terms.
MSMA proposes to measure the score at various noise levels, while our method
sums over the entire diffusion path. It is worth emphasizing that MSMA
requires different models for different ID tasks, unlike our single model
setup.

Reconstruction-based approaches evaluate how well a generative model, trained
on ID data, can reconstruct a test sample. Earlier works utilize autoencoders
[30] and GANs [31] for reconstruction. Recent works have adapted unconditional
DMs to this approach due to its impressive sample quality. A test sample is
first artificially corrupted before being reconstructed using the DM’s
sampling process. DDPM-OOD [7] noises a sample using the forward process and
evaluates the perceptual quality [32] of the reconstructed sample. Projection
Regret [9] adopts a similar approach, but uses a Consistency Model [33] and
introduces an additional projection regret score. LMD [8] corrupts the image
by masking and reconstructs the sample via inpainting. Evidently, DiffPath
differs from these diffusion approaches as we do not utilize reconstructions.
We also stress again that these baselines require different models for
different ID tasks.

##  5 Experiments

Table 3: AUROC scores for various in-distribution vs out-of-distribution
tasks. Higher is better. Bold and underline denotes the best and second best
result respectively. We also show the number of function evaluations (NFE) for
diffusion methods, where lower is better.

| C10 vs | SVHN vs | CelebA vs |  |
---|---|---|---|---|---
Method | SVHN | CelebA | C100 | Textures | C10 | CelebA | C100 | Textures | C10 | SVHN | C100 | Textures | Average | NFE
IC | 0.950 | 0.863 |  0.736 | - | - | - | - | - | - | - | - | - | - | -
IGEBM | 0.630 | 0.700 | 0.500 | 0.480 | - | - | - | - | - | - | - | - | - | -
VAEBM | 0.830 | 0.770 | 0.620 | - | - | - | - | - | - | - | - | - | - | -
Improved CD | 0.910 | - | 0.830 | 0.880 | - | - | - | - | - | - | - | - | - | -
DoS | 0.955 | 0.995 | 0.571 | - | 0.962 | 1.00 | 0.965 | - | 0.949 |  0.997 | 0.956 | - | 0.928 | -
WAIC111Results obtained from Morningstar et al. [6]. | 0.143 | 0.928 | 0.532 | - | 0.802 | 0.991 | 0.831 | - | 0.507 | 0.139 | 0.535 | - | 0.601 | -
TT111Results obtained from Morningstar et al. [6]. | 0.870 | 0.848 | 0.548 | - | 0.970 | 1.00 | 0.965 | - | 0.634 | 0.982 | 0.671 | - | 0.832 | -
LR111Results obtained from Morningstar et al. [6]. | 0.064 | 0.914 | 0.520 | - | 0.819 | 0.912 | 0.779 | - | 0.323 | 0.028 | 0.357 | - | 0.524 | -
Diffusion-based
NLL | 0.091 | 0.574 | 0.521 | 0.609 | 0.990 |  0.999 | 0.992 | 0.983 | 0.814 | 0.105 | 0.786 | 0.809 | 0.689 | 1000
IC | 0.921 | 0.516 | 0.519 | 0.553 | 0.080 | 0.028 | 0.100 | 0.174 | 0.485 | 0.972 | 0.510 | 0.559 | 0.451 | 1000
MSMA | 0.957 | 0.675 | 0.615 | 0.967 | 0.976 | 0.995 |  0.980 | 0.996 | 0.871 |  0.997 | 0.886 | 0.967 | 0.907 | 10
DDPM-OOD | 0.390 | 0.659 | 0.536 | 0.598 | 0.951 | 0.986 | 0.945 | 0.910 | 0.795 | 0.636 | 0.778 | 0.773 | 0.746 | 350
LMD | 0.992 | 0.557 | 0.604 | 0.667 | 0.919 | 0.890 | 0.881 | 0.914 | 0.989 | 1.00 | 0.979 |  0.972 | 0.865 | $10^{4}$
Ours
DiffPath-1D | 0.965 |  0.998 | 0.551 | 0.966 | 0.971 | 1.00 | 0.972 | 0.992 |  0.998 | 1.00 | 0.994 | 0.827 | 0.936 |  50
DiffPath-1D |  0.966 | 0.999 | 0.552 | 0.975 | 0.971 | 1.00 | 0.973 | 0.988 | 0.999 | 1.00 | 0.998 | 0.851 | 0.939 | 100
DiffPath-6D | 0.936 | 0.997 | 0.596 |  0.997 | 0.921 |  0.999 | 0.932 |  0.999 | 0.997 | 1.00 |  0.995 | 0.935 |  0.942 |  50
DiffPath-6D | 0.920 | 0.999 | 0.593 | 0.999 | 0.924 |  0.999 | 0.936 | 1.00 | 0.999 | 1.00 | 0.998 | 0.946 | 0.943 | 100

#### Datasets.

We evaluate pairwise OOD detection using CIFAR10 (C10), SVHN, and CelebA as ID
datasets, and CIFAR100 (C100) and Textures as additional OOD datasets. We
employ an unconditional diffusion model trained on ImageNet at $64\times 64$
resolution, using the Improved DDPM’s [34] open-source
checkpoint222https://github.com/openai/improved-diffusion. All datasets are
resized to $64\times 64$ when evaluating with DiffPath (see Sec. 5.2 for
ablations on resizing techniques).

#### Methodology and Baselines.

Our methodology features two variants of our model, DiffPath-1D using KDE and
DiffPath-6D using a Gaussian Mixture Model for OOD scoring, as outlined in
Sec. 3. We compare against a variety of generative baselines, including
Energy-based Model (EBM) such as IGEBM [28], VAEBM [35] and Improved CD [36],
as well as Input Complexity (IC) [16], Density of States (DOS) [6], Watanabe-
Akaike Information Criterion (WAIC) [5], Typicality Test (TT) [27] and
Likelihood Ratio (LR) [4] applied to the Glow [37] model. Additionally, we
compare against diffusion baselines such as vanilla NLL and IC based on the
DM’s likelihoods and re-implementations of DDPM-OOD [7], LMD [8], and MSMA
[29] based on open-source code for full comparisons.

Table 4: Ablation results when we vary $p_{0}({\mathbf{x}})$, the distribution
the single DM is trained on. We use DiffPath-1D with 50 NFEs.  Bold and
underline denotes the best and second best result.

| C10 vs | SVHN vs | CelebA vs |
---|---|---|---|---
$p_{0}({\mathbf{x}})$ | SVHN | CelebA | C100 | Textures | C10 | CelebA | C100 | Textures | C10 | SVHN | C100 | Textures | Average
C10 | 0.862 | 0.359 |  0.538 |  0.724 | 0.788 | 0.807 | 0.803 | 0.862 | 0.658 | 0.927 | 0.704 |  0.839 | 0.739
SVHN | 0.975 | 0.394 | 0.512 | 0.698 | 0.985 |  0.997 | 0.984 |  0.960 | 0.625 |  0.996 | 0.621 | 0.779 |  0.794
CelebA | 0.915 |  0.55 | 0.519 | 0.693 | 0.881 | 0.779 | 0.872 | 0.840 |  0.772 | 0.902 |  0.789 | 0.850 | 0.780
ImageNet |  0.965 | 0.998 | 0.551 | 0.966 |  0.971 | 1.00 |  0.972 | 0.992 | 0.998 | 1.00 | 0.994 | 0.827 | 0.936

###  5.1 Main Results

Table 3 summarizes our main results. Both variants of DiffPath, with 50 and
100 NFEs, obtain the top four best average results across all 12 tasks, with
DiffPath-6D at 100 NFEs performing the best. Interestingly, dropping to 50
NFEs only reduces performance marginally while improving inference speed by
$2\times$. Specifically, DiffPath-6D at 100 NFEs achieves the best result in
six out of twelve tasks, and second best result on one of them, while at 50
NFEs achieves second best result in four out of twelve tasks. Despite using
only a scalar statistic for OOD detection, DiffPath-1D at 100 NFEs still
obtains the best result in five out of twelve tasks. Thus, our results
validate the efficiency and versatility of DiffPath on diverse OOD tasks.

All variants of DiffPath struggle with the C10 vs C100 task due to
similarities between these datasets, scoring below 0.7 AUROC. However, this
task is also challenging for all baselines, with only Improved CD and IC
achieving AUROC scores above 0.7. Interestingly, DiffPath’s performance is
strongest when either ID or OOD dataset is CelebA. For instance, we achieve
perfect score when distinguishing against SVHN across both variants, and near
perfect score (greater than $0.99$) against C10/C100. This corroborates the
result in Fig. 2, which shows that the histogram of CelebA is well separated
from the other datasets.

Interestingly, DiffPath maintains consistent performance when swapping ID and
OOD datasets. As a single model is used, the statistics extracted for a
dataset is identical (up to differences between train test split) regardless
of whether it is ID or OOD. Other methods do not exhibit such behavior as they
require different generative models for different tasks. For example, in the
C10 vs CelebA/CelebA vs C10 tasks, DiffPath-1D achieves identical scores of
0.999/0.999, while MSMA and LMD achieves scores of 0.675/0.871 and 0.557/0.989
respectively. We believe this symmetry is desirable as it is similar to binary
classification, where performance is unchanged when labels are swapped.

DiffPath requires fewer NFEs per sample than most diffusion models, making it
faster during inference. Only MSMA uses fewer NFEs, but it performs worse and
needs different models for each ID task. Reconstruction-based methods like
DDPM-OOD and LMD require many reconstructions per sample, leading to high NFE
requirements. In particular, LMD requires a hefty $10^{4}$ NFEs as the
original implementation requires 10 reconstructions per sample, with each
utilizing 1000 NFEs. As for NLL and IC, both underperform our method despite
our experiments showing that roughly 1000 NFEs are required to obtain an
accurate estimate of the likelihood.

Figure 4: Plots of AUROC performance for various ablation experiments. Average
performance over all tasks when the number of DDIM steps is varied (left) and
when the image resizing technique is varied (right). We use DiffPath-6D with
50 NFEs for both. Full results in Sec. B.1 of the appendix.

###  5.2 Ablations

#### Choice of Diffusion Training Set.

Here we quantitatively justify our claim in Sec. 3.2 that the data
distribution used to train the DM, $p_{0}({\mathbf{x}})$, has to be large and
diverse for OOD detection. We report ablation results in Table. 4 using DMs
trained on CIFAR10, SVHN, CelebA and ImageNet. The results show that the DM
trained on ImageNet significantly outperforms the rest, achieving the best
result across 8 out of 12 tasks, second best on 3 tasks and the best average
result by a large margin. This is despite the other variants being trained on
samples seen in the experiments, while ImageNet was not. This suggests that it
is unnecessary for generative models to be trained on in-distribution data for
OOD detection; rather, the model should be trained to learn diverse and
general representations.

#### Number of DDIM Steps.

We investigate how the performance of our method varies with the number of
DDIM steps used in Fig. 4. Interestingly, the best results are obtained at 25,
50 and 100 DDIM steps, with near identical performance. As the experiments
with 50 and 100 steps are more consistent throughout tasks (see Sec. B.1 in
appendix for full ablation results), we report results for 50 and 100 NFEs in
our main experiments in Table. 3. We hypothesize that the slight decrease in
performance at larger DDIM steps is due to the use of the finite difference in
Eq. 10, where the difference in consecutive predictions
$\bm{\epsilon}(\cdot,t+\Delta t)-\bm{\epsilon}(\cdot,t)$ shrinks rapidly. We
leave the investigation of this phenomena in greater detail to future work.

#### Image Resizing Technique.

Since we utilize a single DM trained on ImageNet at $64\times 64$ resolution,
samples from image distributions of different resolutions have to be resized
appropriately. As such, we ablate for resizing techniques that are implemented
in the PyTorch library. The results in Fig. 4 show that DiffPath is relatively
robust, with average AUROC approximately 0.9 or higher across all resizing
techniques considered. However, bicubic and bilinear interpolation outperform
the rest, with the latter obtaining the best performance. Thus, we use
bilinear interpolation for the experiments described in this work.

##  6 Conclusion

In this work, we proposed Diffusion Paths (DiffPath), a method of OOD
detection using a single diffusion model by characterizing properties of the
forward diffusion path. In light of the growing popularity of foundational
generative models, our work demonstrates that a generalist model–one that has
learnt diverse representations from data–can also be applied to out-of-
distribution detection. There are several interesting future directions that
arise from our work; for instance, applying DiffPath to other modalities,
given that diffusion modelling has been extended to domains such as video,
audio, language and time series, as well as investigating if higher-order
terms in the Taylor expansion, or leveraging different instantiations of the
PF ODE might lead to better performance. We hope that our work sparks greater
interest among practitioners in utilizing other foundational generative models
beyond DMs for OOD detection, such as large language models.

#### Limitations and Future Work.

We only calculate the simple time derivative and found that it works well
experimentally, although one might compute the full derivative to quantify the
curvature completely. We leave this for future work. Also, we consider a DM
trained on ImageNet, which may not be general enough for specialized
applications such as medical images. For such purposes, one could consider
including domain-specific data when training the DM, or utilize larger DMs
such as Stable Diffusion.

##  7 Acknowledgements

We would like to thank Pranav Goyal for help with the experiments. This
research is supported by the National Research Foundation Singapore and DSO
National Laboratories under the AI Singapore Programme (AISG Award No:
AISG2-RP-2020-017).

## References

  * Nguyen et al. [2015] Anh Nguyen, Jason Yosinski, and Jeff Clune.  Deep neural networks are easily fooled: High confidence predictions for unrecognizable images.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 427–436, 2015.
  * Nalisnick et al. [2018] Eric Nalisnick, Akihiro Matsukawa, Yee Whye Teh, Dilan Gorur, and Balaji Lakshminarayanan.  Do deep generative models know what they don’t know?  In _International Conference on Learning Representations_ , 2018.
  * Yang et al. [2021] Jingkang Yang, Kaiyang Zhou, Yixuan Li, and Ziwei Liu.  Generalized out-of-distribution detection: A survey.  _arXiv preprint arXiv:2110.11334_ , 2021.
  * Ren et al. [2019] Jie Ren, Peter J Liu, Emily Fertig, Jasper Snoek, Ryan Poplin, Mark Depristo, Joshua Dillon, and Balaji Lakshminarayanan.  Likelihood ratios for out-of-distribution detection.  _Advances in neural information processing systems_ , 32, 2019.
  * Choi et al. [2018] Hyunsun Choi, Eric Jang, and Alexander A Alemi.  Waic, but why? generative ensembles for robust anomaly detection.  _arXiv preprint arXiv:1810.01392_ , 2018.
  * Morningstar et al. [2021] Warren Morningstar, Cusuh Ham, Andrew Gallagher, Balaji Lakshminarayanan, Alex Alemi, and Joshua Dillon.  Density of states estimation for out of distribution detection.  In _International Conference on Artificial Intelligence and Statistics_ , pages 3232–3240. PMLR, 2021.
  * Graham et al. [2023] Mark S Graham, Walter HL Pinaya, Petru-Daniel Tudosiu, Parashkev Nachev, Sebastien Ourselin, and Jorge Cardoso.  Denoising diffusion models for out-of-distribution detection.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2947–2956, 2023.
  * Liu et al. [2023] Zhenzhen Liu, Jin Peng Zhou, Yufan Wang, and Kilian Q Weinberger.  Unsupervised out-of-distribution detection with diffusion inpainting.  In _International Conference on Machine Learning_ , pages 22528–22538. PMLR, 2023.
  * Choi et al. [2024] Sungik Choi, Hankook Lee, Honglak Lee, and Moontae Lee.  Projection regret: Reducing background bias for novelty detection via diffusion models.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Xiao et al. [2021] Zhisheng Xiao, Qing Yan, and Yali Amit.  Do we really need to learn representations from in-domain data for outlier detection?  _arXiv preprint arXiv:2105.09270_ , 2021.
  * Brown et al. [2020] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Bommasani et al. [2021] Rishi Bommasani, Drew A Hudson, Ehsan Adeli, Russ Altman, Simran Arora, Sydney von Arx, Michael S Bernstein, Jeannette Bohg, Antoine Bosselut, Emma Brunskill, et al.  On the opportunities and risks of foundation models.  _arXiv preprint arXiv:2108.07258_ , 2021.
  * Song et al. [2020a] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2020a.
  * Song and Ermon [2019] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  _Advances in neural information processing systems_ , 32, 2019.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _Advances in neural information processing systems_ , 33:6840–6851, 2020.
  * Serrà et al. [2019] Joan Serrà, David Álvarez, Vicenç Gómez, Olga Slizovskaia, José F Núñez, and Jordi Luque.  Input complexity and out-of-distribution detection with likelihood-based generative models.  In _International Conference on Learning Representations_ , 2019.
  * Hyvärinen and Dayan [2005] Aapo Hyvärinen and Peter Dayan.  Estimation of non-normalized statistical models by score matching.  _Journal of Machine Learning Research_ , 6(4), 2005.
  * Dockhorn et al. [2022] Tim Dockhorn, Arash Vahdat, and Karsten Kreis.  Genie: Higher-order denoising diffusion solvers.  _Advances in Neural Information Processing Systems_ , 35:30150–30166, 2022.
  * Song et al. [2020b] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  In _International Conference on Learning Representations_ , 2020b.
  * Kloeden and Platen [1992] Peter E. Kloeden and Eckhard Platen.  _Numerical Solution of Stochastic Differential Equations_.  Springer, Berlin, 1992.
  * Khrulkov et al. [2022] Valentin Khrulkov, Gleb Ryzhakov, Andrei Chertkov, and Ivan Oseledets.  Understanding ddpm latent codes through optimal transport.  In _The Eleventh International Conference on Learning Representations_ , 2022.
  * Su et al. [2022] Xuan Su, Jiaming Song, Chenlin Meng, and Stefano Ermon.  Dual diffusion implicit bridges for image-to-image translation.  In _The Eleventh International Conference on Learning Representations_ , 2022.
  * Särkkä and Solin [2019] Simo Särkkä and Arno Solin.  _Applied stochastic differential equations_ , volume 10.  Cambridge University Press, 2019.
  * Denouden et al. [2018] Taylor Denouden, Rick Salay, Krzysztof Czarnecki, Vahdat Abdelzad, Buu Phan, and Sachin Vernekar.  Improving reconstruction autoencoder out-of-distribution detection with mahalanobis distance.  _arXiv preprint arXiv:1812.02765_ , 2018.
  * Tack et al. [2020] Jihoon Tack, Sangwoo Mo, Jongheon Jeong, and Jinwoo Shin.  Csi: Novelty detection via contrastive learning on distributionally shifted instances.  _Advances in neural information processing systems_ , 33:11839–11852, 2020.
  * Sehwag et al. [2021] Vikash Sehwag, Mung Chiang, and Prateek Mittal.  Ssd: A unified framework for self-supervised outlier detection.  _arXiv preprint arXiv:2103.12051_ , 2021.
  * Nalisnick et al. [2019] Eric Nalisnick, Akihiro Matsukawa, Yee Whye Teh, and Balaji Lakshminarayanan.  Detecting out-of-distribution inputs to deep generative models using typicality.  _arXiv preprint arXiv:1906.02994_ , 2019.
  * Du and Mordatch [2019] Yilun Du and Igor Mordatch.  Implicit generation and modeling with energy based models.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Mahmood et al. [2020] Ahsan Mahmood, Junier Oliva, and Martin Andreas Styner.  Multiscale score matching for out-of-distribution detection.  In _International Conference on Learning Representations_ , 2020.
  * Zhou and Paffenroth [2017] Chong Zhou and Randy C Paffenroth.  Anomaly detection with robust deep autoencoders.  In _Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 665–674, 2017.
  * Schlegl et al. [2017] Thomas Schlegl, Philipp Seeböck, Sebastian M Waldstein, Ursula Schmidt-Erfurth, and Georg Langs.  Unsupervised anomaly detection with generative adversarial networks to guide marker discovery.  In _International conference on information processing in medical imaging_ , pages 146–157. Springer, 2017.
  * Zhang et al. [2018] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang.  The unreasonable effectiveness of deep features as a perceptual metric.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 586–595, 2018.
  * Song et al. [2023] Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever.  Consistency models.  In _International Conference on Machine Learning_ , pages 32211–32252. PMLR, 2023.
  * Nichol and Dhariwal [2021] Alexander Quinn Nichol and Prafulla Dhariwal.  Improved denoising diffusion probabilistic models.  In _International conference on machine learning_ , pages 8162–8171. PMLR, 2021.
  * Xiao et al. [2020] Zhisheng Xiao, Karsten Kreis, Jan Kautz, and Arash Vahdat.  Vaebm: A symbiosis between variational autoencoders and energy-based models.  In _International Conference on Learning Representations_ , 2020.
  * Du et al. [2021] Yilun Du, Shuang Li, Joshua Tenenbaum, and Igor Mordatch.  Improved contrastive divergence training of energy-based models.  In _International Conference on Machine Learning_ , pages 2837–2848. PMLR, 2021.
  * Kingma and Dhariwal [2018] Durk P Kingma and Prafulla Dhariwal.  Glow: Generative flow with invertible 1x1 convolutions.  _Advances in neural information processing systems_ , 31, 2018.
  * Song et al. [2021] Yang Song, Conor Durkan, Iain Murray, and Stefano Ermon.  Maximum likelihood training of score-based diffusion models.  _Advances in neural information processing systems_ , 34:1415–1428, 2021.
  * Peyré and Cuturi [2019] Gabriel Peyré and Marco Cuturi.  Computational optimal transport.  _Foundations and Trends in Machine Learning_ , 11(5-6):355–607, 2019.

Supplementary Material for “Out-of-Distribution Detection with a Single
Unconditional Diffusion Model"

##  Appendix A Proofs

###  A.1 Theorem 1

See 1

###### Proof.

The proof is a modification from Song et al. [38]. Let us first state the PF
ODEs of the two distributions explicitly:

$\displaystyle\frac{{\mathrm{d}}{\mathbf{x}}_{t}}{{\mathrm{d}}t}$
$\displaystyle={\bm{f}}({\mathbf{x}}_{t},t)+\frac{g(t)^{2}}{2\sigma_{t}}\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t),\
\ \
\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t)=-\sigma_{t}\nabla_{\mathbf{x}}\log\phi_{t}({\mathbf{x}})$
(11) $\displaystyle\frac{{\mathrm{d}}{\mathbf{x}}_{t}}{{\mathrm{d}}t}$
$\displaystyle={\bm{f}}({\mathbf{x}}_{t},t)+\frac{g(t)^{2}}{2\sigma_{t}}\bm{\epsilon}_{\psi}({\mathbf{x}}_{t},t),\
\ \
\bm{\epsilon}_{\psi}({\mathbf{x}}_{t},t)=-\sigma_{t}\nabla_{\mathbf{x}}\log\psi_{t}({\mathbf{x}}).$
(12)

We make the following assumption about $\phi_{t}$ and $\psi_{t}$:

$\forall t\in[0,T]\ \exists k>0\textrm{ s.t.
}\phi_{t}({\mathbf{x}})=O(e^{-\left\lVert{\mathbf{x}}\right\rVert^{k}_{2}}),\
\
\psi_{t}({\mathbf{x}})=O(e^{-\left\lVert{\mathbf{x}}\right\rVert^{k}_{2}})\textrm{
as }\left\lVert{\mathbf{x}}\right\rVert_{2}\rightarrow\infty.$ (13)

This assumption essentially states that both marginals have at most polynomial
growth at infinity and thus can be normalized.

We start by rewriting the KL divergence between $\phi_{0}$ and $\psi_{0}$ in
integral form:

$\displaystyle\begin{split}D_{\mathrm{KL}}(\phi_{0}\|\psi_{0})&=D_{\mathrm{KL}}(\phi_{0}\|\psi_{0})-D_{\mathrm{KL}}(\phi_{T}\|\psi_{T})+D_{\mathrm{KL}}(\phi_{T}\|\psi_{T})\\\
&=-\int_{0}^{T}\frac{\partial D_{\mathrm{KL}}(\phi_{t}\|\psi_{t})}{\partial
t}{\mathrm{d}}t+D_{\mathrm{KL}}(\phi_{T}\|\psi_{T}).\end{split}$ (14)

As we can treat the PF ODE as a special case of an SDE with zero diffusion
term, we can obtain the Fokker-Planck of the marginal density of the PF ODEs,
also known as the continuity equation, as follows:

$\frac{\partial\phi_{t}}{\partial
t}=\nabla_{\mathbf{x}}\cdot\left(-{\bm{f}}({\mathbf{x}}_{t},t)\phi_{t}({\mathbf{x}})-\frac{g(t)^{2}}{2\sigma_{t}}\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t)\phi_{t}({\mathbf{x}})\right)=\nabla_{\mathbf{x}}\cdot\left({\bm{h}}_{\phi}\phi_{t}({\mathbf{x}})\right)$
(15)

where we define
${\bm{h}}_{\phi}\coloneq-{\bm{f}}({\mathbf{x}}_{t},t)-\frac{g(t)^{2}}{2\sigma_{t}}\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t)$
for simplicity. Similarly, $\frac{\partial\psi_{t}}{\partial
t}=\nabla_{\mathbf{x}}\cdot({\bm{h}}_{\psi}\psi_{t}({\mathbf{x}}))$. Let us
now rewrite the time-derivative $\frac{\partial
D_{\mathrm{KL}}(\phi_{t}\|\psi_{t})}{\partial t}$ in Eq. 14 as follows:

$\displaystyle\frac{\partial D_{\mathrm{KL}}(\phi_{t}\|\psi_{t})}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
t}\int\phi_{t}({\mathbf{x}})\log\frac{\phi_{t}({\mathbf{x}})}{\psi_{t}({\mathbf{x}})}{\mathrm{d}}{\mathbf{x}}$
$\displaystyle=\int\frac{\partial\phi_{t}}{\partial
t}\log\frac{\phi_{t}({\mathbf{x}})}{\psi_{t}({\mathbf{x}})}{\mathrm{d}}{\mathbf{x}}+\underbrace{\int\frac{\partial\phi_{t}({\mathbf{x}})}{\partial
t}{\mathrm{d}}{\mathbf{x}}}_{=0}-\int\frac{\phi_{t}({\mathbf{x}})}{\psi_{t}({\mathbf{x}})}\frac{\partial\psi_{t}({\mathbf{x}})}{\partial
t}{\mathrm{d}}{\mathbf{x}}$
$\displaystyle=\int\nabla_{\mathbf{x}}\cdot\left({\bm{h}}_{\phi}({\mathbf{x}},t)\phi_{t}({\mathbf{x}})\right)\log\frac{\phi_{t}({\mathbf{x}})}{\psi_{t}({\mathbf{x}})}{\mathrm{d}}{\mathbf{x}}-\int\frac{\phi_{t}({\mathbf{x}})}{\psi_{t}({\mathbf{x}})}\nabla_{\mathbf{x}}\cdot({\bm{h}}_{\psi}({\mathbf{x}},t)\psi_{t}({\mathbf{x}})){\mathrm{d}}{\mathbf{x}}$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\int\phi_{t}({\mathbf{x}})[{\bm{h}}_{\phi}^{\top}({\mathbf{x}},t)-{\bm{h}}_{\psi}^{\top}({\mathbf{x}},t)][\nabla_{\mathbf{x}}\log\phi_{t}({\mathbf{x}})-\nabla_{x}\log\psi_{t}({\mathbf{x}})]{\mathrm{d}}{\mathbf{x}}$
$\displaystyle=-\frac{1}{2}\int\phi_{t}({\mathbf{x}})\frac{g(t)^{2}}{\sigma_{t}}\left\lVert\bm{\epsilon}_{\phi}({\mathbf{x}},t)-\bm{\epsilon}_{\psi}({\mathbf{x}},t)\right\rVert_{2}^{2}{\mathrm{d}}{\mathbf{x}},$

where (i) is due to integration by parts and the fact that
$\lim_{{\mathbf{x}}\rightarrow\infty}{\bm{h}}_{\phi}({\mathbf{x}},t)\phi_{t}({\mathbf{x}})=0$
and
$\lim_{{\mathbf{x}}\rightarrow\infty}{\bm{h}}_{\psi}({\mathbf{x}},t)\psi_{t}({\mathbf{x}})=0$
due to assumption 13. Combining with Eq. 14 gives us the desired result:

$D_{\mathrm{KL}}(\phi_{0}\|\psi_{0})=\frac{1}{2}\int_{0}^{T}\mathbb{E}_{x\sim\phi_{t}}\frac{g(t)^{2}}{\sigma_{t}}\left\lVert\bm{\epsilon}_{\phi}({\mathbf{x}}_{t},t)-\bm{\epsilon}_{\psi}({\mathbf{x}}_{t},t)\right\rVert^{2}_{2}{\mathrm{d}}t+D_{\mathrm{KL}}(\phi_{T}\|\psi_{T}).$
(16)

∎

###  A.2 OT Toy Example

We derive here in detail the toy example discussed in Sec. 3.4, where we will
prove the optimal transport map between a multivariate Gaussian and standard
normal is identical to the diffusion PF ODE path. We consider our source
distribution as
$p_{0}({\mathbf{x}})\sim\mathcal{N}({\mathbf{a}},{\mathbf{I}})$,
${\mathbf{a}}\in\mathbb{R}^{d}\textrm{ and }{\mathbf{I}}\in\mathbb{R}^{d\times
d}$. We choose our forward SDE to be parameterized as:

${\mathrm{d}}{\mathbf{x}}_{t}=-{\mathbf{x}}{\mathrm{d}}t+\sqrt{2}{\mathrm{d}}{\mathbf{w}}_{t},$
(17)

which is the same Ornstein–Uhlenbeck process as the DDPM forward SDE Eq. 4
with a constant noise schedule $\beta(t)=2$. This is also commonly referred to
as the Langevin equation.

As Eq. 17 has affine drift coefficients and a starting distribution which is
normal, we know that the marginal distributions at intermediate times are also
normal, $p_{t}({\mathbf{x}})\sim\mathcal{N}(\bm{\mu}(t),\bm{\Sigma}(t))$.
Furthermore, we can calculate the means and variances analytically using Eq.
5.50 and Eq. 5.51 of Särkkä and Solin [23]:

$\frac{{\mathrm{d}}\bm{\mu}(t)}{{\mathrm{d}}t}=-\bm{\mu}(t),\ \ \ \ \
\frac{{\mathrm{d}}\bm{\Sigma}(t)}{{\mathrm{d}}t}=-2\bm{\Sigma}(t)+2$ (18)

with solutions

$\bm{\mu}(t)=\bm{\mu}(0)e^{-t}={\mathbf{a}}e^{-t},\ \ \ \
\bm{\Sigma}(t)={\mathbf{I}}+e^{-2t}(\bm{\Sigma}(0)-{\mathbf{I}})={\mathbf{I}}.$
(19)

Thus, the marginal density has the form
$p_{t}({\mathbf{x}})=\mathcal{N}({\mathbf{a}}e^{-t},{\mathbf{I}})$, from which
we compute the score as $\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}})=-{\mathbf{x}}+{\mathbf{a}}e^{-t}$. We can substitute this
into the corresponding PF ODE to obtain:

$\displaystyle\begin{split}\frac{{\mathrm{d}}{\mathbf{x}}_{t}}{{\mathrm{d}}t}&=-{\mathbf{x}}-\nabla_{\mathbf{x}}\log
p_{t}({\mathbf{x}})\\\ &=-{\mathbf{a}}e^{-t}.\end{split}$ (20)

The optimal transport map denoted $E_{p_{0}}({\mathbf{x}})$ is obtained by
solving Eq. 20 to get
${\mathbf{x}}_{t}={\mathbf{x}}_{0}+{\mathbf{a}}(e^{-t}-{\mathbf{I}})$ and
taking the limit $t\rightarrow\infty$. This gives us
$E_{p_{0}}({\mathbf{x}})={\mathbf{x}}-{\mathbf{a}}$, which is precisely the OT
map between $p_{0}({\mathbf{x}})\sim\mathcal{N}({\mathbf{a}},{\mathbf{I}})$
and $\mathcal{N}(\bm{0},{\mathbf{I}})$ (cf. Eq. 2.40 in Peyré and Cuturi
[39]).

##  Appendix B Experimental Details

#### DiffPath

As mentioned in Sec. 5, we utilize a single unconditional diffusion model
trained on ImageNet at $64\times 64$ resolution. The checkpoint that we use is
trained using Improved DDPM’s $L_{\textrm{hybrid}}$ objective (Eq. 16 of
Nichol and Dhariwal [34]). The model uses a cosine noise schedule with a total
of 4000 diffusion steps. The checkpoint is publicly available on the official
GitHub repository333https://github.com/openai/improved-diffusion. For
DiffPath-1D, we fit a KDE using a Gaussian kernel with a bandwith of 5. For
DiffPath-6D, we fit a GMM using 50 mixture components, each with diagonal
covariance. Both are implemented using the sklearn library. We calculate the
diffusion statistics using a single A5000 GPU.

#### Diffusion Baselines

For all diffusion baselines, we rely on the official GitHub repositories and
open-source checkpoints where possible. The repositories are listed as
follows: MSMA444https://github.com/ahsanMah/msma, DDPM-
OOD555https://github.com/marksgraham/ddpm-ood,
LMD666https://github.com/zhenzhel/lift_map_detect. For NLL and IC, we either
use the pre-trained CIFAR10 checkpoint from Improved DDPM or train separate
models for SVHN and CelebA using the same hyperparameters at $32\times 32$
resolution. We calculate the NLL using the default implementation in Improved
DDPM, while we compute IC using code from the LR
repository777https://github.com/XavierXiao/Likelihood-Regret due to lack of
official code from the IC authors. We train all baselines using 1-3 A5000
GPUs.

###  B.1 Full Results for Ablations

We provide full results of the ablations presented in Fig. 4 in Table 5 and
Table 6. When varying the NFE in Table 5, we observe that certain tasks are
more sensitive than others, although the average across all tasks is much more
stable. For instance, the worst performance between CIFAR10 vs SVHN has AUROC
of 0.834 at 500 steps, while the best performance has AUROC of 0.936 at 50
steps. However, other tasks remain unchanged when NFE is varied, such as SVHN
vs Textures or CelebA vs SVHN. We generally find 50 and 100 steps to the sweet
spot.

When the resizing method is varied in Table 6, we find bicubic and bilinear
interpolation to perform significantly better than the rest, with bilinear
achieving the best scores in most tasks. As such, we use bilinear in all other
experiments in the main text.

Table 5: Full ablation results when we vary the number of DDIM steps (NFE). We
use DiffPath-6D.  Bold and underline denotes the best and second best result.

| C10 vs | SVHN vs | CelebA vs |
---|---|---|---|---
NFE | SVHN | CelebA | C100 | Textures | C10 | CelebA | C100 | Textures | C10 | SVHN | C100 | Textures | Average
25 | 0.917 | 0.996 | 0.587 | 0.997 | 0.919 |  0.999 | 0.929 | 1 | 0.995 | 1 | 0.993 | 0.970 |  0.942
50 | 0.936 | 0.997 | 0.596 | 0.997 |  0.921 |  0.999 |  0.932 |  0.999 |  0.997 | 1 | 0.995 | 0.935 |  0.942
100 |  0.920 |  0.999 |  0.593 | 0.999 | 0.924 |  0.999 | 0.936 | 1 | 0.999 | 1 |  0.998 | 0.946 | 0.943
200 | 0.882 | 1 | 0.588 | 0.999 | 0.912 | 1 | 0.926 | 1 | 0.999 | 1 | 0.999 |  0.973 | 0.940
500 | 0.834 |  0.999 | 0.583 | 0.999 | 0.850 | 1 | 0.870 | 1 | 0.999 | 1 |  0.998 | 0.981 | 0.926

Table 6: Full ablation results when we vary the resizing method of the images.
We use DiffPath-6D.  Bold and underline denotes the best and second best
result.

| C10 vs | SVHN vs | CelebA vs |
---|---|---|---|---
Resizing Method | SVHN | CelebA | C100 | Textures | C10 | CelebA | C100 | Textures | C10 | SVHN | C100 | Textures | Average
Lanczos | 0.819 | 0.997 |  0.593 |  0.994 | 0.891 |  0.996 | 0.910 |  0.996 | 0.997 | 1 | 0.996 | 0.937 | 0.927
Nearest | 0.815 |  0.995 | 0.575 | 0.973 | 0.836 | 0.999 | 0.865 | 0.993 | 0.979 |  0.995 | 0.983 | 0.983 | 0.916
Bicubic |  0.900 | 0.988 | 0.589 | 0.987 |  0.899 | 0.994 |  0.916 |  0.996 |  0.986 | 1 | 0.979 | 0.930 |  0.930
Hamming | 0.823 | 0.898 | 0.575 | 0.937 | 0.852 | 0.990 | 0.877 | 0.985 | 0.835 | 0.986 | 0.862 |  0.938 | 0.880
Box | 0.814 | 0.974 | 0.574 | 0.957 | 0.835 | 0.995 | 0.866 | 0.991 | 0.951 | 0.987 | 0.962 |  0.938 | 0.904
Bilinear | 0.936 | 0.997 | 0.596 | 0.997 | 0.921 | 0.999 | 0.932 | 0.999 | 0.997 | 1 |  0.995 | 0.935 | 0.942


# Using Surprise Index for Competency Assessment in Autonomous Decision-Making

Akash Ratheesh 111Professional Research Assistant, Ann and H.J. Smead
Aerospace Engineering Sciences, University of Colorado Boulder University of
Colorado Boulder, 429 UCB, Boulder, CO, 80309  Ofer Dagan 222Graduate Research
Assistant, Ann and H.J. Smead Aerospace Engineering Sciences, University of
Colorado Boulder University of Colorado Boulder, 429 UCB, Boulder, CO, 80309
Nisar R. Ahmed 333Associate Professor, Ann and H.J. Smead Aerospace
Engineering Sciences, University of Colorado Boulder Natasha Bosanac
444Assistant Professor, Ann and H.J. Smead Aerospace Engineering Sciences,
University of Colorado Boulder University of Colorado Boulder, 429 UCB,
Boulder, CO, 80309  Jay McMahon‡ University of Colorado Boulder, 429 UCB,
Boulder, CO, 80309

###### Abstract

This paper considers the problem of evaluating an autonomous system’s
competency in performing a task, particularly when working in dynamic and
uncertain environments. The inherent opacity of machine learning models, from
the perspective of the user, often described as a ‘black box’, poses a
challenge. To overcome this, we propose using a measure called the Surprise
index, which leverages available measurement data to quantify whether the
dynamic system performs as expected. We show that the surprise index can be
computed in closed form for dynamic systems when observed evidence in a
probabilistic model if the joint distribution for that evidence follows a
multivariate Gaussian marginal distribution. We then apply it to a nonlinear
spacecraft maneuver problem, where actions are chosen by a reinforcement
learning agent and show it can indicate how well the trajectory follows the
required orbit.

##  1 Introduction

Over the past few years, there has been increasing use of autonomous systems
for various applications, from autonomous driving to aerospace systems. These
systems are often required to operate in complex and uncertain environments,
e.g., the highway or in space. For successful deployment and useful operation
of these systems around and by humans, there is a need to evaluate the
autonomous system’s competency in performing its task. For this reason [1, 2,
3] present a _Factorized Machine Self-Confidence (FaMSeC)_ framework, aimed at
allowing a decision-making autonomous system to self-asses its competency.
This paper focuses on one of the factors considered by that framework, the
model validity, which is meant to assess how well the underlying model, used
to make decisions, reflects reality. One key challenge is that the model,
often generated by some machine learning (ML) algorithm, might be a “black
box” with respect to the end-user. An intuitive solution is to use available
data in the form of measurements and assess whether what we see is more
surprising than what we expect to see. The surprise index (SI) [4] is a
measure that aims to quantify this intuition. For a set of measurements, it
sums up the probability of all other _possible_ measurement instances that
have a probability less than the observed ones. A very probable (less
surprising) set of measurements will then have more probability mass summed
‘behind’ it.

An important point to note here is that this automatically gives a notion of
an equivalence class for the instances of observed measurements based on their
probability. That is, there are ‘typical sets’ of observations or evidence
that will not appear surprising to the model, and there will be other
‘atypical’ sets of observations or evidence that will appear to be more
surprising. The (a)typicality of a set of evidence values can be denoted by
the probability of seeing available evidence. Since it is generally difficult
to tell apart outliers and extreme events (e.g., for long-tailed
distributions) from anomalies that are not explained well by the model, the
surprise index allows us to see if a model does not do a good job of
explaining the evidence.

This has much in common with NEES/NIS Monte Carlo testing for Kalman filter
consistency checking [5]. That is, we are interested in asking the same
question: does the estimated uncertainty of the filter/probabilistic model
correctly describe the actual uncertainty encountered in the ‘real world’
implementation for the problem the model is constructed against? However, we
will show later in this paper that there are subtle differences in how we go
about answering that question with the surprise index vs. NEES/NIS tests.

As noted in [6, 7], the SI is not easy to compute, and it is generally
intractable since it is exponential in the number of evidence variables. For
finitely countable discrete random variable distributions, the SI could be (in
principle) evaluated by assigning a total ordering to the outcomes of the
random variable, according to the probability of outcome occurrence. Then the
corresponding cumulative density function (cdf) gives the desired SI. In Bayes
Nets and other probabilistic models, this may not be so easy to do since the
number of observations and possible outcomes of each observation could be
prohibitively large. Furthermore, sorting by the probabilities requires
knowing the probabilities of the joint evidence outcomes in the first place,
which is typically not available in closed form – indeed, this is what makes
Bayesian inference approximations necessary in the first place for such
models. To date and to the best of our knowledge, there have not been any
results describing how the SI can be exactly computed, even for special cases.
The main contributions of this paper are:

  1. 1.

We develop a closed-form solution of SI for a multivariate Gaussian
observation model and dynamical system.

  2. 2.

We test our closed-form solution on a $2D$ GPS localization example to
validate it against the empirical SI [6].

  3. 3.

We demonstrate the application of SI to assess the validity of a new
reinforcement learning-based nonlinear spacecraft maneuver planner [8, 9].

The rest of the paper is organized as follows: in Sec. 2, we motivate the use
of SI in the context of “black box" control systems and give the necessary
background on SI. Sec. 3 develops a closed-form solution of SI, validates it
against an empirical solution, and extends to linear dynamic systems. We then
apply the solution for a nonlinear spacecraft maneuver scenario in Sec. 4, and
draw conclusions regarding the performance of a reinforcement learning-based
control policy. Sec. 5 then presents conclusions and details future research
directions.

##  2 Background

Reinforcement learning (RL) based methods are being increasingly explored in
the areas of spacecraft control and trajectory planning. In real-world
spacecraft operations, the true states of the spacecraft are unknown to the
user/‘trajectory planner,’ and maneuvering operations are planned based on an
estimate of the spacecraft states. Previous works [8, 9] have used RL-based
techniques to perform maneuver planning for periodic orbits in the Earth-Moon
System. These orbits are unstable and are very sensitive to observation and
navigation errors. Thus, any errors in the maneuvring action generated by
these RL-based planners can have a significant impact on the resulting
trajectory. Also, the RL-based planners are trained in a model-free manner,
and true dynamical models are not available to the planner. Therefore it is
crucial to evaluate these RL agents to validate that the actions generated are
‘what is expected’ from the system. Secondly, although a reward function can
sufficiently indicate whether a policy is ‘well-trained,’ it may not capture
any off-nominal behaviors that can happen during the execution of the
trajectory, e.g., deviations due to thruster misfire/malfunction.

The surprise index can be useful in these scenarios to evaluate these
trajectories in an online manner to assess whether the observations of the
trajectory resulting from the action match the expected behavior from the
system in a probabilistic sense. Suppose we have some probabilistic model, for
instance a Bayes Net (BN) or Markov random field (MRF) or something else which
stipulates a joint probability density function (pdf) of the variables of
interest, with observed evidence variables
$E_{1:N}=\\{E_{1},E_{2},...,E_{N}\\}$. Let the specific observed instantiation
of those evidence variables be denoted
$e^{o}_{1:N}=\\{e^{0}_{1},e^{0}_{2},...,e^{0}_{N}\\}$, where
$e^{o}_{i}\in\text{dom}(E_{i})$ and dom is the domain of $E_{i}$. We will use
the superscript _o_ notation to denote the observed instance and remove the
_o_ to talk about other hypothetical instances
$e_{1:N}=\\{e_{1},e_{2},...,e_{N}\\}$ which ‘could’ have been observed in
principle (this is a key distinction). Note that it can be the case that the
evidence variables are either discrete, continuous or some mix of both.

The surprise index of an observed instance of evidence variables $e^{o}_{1:N}$
is defined as the sum of the probability of all other events that are less
likely to occur,

$\mathcal{SI}(e^{0}_{1:N})=\sum_{e_{1:N}:P(e_{1:N})<P(e^{0}_{1:N})}P(e_{1:N}).$
(1)

Note that since it sums the probability of all other possible instances
$e_{1:N}$ of $E_{1:N}$ which have a probability less than the observed
instance $e^{o}_{1:n}$, it automatically defines the range of the surprise
index to be on the interval [0,1], with a higher value indicating less
surprise. In other words, the more probability mass that an observed instance
has sitting ‘behind’ it, the less surprising it appears to be, according to
the probabilistic model.

The problem is that, in practice, computing the surprise index requires
iterating through all combinations of $e_{1:N}$, which is exponential in the
number of evidence variables. Thus it is intractable, as noted in [6] and in
[7]. In [6], Zagorecki et al. proposed an asymptotic approximation to the SI,
which makes use of the fact that the joint distribution of the evidence
follows an approximately log-normal distribution (i.e., the logarithm of the
joint evidence probabilities has a Gaussian distribution, owing to the central
limit theorem). Zagorecki et al. claim the approximation is particularly
useful when very large $N$ is encountered, i.e., for large networks. In [7],
Laskey mentions the SI in the context of original work by Habbema [4], but
essentially goes no further than saying that it is intractable to compute,
proceeding instead to look at different heuristics such as Jensen et al.’s
conflict score [10], which is closely related to the mutual information of the
evidence variables.

Nevertheless, as discussed in Zagorecki et al.[6], the surprise index has
several useful properties that make it a useful measure for model quality
assessment in machine competency self-assessment. These properties also make
it preferable in many respects to alternative goodness of fit type measures
for probabilistic models, such as log-likelihood scores (which can become
unboundedly large for a large number of observations) or entropy (for which
minimum and maximum values depend on the number of outcomes described by a
probability distribution). Unlike these measures, the surprise index retains
an interpretable meaning as a probability and so always takes on the range of
0 to 1 regardless of the number of outcomes. Moreover, the surprise index has
a very close connection to the concept of p-values from statistics, which are
also widely used for goodness of fit testing. Unlike p-values, however, the
surprise index does not require hypothesis tests to attach a specific meaning
to its calculation.

In the next section, we show that the surprise index can actually be computed
in closed form for observed evidence in a probabilistic model if the joint
distribution for that evidence follows a multivariate Gaussian marginal
distribution. The significance of this is to connect the surprise index to the
well-known chi-square hypothesis tests for things like Kalman filters (NIS
tests for filter consistency, specifically) [5]. As we’ll explore later, the
surprise index may provide a useful complementary way of assessing ‘goodness
of fit’ (g.o.f.) for a probabilistic model that does not rely on the ability
of that model to actually process data (i.e., do inference, which the NIS test
for Kalman filters does require). Hence, it could be a useful alternative
metric for tuning things like Kalman filters, extended Kalman filters, and
other probabilistic inference tools for estimation and control (e.g., LQG,
RL).

##  3 Approach and Methodology

Consider a probabilistic model where the evidence/observation distribution is
a continuous N-variate Gaussian pdf. Suppose we observe random vectors
$E_{i}=Y_{i}\in\mathbb{R}^{p\times 1}$, so that $E_{1:N}=Y_{1:N}$ can be
represented as a $p\times N$ array of vectors, where

$Y_{1:N}\sim\mathcal{N}_{y_{1:N}}(\mu,\Sigma).$ (2)

We can consider the realization, $y_{1:N}$ to be a concatenated vector of
$p\times 1$ random vector stacked in to $N$ blocks, where
$\mu\in\mathbb{R}^{Np\times 1}$ and $\Sigma\in\mathbb{R}^{Np\times Np}$. This
describes the joint distribution for the possible sets of $y_{1:N}$ instances
of $Y_{1:N}$. Similarly, the pdf for the realized (observed) instances
$y^{o}_{1:N}$ would be

$P(y^{o}_{1:N})=\mathcal{N}_{y^{o}_{1:N}}(\mu,\Sigma).$ (3)

In order to calculate SI according to Eq. (1), we need to order all possible
sets of observations $E_{1:N}\in S$ and sum/integrate the ones that have less
probability than the observed set $E^{o}_{1:N}\in S^{o}$. For example, if we
consider two non-overlapping regions for integration $S_{a}$ and $S_{b}$,
where the realizations for $Y_{1:N}$ in $S_{a}$ are equally more probable
(less surprising) than the realizations in $S_{b}$ (more surprising), then

$p(y_{1:N}\in
S_{a})=\int_{S_{a}}\mathcal{N}_{y_{1:N}}(\mu,\Sigma)dy_{1:N}>p(y_{1:N}\in
S_{b})=\int_{S_{b}}\mathcal{N}_{y_{1:N}}(\mu,\Sigma)dy_{1:N}.$ (4)

Fortunately, for multivariate Gaussian pdf, such ordered regions do, in fact,
exist, and their boundaries are given by the level sets defined by the
Mahalanobis distance. The Mahalanobis distance $d_{M}$ for the joint
realization vector $y_{1:N}$ is given by

$d_{M}=\frac{1}{2}(y_{1:N}-\mu)^{T}\Sigma^{-1}(y_{1:N}-\mu).$ (5)

To demonstrate, suppose that $d_{1}>0$ and $d_{2}=d_{1}+\delta$ for some small
$\delta\in[0,d_{1}]$, we can define $S_{a}$ and $S_{b}$ as follows

$\begin{split}&S_{a}=\\{y_{1:N}:d_{M}\leq d_{1}\\}\\\
&S_{b}=\\{y_{1:N}:d_{1}\leq d_{M}\leq d_{2}\\},\end{split}$ (6)

then for $N=2$ computing the probabilities according to Eq. (4) leads directly
to the dark and lighter purple regions shown in Fig. 1, corresponding to
$S_{a}$ and $S_{b}$, respectively. Note that this readily extends to higher
dimensions, where $S_{a}$ and $S_{b}$ are defined as the contours of hyper-
ellipses in $\mathbb{R}^{N_{p}\times 1}$ about the mean $\mu$, with
orientation defined by $\Sigma$.

Figure 1: Ordered probabilistic equivalence sets for Gaussian pdf with $N=2$.

The question now is how to compute the probabilities and resulting SI for a
specific set of observations $y_{1:N}^{o}$ that have a joint Gaussian pdf. We
know that the quadratic form for Gaussian random variables given by the
Mahalanobis distance results in a chi-square distributed scalar random
variable with $Np$ degree of freedom [5],

$\epsilon=(y_{1:N}-\mu)^{T}\Sigma^{-1}(y_{1:N}-\mu)\sim\chi^{2}_{\epsilon}(Np).$
(7)

Notice that the chi-square variable $\epsilon$ and the Mahalanobis distance
$d_{M}$ (Eq. (5)) are related through $d_{M}=\frac{1}{2}\epsilon$, or
$\epsilon=2d_{M}$. Thus, all possible observation sets, $S_{r}$, that are less
probable than the observed set $y_{1:N}^{o}$, will have larger Mahalanobis
distance, i.e. larger $\epsilon$. But since $S_{r}$ is unbounded, we can
instead look at the complimentry555This evokes the commonly used ‘observation
gating probability’ found in the data association literature, e.g., for
probabilistic data association methods in multi-target tracking [11], more
probable, set $\bar{S}_{O}=\\{y_{1:N}:d_{M}(y_{1:N})\leq
d_{M}(y^{o}_{1:N})\\}=\mathbb{R}^{Np\times 1}-\bar{S_{r}}$. In Fig. 1, if
$d_{M}(y^{o}_{1:N})=d_{O}=d_{1}$, this means that instead of considering all
level sets that lie outside of $d_{1}$, i.e., for different $S_{b}$, we
consider the ones inside $S_{a}=\bar{S}_{O}$. It follows that

$\begin{split}\mathcal{SI}(y^{o}_{1:N})&=P(y_{1:N}\in
S_{r})=\int_{S_{r}}\mathcal{N}_{y_{1:N}}(\mu,\Sigma)dy_{1}...dy_{N}=1-P(y_{1:N}\in
S_{O})\\\
&=1-\int_{S_{O}}\mathcal{N}_{y_{1:N}}(\mu,\Sigma)dy_{1}...dy_{N}=1-P(\epsilon<\epsilon_{O})=1-\int^{\epsilon_{O}}_{0}\chi^{2}_{\epsilon}(Np)d\epsilon.\end{split}$
(8)

Note that this is just the complementary cdf of the chi-square distribution.
Also note that the maximum value of $\mathcal{SI}=1$ (least surprising) is
achieved at $y^{o}=\mu$, whereas the minimum value of $\mathcal{SI}=0$ (most
surprising) is achieved asymptotically as the elements of $y^{o}$ approach
$\pm\infty$.

At this point, it is worth re-stating the interpretation of the surprise
index: it denotes the probability of obtaining any realization for $Y_{1:N}$
that is less probable than the observed realization sequence $y^{o}_{1:N}$
under the assumed model. Hence, for a jointly Gaussian distributed set of
observations, the probability of obtaining any other possible observation
realization that is less probable than a given observation decreases outwardly
from the mean, in agreement with intuition.

###  3.1 Validation Example - GPS Localization

Suppose we have a 2-D GPS-like localization sensor that provides noisy
measurements $y_{k}\in\mathbb{R}^{2}$ for a stationary platform located at
position $x\in\mathbb{R}^{2}$. Consider that the measurements generated by the
sensor are defined by the following prior probabilistic model

$y_{k}=x+v_{k},\quad x\sim\mathcal{N}(\mu_{0},\Sigma_{0}),\quad
v_{k}\sim\mathcal{N}(0,R).$ (9)

Here, $\mu_{0}$ is the prior mean of the platform’s position, $\Sigma_{0}$ is
the prior covariance in the platform position, and $R$ is the true discrete-
time additive white Gaussian noise (AWGN) covariance of the sensor.

Now suppose we record a fixed-length realization of $Y_{1:N}$,
$y^{o}_{1:N}=\\{y^{o}_{1},...,y^{o}_{N}\\}$. In order to calculate the
Surprise Index for the successive sequence $y^{o}_{1:k}$ for $1\leq k\leq N$,
we have to obtain the prior probability of observing sequences $y_{1:k}$ that
are less probable than $y^{o}_{1:k}$. Figure 2 shows the comparison between
the SI calculated using a grid-based brute force approximation – not
dissimilar to the empirical method described in [6] – and the analytical
method we developed (8). We can see good agreement between the two methods,
thus giving initial validation to our new analytical solution.

Note that by the definition of the SI, we do not need to filter or process the
data at all to produce a state estimate for $x$. Rather, we ‘just’ need to
assess the a priori probability of obtaining sequences $y_{1:k}$ that are less
probable than $y^{o}_{1:k}$. This is important to consider in the context of
assessing models for processing data with techniques like batch least-squares
or Kalman filtering since typical ‘residual’ statistics for these methods are
based on measures of a posteriori fit or surprise (NEES/NIS).

Figure 2: Empirical grid approximation of SI according to [6] (left) and the
new analytical SI (right).

###  3.2 Generalization to Dynamic Filtering Problems

We derive the general expressions for computing the SI with a discrete-time
linear Gaussian state space system, i.e., a linear Gauss-Markov process.
Consider a discrete-time linear Gaussian state space system with known single-
step state transition matrix $F_{k}$, control parameter $G_{k}$, process
injection parameter $\Gamma_{k}$ and observation matrix $H_{k}$,

$\begin{split}&x_{k+1}=F_{k}x_{k}+G_{k}u_{k}+\Gamma_{k}w_{k},\quad
w_{k}\sim\mathcal{N}(0,Q_{k})\\\ &y_{k+1}=H_{k+1}x_{k+1}+v_{k+1},\quad
v_{k+1}\sim\mathcal{N}(0,R_{k+1}),\end{split}$ (10)

where $x$ is the state, $y$ is the sensor measurements, $u_{k}$ are the known
control inputs at time step $k$, $w_{k}$ and $v_{k+1}$ are AWGN process noise
and measurement noise with covariance matrix $Q_{k}$ and $R_{k+1}$
respectively. Suppose the initial states of the system are unknown and modeled
as Gaussian prior $x_{0}=\mathcal{N}(\mu_{x_{0}},\Sigma_{x_{0}})$ where
$\mu_{x_{0}}$ is the mean and $\Sigma_{x_{0}}$ is the apriori covariance.

We need to compute the marginal joint pdf of the observation sequence vector
$\bar{y}_{1:N}\equiv[\bar{y}^{T}_{1},...,\bar{y}^{T}_{N}]$,

$p(\bar{y}_{1},...,\bar{y}_{N})=\mathcal{N}_{y_{1:N}}\left(\mu_{y_{1:N}},\Sigma_{y_{1:N}}\right).$
(11)

If we assume that the system does not have any process noise, we can compute
the expected measurement $\mu_{y_{k}}\in\mu_{y_{1:N}}$ for $1\leq k\leq N$ by,

$\mu_{y_{k}}=H_{k}x_{k}=H_{k}[\Phi(0,k)\mu_{x_{0}}]$ (12)

where $\Phi$ is the state transition function,

$\Phi(k-T,k)=\prod^{T}_{i=1}F_{k-i},\quad\Phi(k,k)=I$ (13)

In order to compute the marginal joint pdf of the observation sequence, we
also need to find the covariance matrix $\Sigma_{y_{1:N}}$. The block diagonal
and block off-diagonal term of $\Sigma_{y_{1:N}}$ is given by

$\begin{split}\mathbb{E}\left[(y_{i}-\mu_{i})(y_{i}-\mu_{i})^{T}\right]=\mathbb{E}[y_{i}y_{i}^{T}]-\mu_{yi}\mu^{T}_{yi}=R_{i}+H_{i}\left[\Phi(0,i)\Sigma_{x_{0}}\Phi^{T}(0,i)\right]H^{T}_{i},\end{split}$
(14)
$\begin{split}\mathbb{E}\left[(y_{i}-\mu_{i})(y_{j}-\mu_{j})^{T}\right]=\mathbb{E}[y_{i}y_{j}^{T}]-\mu_{yi}\mu^{T}_{yj}=H_{i}\left[\Phi(0,i)\Sigma_{x_{0}}\Phi^{T}(0,j)\right]H^{T}_{j}.\end{split}$
(15)

Now, if the system is modeled with process noise, Eqs. 12,14, 15 becomes

$\mu_{y_{k}}=H_{k}\left[\Phi(0,k)\mu_{x_{0}}+\sum_{i=1}^{k}\Phi(k+1-i,k)G_{k-1}u_{k-1}\right],$
(16)

$\mathbb{E}\left[(y_{i}-\mu_{i})(y_{i}-\mu_{i})^{T}\right]=R_{i}+H_{i}\left[[\Phi(0,i)\Sigma_{x_{0}}\Phi^{T}(0,i)+\sum_{d=0}^{i-1}\Phi(d+1,i)\Gamma_{d}Q_{d}\Gamma_{d}^{T}\Phi^{T}(d+1,i)\right]H_{i}^{T},$
(17)

$\mathbb{E}\left[(y_{i}-\mu_{i})(y_{j}-\mu_{j})^{T}\right]=H_{i}\left[\Phi(0,i)\Sigma_{x_{0}}\Phi^{T}(0,j)+\sum_{d=0}^{m=min(i-1,j-1)}\Phi(d+1,i)\Gamma_{d}Q_{d}\Gamma_{d}^{T}\Phi^{T}(d+1,j)\right]H^{T}_{j}.$
(18)

With this derivation, we are able to analytically compute the SI for dynamic
systems and apply it to validate complex models.

##  4 Example Problem - Spacecraft Maneuver

A surprise index can be useful for detecting off-nominal behavior from the
system. While previous works on SI rely on empirical/sampling-based methods to
approximate the SI, which tends to be computationally intractable, we have
shown here that when the measurement model follows a multivariate Gaussian
probabilistic model, we can find a closed-form solution for SI.

We now show how SI can be used to analyze the decision-making or policy of an
RL-based model. More specifically, we test it for spacecraft stationkeeping
scenarios around halo orbits in the Sun-Earth system. In these scenarios, the
RL policy is trained to correct any deviation caused by perturbations and
other uncertainties, e.g., observation and navigation errors [9]. The RL
policy is trained in the circular restricted three-body problem (CR3BP); for
details on the dynamical model, training process, and required parameters, see
[9]. We simulate spacecraft trajectories, operating under policy $\pi_{2}$,
where the input to the policy is a noisy measurement of the state, and the
sensor uncertainty is accounted for in the training process. The trajectories
are perturbed around the working point, where the agent was trained vs.
further from this working point, to test the robustness of the policy. In
reality, these perturbations may occur due to imperfect models/data used to
train the RL policy or due to other errors and uncertainties that might occur
during the actual execution of the maneuvers.

SI can be used to evaluate these policies by leveraging the generalized
dynamic formulation described in Section 3.2. As in [9], we assume that the
system has no process noise and that we have a sensor with full observability
of the spacecraft. The (linear) measurement model is defined according to eq.
(10), with measurement matrix $H=\mathcal{I}_{6\times 6}$. With these
assumptions, equations (17)-(18) become,

$\mathbb{E}\left[(y_{i}-\mu_{i})(y_{i}-\mu_{i})^{T}\right]=\Sigma_{y_{i}}=R_{i}+\Phi(0,i)\Sigma_{x_{0}}\Phi^{T}(0,i)=R_{i}+\Sigma_{x_{i}},$
(19)

$\mathbb{E}\left[(y_{i}-\mu_{i})(y_{j}-\mu_{j})^{T}\right]=\Sigma_{y_{i},y_{j}}=\Phi(0,i)\Sigma_{x_{0}}\Phi^{T}(0,j)=\Sigma_{x_{i},x_{j}}.$
(20)

The challenge, in this case, is that the system’s dynamics are highly
nonlinear, thus, a closed-form expression for the state transition matrix
$\Phi(0,i)$ does not exist. Instead, the dynamics of the system are defined by
the nonlinear function $f(x,u)$, and while an approximation for $\Phi(0,i)$
can be calculated through the Jacobian of $f(x,u)$, it provides a poor
approximation for uncertainty propagation. Another method to propagate
uncertainty through nonlinear systems is by the sigma point or unscented
transformation [12].

###  4.1 Uncertainty Propagation via Sigma Points

The sigma point transformation allows for approximating the first two moments
(mean and covariance) of the distribution by propagating $2N+1$ ‘sigma points’
through the nonlinear dynamics, where $N$ is the number of states. After
propagating the points, the new mean and covariance can be calculated using,

$\begin{split}\mu_{x}&=\sum_{n=0}^{2N}\alpha_{n}\chi_{n}\\\
\Sigma_{x}&=\sum_{n=0}^{2N}\alpha_{n}(\chi_{n}-\mu_{x})(\chi_{n}-\mu_{x})^{T},\end{split}$
(21)

where $\chi_{n}$ is the $n$-th sigma point, and $\alpha_{n}$ is a weighting
factor, such that $\sum\alpha_{n}=1$. For more details on calculating the
sigma points and the weighting factor, see [12]. While the above equation can
be used to calculate the covariance in eq. (17), this accounts for only the
main diagonal blocks of the joint pdf. To calculate the off-diagonal blocks or
correlations between measurements taken in different time steps,
$\Sigma_{y_{i},y_{j}}$, according to eq. (18), we compute the following,

$\Sigma_{x_{i},x_{j}}=\sum_{n=0}^{2N}\alpha_{n}(\chi_{n_{i}}-\mu_{x_{i}})(\chi_{n_{j}}-\mu_{x_{j}})^{T}.$
(22)

Here $\chi_{n_{i}}$ ($\chi_{n_{j}}$) and $\mu_{i}$ ($\mu_{j}$) are the sigma
points and means at time step $i$ ($j$), respectively.

###  4.2 Test Case 1 - Nominal Conditions

To verify that the SI calculation method makes sense in the case of a
nonlinear dynamic system, we first apply it to a nominal case. In this
scenario, the spacecraft exactly follows the reference trajectory depicted in
Fig. 3.

Figure 3: (a) Reference orbit, with red indicating the segment of the orbit
that the SI was calculated over. (b) SI results for different true sensor
uncertainty values relative to the assumed sensor model.

We then perform 5 different 100 Monte-Carlo (MC) simulations, where we change
the sensor model statistics $R$ for each batch of simulations. The presented
SI is the cumulative SI, based on all measurements from time 0 to time step
$k$. The results demonstrate the expected, intuitive SI for the different
models: when the ‘true’ sensor is better, i.e., has a lower $R$ value, than
the assumed sensor model, as in the case of $R=0.5R_{model}$ and
$R=0.75R_{model}$, SI quickly converges to 1, meaning that the results are not
surprising. On the other hand, when the true sensor is worse than the assumed
model, as in the case of $R=1.5R_{model}$ and $R=2R_{model}$, the system sees
surprising measurements, and SI goes to zero. The most interesting result is
when the true sensor and the model agree, $R=R_{model}$, as shown by the green
line. In this case, the distribution for the SI value at any given time step
should be uniform on the range $[0,1]$. Thus, the mean should be $0.5$, as
achieved from the 100 MC simulations.

###  4.3 Test Case 2 - RL Policy Evaluation

In this test case, we evaluate the SI using the RL policy $\pi_{2}$ developed
in [9]. We perform 3 sets of 100 MC simulations, using sensor statistics used
for training and sampling an initial state in simulation. In each simulation,
the initial state is sampled from a normal distribution with a mean equal to
the nominal trajectory and a different covariance for each set of simulations.
After the initial state is sampled, a noisy measurement $y_{0}$ is taken
according to eq. (10). The noisy measurement is given as input to the RL
policy $\pi_{2}$, which provides an action (a change in velocity). This
process of taking a noisy measurement and using it as input to the policy
$\pi_{2}$, repeats 15 times along one time period of the halo orbit.

Results of the MC simulations are given in Fig. 4. Where Fig. 4(a) uses a
covariance equal to the trained initial covariance, $P_{0}=P_{trained}$.
Figures 4(b)-(c) use a larger initial covariance of $3P_{0}$ and $5P_{0}$,
respectively, i.e., sampling a more deviated initial state, outside of
distribution of the trained policy. The upper row of Fig. 4 shows the
reference (nominal) orbit in black and 100 MC trajectories in blue. We can see
that in all cases, the trajectory follows the reference orbit for more than
$2/3$ of the period. The SI in the lower row of the plots shows the mean (in
blue) and $1\sigma$ and $2\sigma$ regions of the 100 simulations in shaded
gray and pink, respectively. Since evaluating the SI at every time step has a
significant computational burden due to inversion of the joint covariance
matrix (eq. (7), the SI was evaluated using 100 measurements equally spaced
along the trajectory. The SI results show that when the trajectory aligns with
the expected orbit, the mean is about 0.5, as we saw in the nominal case
(green line in Fig. 3). However, when the trajectories start to diverge, the
SI value drops, which happens earlier when the initial uncertainty is higher.

Figure 5 shows two of the 100 trajectories in Fig. 4(a), i.e. with initial
uncertainty $P_{0}$. We can see that the SI value indicates deviation from the
reference trajectory, where the green trajectory follows the orbit most of the
time, with $SI>0.8$, and when it starts diverging at the end, the SI plummets.
The red trajectory starts diverging much earlier, which is also shown by the
values of the SI.

Figure 4: Trajectories (upper row) and SI evaluation (lower row) for
stationkeeping scenario with different initial uncertainty. (a) Initial
uncertainty is equivalent to what the policy was trained on $P_{0}$. (b) and
(c): $3$ and $5$ times the trained initial uncertainty, respectively.  Figure
5: Trajectory comparison – comparing good performance (green trajectory) to a
deviated trajectory (red). (a) shows the trajectories vs. the reference
(black) stationkeeping orbit. (b) shows the corresponding SI.

##  5 Summary

As autonomous systems’ decision-making relies more on ‘black-box’ policies
generated by a learning algorithm, the surprise index provides a possible
avenue to assess the system’s competence. That is, assess whether the
performance of the autonomous system agrees with what we expect it to do. In
this paper, we derived a closed-form solution to calculate the SI in the case
of multivariate Gaussian measurement models. We applied it to the nonlinear
system of the spacecraft stationkeeping maneuver and showed that it provides a
good measure of the spacecraft’s performance against a reference model of its
expected behavior.

However, computing the SI value for large systems or over a long time horizon
might be computationally intractable since it requires inverting the joint
covariance matrix over all measurements. Thus future work will explore
approximation and/or efficient calculation methods for SI.

## Acknowledgments

This work is supported by an Early Stage Innovations grant from NASA’s Space
Technology Research Grants Program, under NASA grant 80NSSC19K0222.

## References

  * Aitken [2014] Aitken, M., “Assured Human-Autonomy Interaction through Machine Self-Confidence,” Master’s thesis, University of Colorado at Boulder, 2014.  ISBN: 9781369177145.
  * Israelsen and Ahmed [2022] Israelsen, B. W., and Ahmed, N., “A Factor-Based Framework for Decision-Making Competency Self-Assessment,” , Mar. 2022.  URL http://arxiv.org/abs/2203.11981, arXiv:2203.11981 [cs, eess].
  * Conlon et al. [2023] Conlon, N., Ahmed, N. R., and Szafir, D., “Dynamic Competency Self-Assessment for Autonomous Agents,” _ACM/IEEE HRI 2023 workshop on Variable Autonomy for human-robot Teaming (VAT-HRI)_ , arXiv, 2023.  10.48550/arXiv.2303.01646, arXiv:2303.01646 [cs].
  * Habbema [1976] Habbema, J., “Models for diagnosis and detection of combinations of diseases,” _Decision Making and Medical Care_ , 1976, pp. 399–411.
  * Bar-Shalom et al. [2001] Bar-Shalom, Y., Li, X. R., and kirubarajan, T., “Linear estimation in static systems,” _Estimation with Applications to Tracking and Navigation_ , John Wiley & Sons, Ltd, 2001, pp. 121–177.
  * Zagorecki et al. [2015] Zagorecki, A., Kozniewski, M., and Druzdzel, M. J., “An Approximation of Surprise Index as a Measure of Confidence,” 2015, p. 3.
  * Laskey [1991] Laskey, K. B., “Conflict and Surprise: Heuristics for Model Revision,” _Uncertainty Proceedings 1991_ , edited by B. D. D’Ambrosio, P. Smets, and P. P. Bonissone, Morgan Kaufmann, San Francisco (CA), 1991, pp. 197–204.  10.1016/B978-1-55860-203-8.50029-2, URL https://www.sciencedirect.com/science/article/pii/B9781558602038500292.
  * Bonasera et al. [2023] Bonasera, S., Bosanac, N., Sullivan, C. J., Elliott, I., Ahmed, N., and McMahon, J. W., “Designing Sun–Earth L2 Halo Orbit Stationkeeping Maneuvers via Reinforcement Learning,” _Journal of Guidance, Control, and Dynamics_ , Vol. 46, No. 2, 2023, pp. 301–311.  10.2514/1.G006783, URL https://doi.org/10.2514/1.G006783, publisher: American Institute of Aeronautics and Astronautics _eprint: https://doi.org/10.2514/1.G006783.
  * Boone et al. [2022] Boone, S., Bonasera, S., McMahon, J. W., Bosanac, N., and Ahmed, N. R., “Incorporating Observation Uncertainty into Reinforcement Learning-Based Spacecraft Guidance Schemes,” _AIAA SCITECH 2022 Forum_ , American Institute of Aeronautics and Astronautics, San Diego, CA & Virtual, 2022.  10.2514/6.2022-1765, URL https://arc.aiaa.org/doi/10.2514/6.2022-1765.
  * Jensen et al. [1990] Jensen, F. V., Chamberlain, B., Jensen, F., and Nordahl, T., “Analysis in HUGIN of Data Conflict,” _Proceedings of the Sixth Conference on Uncertainty in Artificial Intelligence_ , arXiv, 1990.  10.48550/arXiv.1304.1146, URL http://arxiv.org/abs/1304.1146, arXiv:1304.1146 [cs].
  * Bar-Shalom et al. [2009] Bar-Shalom, Y., Daum, F., and Huang, J., “The probabilistic data association filter,” _IEEE Control Systems Magazine_ , Vol. 29, No. 6, 2009, pp. 82–100.  10.1109/MCS.2009.934469, URL https://ieeexplore.ieee.org/document/5338565, conference Name: IEEE Control Systems Magazine.
  * Julier et al. [1995] Julier, S. J., Uhlmann, J. K., and Durrant-Whyte, H. F., “A new approach for filtering nonlinear systems,” _Proceedings of 1995 American Control Conference - ACC’95_ , Vol. 3, 1995, pp. 1628–1632 vol.3.  10.1109/ACC.1995.529783.