# Ren-integrable and ren-symmetric integrable systems

S. Y. Lou  School of Physical Science and Technology, Ningbo University,
Ningbo, 315211, China

###### Abstract

A new type of symmetry, ren-symmetry describing anyon physics and the
corresponding topological physics, is proposed. Ren-symmetry is a
generalization of super-symmetry which is widely applied in super-symmetric
physics such as the super-symmetric quantum mechanics, super-symmetric
gravity, super-symmetric string theory, super-symmetric integrable systems and
so on. The super-symmetry and Grassmann-number are, in some sense, the dual
conceptions, which turns out that these conceptions coincide for the ren
situation, that is, a similar conception of ren-number is devised to ren-
symmetry. In particular, some basic results of the ren-number and ren-symmetry
are exposed which allow one to derive, in principle, some new types of
integrable systems including ren-integrable models and ren-symmetric
integrable systems. Training examples of ren-integrable KdV type systems and
ren-symmetric KdV equations are explicitly given.

##  I Introduction

The idea of symmetry originates in natural scientific fields, and its
importance there is well known. Symmetry considerations belong to the most
universal and astonishing methods by which scientists have successfully solved
the problems in building new solutions from known ones Sol , doing dimensional
reductions of nonlinear partial differential equations Sol1 ; Sol2 ; Sol3 ;
Sol4 , getting new integrable systems Sol5 ; Sol6 ; Sol7 ; Sol8  and even
constructing all solutions for certain nonlinear systems Sol9 .

By using the SU(3)$\times$SU(2)$\times$U(1) symmetry, three fundamental
interactions, the strong, weak and electromagnetic interactions, have been
unified into the so-called standard model. However, in order to unify the
gravitational interaction to the standard model, one has to introduce a new
type of symmetry, say the super-symmetry between bosons and fermions. New
areas of physical fields including the super-symmetric gravity SG ; SG1 ,
super-symmetric quantum mechanics SM , super-symmetric string theory SS  and
super-symmetric integrable systems SusyI ; SusyII ; SusyIII ; SusyIV ; SusyV ;
SusyVI  have been developed which are highly motivated by super-symmetries, in
the belief that they possess a high potential for future development.

In super-symmetry theory, it is essential to introduce the Grassmann variable
$\theta$ GA ; GA1  and the super-symmetric derivative ${\cal{D}}$ with the
properties

$\displaystyle\theta^{2}=0,\ \theta_{1}\theta_{2}=-\theta_{2}\theta_{1},$ (1)
$\displaystyle{\cal{D}}=\partial_{\theta}+\theta\partial_{x},\
{\cal{D}}^{2}=\partial_{x}.$ (2)

The super-symmetric derivative ${\cal{D}}$ is invariant under the super-
symmetric transformation

$\displaystyle\theta\rightarrow\theta+\eta,\ x\rightarrow x-\theta\eta.$ (3)

Recently, different to the bosons and fermions, anyons with fractional
charges, spin and statistics in two dimensions have been attracted high
sufficient attentions by many scientists any1 ; any2 ; any3 ; any4 . Anyons
can be used to describe some kinds of quasi-particles (the low-energy
excitations in Hamiltonian systems) including the fractional quantum Hall
states qHall , vortices in topological superconductors topo  and Majorana zero
modes in semiconductors proximitized by superconductors Majo . By analogy with
fermion case in which fermions can be described by Grassmann fields, some new
fields endowed to describe anyons, we call anyon-fields and/or ren-fields. We
shall use the adjective “ren” to stress the arbitrary of $\alpha$ and to avoid
the confusion on “arbitrary symmetry” or “any symmetry”. “Ren” means
“arbitrary” in Chinese.

A comparison of the Grassmann number $\theta=\sqrt{0}$ and the super-symmetric
derivative ${\cal{D}}=\sqrt{\partial_{x}}$, corresponding to ren point of
view, suggests one to use, as the ren-number and the ren-symmetric derivative,
respectively, the following radical generalization of the formulae

$\theta_{\alpha}\equiv\theta=\sqrt[\alpha]{0},~{}~{}{\cal{R}}=\sqrt[\alpha]{\partial_{x}}$

with $\alpha$ being arbitrary.

The introduction of the G-number and the super-symmetric derivative yields
some significant novel mathematical and physical fields such as the Grassmann
algebra GA , the super-symmetric quantum mechanics SM , the super-symmetric
string theory SS , the super-symmetric gravity SG , the super/Kuper-integrable
systems SI ; SI1 ; SI2 ; SI3  and super-symmetric integrable theories SusyI .
Therefore, we hope that the introduction of the ren-number and the ren-
symmetric derivative, may successfully create some different mathematical and
physical fields such as the ren-algebra, ren-calculus, ren integrable models
and ren-symmetric integrable systems. The usual G-number, G-algebra, super-
symmetric theory, super-integrable systems and super-symmetric integrable
systems just correspond to the ren-case for $\alpha=2$.

In Sec. II of this paper, the concept of the R-number $\theta$ for an
arbitrary positive integer $\alpha$ is defined with the aim of deriving ren-
algebra, ren-derivative, and ren-symmetric derivative. Then, we deal in Sec.
III with the problem of finding some types of ren-integrable systems by
coupling the usual boson fields and the anyon fields. When $\alpha$ is fixed
to $\alpha=2$, the ren-integrable system is just the known super- or Kuper-
integrable systems SI . By means of the ren-symmetric derivative, we study the
ren-symmetric integrable systems in Sec. IV, the ren-symmetric KdV systems for
$\alpha=3,\ 4$ are explicitly given. The well known super-symmetric integrable
systems are just the special cases of the ren-symmetric integrable systems
with $\alpha=2$. The last section includes a short summary and some
discussions.

##  II Ren-algebra, ren-derivative and ren-symmetric derivative

Definition 1. A ren-number $\theta\equiv\theta_{\alpha}$ is defined as a
number possessing the properties

$\theta^{\alpha}=0,\ \theta^{i}\neq 0,\ i=1,\ 2,\ \ldots,\ \alpha-1,$ (4)

where $\alpha$ is an arbitrary positive integer.

This is, a ren-number can be spelled out as a non-zero $\alpha$ root of zero

$\theta=\sqrt[\alpha]{0}\neq 0.$ (5)

It is clear that there are $\alpha-1$ solutions of (5), $\\{\theta,\
\theta^{2},\ \ldots,\ \theta^{\alpha-1}\\}$. Such a definition, exists in one
important special case $\alpha=2$, the usual G-number $\theta=\theta_{2}$.

For $\alpha=2$, we know that if $a$ and $b$ are Grassmann numbers then the
combination is still a G-number when the anti-communication relation

$ab=-ba$ (6)

holds.

Similarly, for $\alpha\geq 2$, if $a$ and $b$ are ren-numbers with the
$q$-commutation relation ($\mbox{\rm i}\equiv\sqrt{-1}$),

$\displaystyle ab=q_{j}ba,\quad q_{j}^{\alpha}=1,\ q_{j}=q^{j},\
q=\exp\left(\frac{2\pi\mbox{\rm i}}{\alpha}\right),$ $\displaystyle\qquad
j\in{\cal{J}_{\alpha}}=\\{1,\ 2,\ \ldots,\
\alpha-1\\}/\\{n_{m}p_{m}<\alpha\\},$ (7)

then so is the combination $a+b$ of $a$ and $b$, where
$\\{n_{m}p_{m}<\alpha\\}$ is a set with $p_{m},\ m=1,\ \ldots,\ M$ being the
prime factors of $\alpha$ and $n_{m}=1,\ \ldots,\ N_{m}$ being integers with
$N_{m}p_{m}<\alpha$. For $\alpha=3,\ 4,\ \ldots,\ 10$, we have

$\displaystyle{\cal{J}}_{3}=\\{1,\ 2\\},\ {\cal{J}}_{4}=\\{1,\ 3\\},\
{\cal{J}}_{5}=\\{1,\ 2,\ 3,\ 4\\},\ {\cal{J}}_{6}=\\{1,\ 5\\},\
{\cal{J}}_{7}=\\{1,\ 2,\ \ldots,\ 6\\},\ $ $\displaystyle{\cal{J}}_{8}=\\{1,\
3,\ 5,\ 7\\},\ {\cal{J}}_{9}=\\{1,\ 2,\ 4,\ 5,\ 7,\ 8\\},\
{\cal{J}}_{10}=\\{1,\ 3,\ 7,\ 9\\}.$ (8)

From the expression of $q$ given in (II), we know that the usual number (boson
number) is related to $\alpha=\infty$ and the Grassmann number (fermion
number) corresponds to $\alpha=2$.

Definition 2. The degree, $\beta\ (\mbox{\rm mod}(\alpha))$, of a ren-number
$\gamma$ is defined as

$\gamma\theta=q^{\beta}\theta\gamma,\\\ $ (9)

where the degree of $\theta$ is always fixed as one in this paper. $\gamma$
with (9) is said also a $\beta$ order ren-number.

If the ren-numbers $\gamma_{1}$ and $\gamma_{2}$ possess the degrees
$\beta_{1}$ and $\beta_{2}$, respectively, then we have the commutation
relation

$\gamma_{1}\gamma_{2}=q^{\beta_{1}\beta_{2}}\gamma_{2}\gamma_{1},\\\ $ (10)

which is consistent with (9) when $\beta_{2}=1$ and $\beta_{1}\beta_{2}=0$.

Definition 3. Ren-derivative, $\frac{\mbox{\rm d}}{\mbox{\rm d}\theta}$, is a
derivative with respect to the ren-variable $\theta$,

$\frac{\mbox{\rm d}f(\theta)}{\mbox{\rm
d}\theta}=\left.\frac{f(\theta_{1})-f(\theta)}{\theta_{1}-\theta}\right|_{\theta_{1}\rightarrow\theta}=\frac{f(\theta)-f(q\theta)}{(1-q)\theta}.$
(11)

Similar to the Grassmann case, because the definition of the ren-number (4),
an arbitrary function of ren-variable, $f(\theta)$ can be written as

$f(\theta)=\sum_{i=0}^{\alpha-1}\theta^{i}g_{i}=\sum_{i=0}^{\alpha-1}f_{i}\theta^{i}.$
(12)

If $f(\theta)$ is a $\beta$ order ren-number, then $f_{i}$ and $g_{i}$ in (12)
are $\beta-i$ order ren-numbers with $f_{i}=q^{i\beta-i^{2}}g_{i}$.

According to the property (12), it is enough to find all the possible ren-
derivatives for an arbitrary ren-function $f(\theta)$ by calculating

$\frac{\mbox{\rm d}\theta^{i}}{\mbox{\rm d}\theta},\ i=1,\ 2,\ \ldots,\
\alpha-1.$

Based on the commutation relation (II) and the definition of the ren-
derivative (11), it is readily to prove that

$\frac{\mbox{\rm d}\theta^{i}}{\mbox{\rm
d}\theta}=\sum_{k=0}^{i-1}q^{k}\theta^{i-1}=\frac{1-q^{i}}{1-q}\theta^{i-1}\equiv
i_{q}\theta^{i-1},$ (13)

where $i_{q}$ is defined as $i_{q}=1+q+\cdots+q^{i-1}$, say, $2_{q}=1+q,\
3_{q}=1+q+q^{2}$ and so on.

Thus, for the ren-function $f(\theta)$ with degree $\beta$, we have

$\frac{\mbox{\rm d}f(\theta)}{\mbox{\rm
d}\theta}=\sum_{i=0}^{\alpha-1}i_{q}\theta^{i-1}g_{i}=\sum_{i=0}^{\alpha-1}\frac{(q^{\beta-i}-q^{\beta})}{1-q}f_{i}\theta^{i-1}=\sum_{i=0}^{\alpha-1}q^{\beta-i}i_{q}f_{i}\theta^{i-1}.$
(14)

The ren-integration may be defined as an inverse operator of the ren-
derivative for $\theta^{k},\ k<\alpha-1$, however, for $\theta^{\alpha-1}$ the
inverse operator of the ren-derivative is not well defined. A different
integration operator can be defined Majid . For $\alpha=2$, the Berezin
integral has been defined Berezin ; Bin . In this paper, we will not discuss
this problem though the similar Berezin integral may be introduced under the
requirement of the translation invariance Majid .

Definition 4. A ren-symmetric derivative $\cal{R}\equiv\cal{R}_{\alpha}$ is
defined as an $\alpha$ root of the usual space derivative $\partial_{x}$,
i.e.,

${\cal{R}^{\alpha}}=\partial_{x},\ {\cal{R}}=\sqrt[\alpha]{\partial_{x}}.$
(15)

It is interesting that in terms of the ren-number $\theta$, the ren-symmetric
derivative $\cal{R}$ can be explicitly written as

${\cal{R}}=\partial_{\theta}+\frac{1}{[(\alpha-1)!]_{q}}\theta^{\alpha-1}\partial_{x},$
(16)

where $[n!]_{q}$ is defined as

$\displaystyle[n!]_{q}$ $\displaystyle\equiv$
$\displaystyle\prod_{i=1}^{n}\frac{1-q^{i}}{1-q}\equiv\prod_{i=1}^{n}i_{q}\equiv
1_{q}\cdot 2_{q}\cdots(n-1)_{q}\cdot n_{q},$ (17)

whose particular case $q=1$ is the usual $n!$.

It is reasonable that the ren-symmetric derivative (16) will reduce back to
the known super-symmetric derivative ${\cal R}_{2}\equiv{\cal
D}=\partial_{\theta}+\theta\partial_{x}$ when $\alpha=2$. The ren-symmetric
derivatives for $\alpha=3,\ 4,\ 5,\ 6$ and $7$ are given by the formulae as a
straightforward computation,

$\displaystyle{\cal
R}_{3}=\partial_{\theta}+\frac{1}{[2!]_{q}}\theta^{2}\partial_{x}=\partial_{\theta}-q\theta^{2}\partial_{x},\
q=\mbox{\rm e}^{\left({2\pi\mbox{\rm\tiny
i}}/3\right)}=\frac{1}{2}(\sqrt{3}\mbox{\rm i}-1),$ $\displaystyle{\cal
R}_{4}=\partial_{\theta}+\frac{1}{[3!]_{q}}\theta^{3}\partial_{x}=\partial_{\theta}-\frac{1}{2}(1+q)\theta^{3}\partial_{x},\
q=\mbox{\rm e}^{\left({2\pi\mbox{\rm\tiny i}}/4\right)}=\mbox{\rm i},$
$\displaystyle{\cal
R}_{5}=\partial_{\theta}+\frac{1}{[4!]_{q}}\theta^{4}\partial_{x}=\partial_{\theta}+\frac{q^{3}}{(1+q)^{2}}\theta^{4}\partial_{x},\
q=\mbox{\rm e}^{\left({2\pi\mbox{\rm\tiny i}}/5\right)},$ $\displaystyle{\cal
R}_{6}=\partial_{\theta}+\frac{1}{[5!]_{q}}\theta^{5}\partial_{x}=\partial_{\theta}+\frac{q}{6}\theta^{5}\partial_{x},\
q=\mbox{\rm e}^{\left({2\pi\mbox{\rm\tiny
i}}/6\right)}=\frac{1}{2}(1+\sqrt{3}\mbox{\rm i}),$ $\displaystyle{\cal
R}_{7}=\partial_{\theta}+\frac{1}{[6!]_{q}}\theta^{6}\partial_{x}=\partial_{\theta}-\frac{q^{6}}{[3!]_{q}^{2}}\theta^{6}\partial_{x},\
q=\mbox{\rm e}^{\left({2\pi\mbox{\rm\tiny i}}/7\right)}.$ (18)

It is not difficult to prove that the ren-symmetric derivative ${\cal{R}}$
possesses the following ren-symmetric transformation

$\theta\rightarrow\theta+\eta,\ x\rightarrow x-f(\theta,\ \eta),$ (19)

with

$f=f(\theta,\
\eta)=\sum_{k=1}^{\alpha-1}\frac{1}{[k!]_{q}[(\alpha-k)!]_{q}}\theta^{k}\eta^{\alpha-k}.$
(20)

With the stress on the first few $f$, say, $\alpha=2,\ 3,\ 4,\ 5$ and $6$, we
have

$\displaystyle f=\theta\eta,$ $\displaystyle\alpha=2,$ (21) $\displaystyle
f=\frac{1}{[2!]_{q}}(\theta\eta^{2}+\theta^{2}\eta),$ $\displaystyle\alpha=3,$
$\displaystyle
f=\frac{1}{[3!]_{q}}\left(\theta\eta^{3}+\frac{3_{q}}{2_{q}}\theta^{2}\eta^{2}+\theta^{3}\eta\right),$
$\displaystyle\alpha=4,$ $\displaystyle
f=\frac{1}{[4!]_{q}}\left(\theta\eta^{4}+\frac{4_{q}}{2_{q}}\theta^{2}\eta^{3}+\frac{4_{q}}{2_{q}}\theta^{3}\eta^{2}+\theta^{4}\eta\right),$
$\displaystyle\alpha=5,$ $\displaystyle
f=\frac{1}{[5!]_{q}}\left(\theta\eta^{5}+\frac{5_{q}}{2_{q}}\theta^{2}\eta^{4}+\frac{5_{q}4_{q}}{3_{q}2_{q}}\theta^{3}\eta^{3}+\frac{5_{q}}{2_{q}}\theta^{4}\eta^{2}+\theta^{5}\eta\right),$
$\displaystyle\alpha=6.$

##  III Ren-integrable systems

In the limit $\alpha=2$, ren-integrable models are just the known super- or
Kuper- integrable models which are first proposed by Kupershmidt in SI . It is
known that the usual bosonic KdV equation,

$u_{t}=(-u_{xx}+3u^{2})_{x},$ (22)

is determined as a compatibility condition $[L,\ S]=LS-SL=0$ of equations from
a Lax pair

$\displaystyle\psi_{xx}-(u+\lambda)\psi\equiv L\psi=0,$
$\displaystyle\psi_{t}-3u_{x}\psi-6u\psi_{x}+4\psi_{xxx}\equiv S\psi=0.$ (23)

Usually, the spectral function $\psi$ is considered as a boson function. In
fact, because the Lax pair (23) is linear, the spectral function $\psi$ may be
a fermion function and even a ren function.

It is known that if $\sigma$ is a symmetry of an integrable evolution equation

$u_{t}=K(u),$ (24)

i.e., a solution of

$\sigma_{t}=K^{\prime}\sigma\equiv\lim_{\epsilon=0}\partial_{\epsilon}K(u+\epsilon\sigma),$
(25)

then

$u_{t}=K(u)+\sigma,$ (26)

is also an integrable model.

Furthermore, if $\sigma=\sigma(\psi)$, where $\psi$ is a spectral function of
the Lax pair,

$\hat{L}\psi=0,\qquad\hat{S}\psi=0,$ (27)

of (24), then the first type of source equation

$\displaystyle u_{t}=K(u)+\sigma(\psi),$ $\displaystyle\hat{L}\psi=0,$ (28)

and the second type of source equation

$\displaystyle u_{t}=K(u)+\sigma(\psi),$ $\displaystyle\hat{S}\psi=0,$ (29)

may all be integrable Source1 ; Source2 ; Source3 ; Source4 ; Source5 ;
Source6 ; Source7 ; Source8 ; Source9 .

Usually, the spectral functions $\psi$ studied in the integrable models are
restricted as bosonic functions. For instance, for the KdV equation (22) the
first and second types of integrable bosonic source equations possess the
forms

$\displaystyle
u_{t}=(-u_{xx}+3u^{2}+\langle\phi|\phi\rangle)_{x},\qquad\langle\phi|\phi\rangle\equiv\sum_{i=1}^{n}\phi_{i}^{2},$
$\displaystyle(\partial_{x}^{2}-u-\lambda_{i})\phi_{i}=0,\quad i=1,\ 2,\
\ldots,\ n,$ (30)

and

$\displaystyle u_{t}=(-u_{xx}+3u^{2}+\langle\phi|\phi\rangle)_{x},$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})\phi_{i}=0,\quad
i=1,\ 2,\ \ldots,\ n,$ (31)

respectively.

Now, if we extend the spectral function of the Lax pair (23) to a ren-
function, $\xi$, then we have a trivial symmetry,
$\sigma=\xi^{\alpha-1}\xi_{xx}$. Applying this symmetry to the second type of
source equation, we can find some coupled ren systems

$\displaystyle
u_{t}=(-u_{xx}+3u^{2})_{x}+\sum_{i=1}^{k}\langle\xi_{i}^{\alpha_{i}-1}|\xi_{i,xx}\rangle,\qquad\langle\xi_{i}^{\alpha_{i}-1}|\xi_{i,xx}\rangle\equiv\sum_{j=1}^{n_{i}}\xi_{ij}^{\alpha_{i}-1}\xi_{ij,xx},$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})\xi_{ij}=0,\quad
j=1,\ 2,\ \ldots,\ n_{i},\ i=1,\ 2,\ \ldots,\ k,\ \xi_{ij}^{\alpha_{i}}=0$
(32)

with one boson field $u$ and $k\times n_{i}$ ren fields $\xi_{ij}$, where $k,\
n_{i}$ and $\alpha_{i}$ are all arbitrary integers.

The integrability of (32) with $\alpha_{i}=2$ for all $i=1,\ 2,\ \ldots,\ k$
is known because the models reduce back to the so-called super/Kuper-
integrable systems SI ; SI1 ; SI2 ; SI3 . Before studying the integrability of
(32) with $\alpha_{i}\neq 2$ for some $2<i\leq k$, we directly write down a
more general nontrivial symmetry of the KdV equation (22),

$\displaystyle\sigma=(\xi_{1}\xi_{2x}-\xi_{1x}\xi_{2})_{x},$ (33)

where $\xi_{1}$ and $\xi_{2}$ are ren-spectral functions of the usual KdV
equation with the same spectral parameter $\lambda=\lambda_{1}=\lambda_{2}$
but with different degrees, $\beta$ and $\alpha-\beta$, respectively.

The simplest second type of source equation related to (33) possesses the form

$\displaystyle u_{t}=[3u^{2}-u_{xx}+12(\xi_{1}\xi_{2x}-\xi_{1x}\xi_{2})]_{x},$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})\xi_{i}=0,\
i=1,\ 2.$ (34)

Theorem. The model (34) is Lax integrable with the Lax pair

$\displaystyle\psi_{xx}-(u+\lambda)\psi+\left(\int\xi_{1}\psi\mbox{\rm
d}x\right)\xi_{2}-\xi_{1}\left(\int\xi_{2}\psi\mbox{\rm
d}x\right)\equiv\hat{L}\psi=0,$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})\psi\equiv\hat{S}\psi=0.$
(35)

Proof. To complete the proof of the theorem, it suffices to prove that the
compatibility condition

$\displaystyle[\hat{L},\ \hat{S}]f=(\hat{L}\hat{S}-\hat{S}\hat{L})f$ (36)

is valid for arbitrary $f$ if (34) is satisfied.

Expanding the expression (36) with the operators $\hat{L}$ and $\hat{S}$
defined in (35),

$\displaystyle\left(\int\xi_{1}f\mbox{\rm
d}x\right)(6u\xi_{2x}+3\xi_{2}u_{x}-4\xi_{2xxx}-\xi_{2t})+(\xi_{1t}+4\xi_{1xxx}-6u\xi_{1x}-3\xi_{1}u_{x})\int\xi_{2}f\mbox{\rm
d}x$
$\displaystyle\qquad+\xi_{1}\left[\int(\xi_{2t}f-4\xi_{2}f_{xxx}+6\xi_{2}uf_{x}+3\xi_{2}u_{x}f)\mbox{\rm
d}x-8\xi_{2xx}f-4\xi_{2x}f_{x}\right]$
$\displaystyle\qquad-\left[\int(\xi_{1t}f-4\xi_{1}f_{xxx}+6\xi_{1}uf_{x}+3\xi_{1}u_{x}f)\mbox{\rm
d}x-8\xi_{1xx}f-4\xi_{1x}f_{x}\right]\xi_{2}$
$\displaystyle\qquad+(u_{t}-6uu_{x}+u_{xxx})f=0,$ (37)

and simplifying the result with the formulae of integration by parts

$\displaystyle\int\xi_{i}f_{xxx}\mbox{\rm
d}x=\xi_{i}f_{xx}-\xi_{ix}f_{x}+\xi_{ixx}f-\int\xi_{ixxx}f\mbox{\rm d}x,$
$\displaystyle\int\xi_{i}uf_{x}\mbox{\rm d}x=\xi_{i}uf-\int
f(u\xi_{i})_{x}\mbox{\rm d}x,$ (38)

(37) is changed to

$\displaystyle\int\xi_{1}f\mbox{\rm d}x\
(6u\xi_{2x}+3\xi_{2}u_{x}-4\xi_{2xxx}-\xi_{2t})+(\xi_{1t}+4\xi_{1xxx}-6u\xi_{1x}-3\xi_{1}u_{x})\int\xi_{2}f\mbox{\rm
d}x$
$\displaystyle+\xi_{1}\int(\xi_{2t}+4\xi_{2xxx}-6u\xi_{2x}-3\xi_{2}u_{x})f\mbox{\rm
d}x-\left[\int(\xi_{1t}+4\xi_{1xxx}-6u\xi_{1}-3\xi_{1}u_{x})f\mbox{\rm
d}x\right]\xi_{2}$
$\displaystyle+(u_{t}-6uu_{x}+u_{xxx}+12\xi_{1xx}\xi_{2}-12\xi_{1}\xi_{2xx})f=0.$
(39)

Because of the arbitrariness of $f$, (39) is correct only after joining to it
the equation (34). The theorem is proved.

Remark. From the proof procedure of the theorem, it is known that we have not
used any commutation relation on $\xi_{1}$ and $\xi_{2}$. That means (34) is
Lax integrable no matter the fields $\xi_{1}$ and $\xi_{2}$ are boson fields,
fermion fields and/or ren-fields with arbitrary $\alpha$.

##  IV Ren-symmetric integrable systems

In Sec. II of this paper, we have defined the ren-symmetric derivative
$\cal{R}$. By means of the ren-symmetric derivative, the usual bosonic
integrable systems can be extended to ren-symmetric integrable ones. Before
discussing ren-symmetric integrable systems, we list some special cases for
$\alpha=2$, i.e., the super-symmetric integrable models.

###  IV.1 Super-symmetric integrable KdV systems

The most general $N=1$ symmetric form of the KdV equation (22) is generated by
the fermionic super-field $\Phi$ with an arbitrary constant $a$,

$\Phi_{t}+\Phi_{xxx}+a({\cal{D}}\Phi_{x})\Phi+(6-a)({\cal{D}}\Phi)\Phi_{x}=0.$
(40)

Mathieu had proven that the super-symmetric KdV equation (40) is Painlevé
integrable only for $a=0$ and $3$ Mathieu . Although the super-symmetric KdV
system (40) is not Painlevé integrable for arbitrary $a$, it does possess
multiple soliton solutions KdVa . In (40), the super-field
$\Phi\equiv\xi+\theta u$ is a fermionic super-field with a fermion field $\xi$
and a boson field $u$.

For the coupled KdV equation, we have an interacting model between a susy-
boson field $U$ and a susy-fermion field $\Phi$

$\displaystyle\Phi_{t}=(-\Phi_{xx}+3\Phi{\cal D}\Phi+6U\Phi)_{x},$
$\displaystyle U_{t}=(-U_{xx}+3U^{2}+3\Phi{\cal D}U)_{x},$ (41)

which is Lax integrable. Incorporating $U=0$, (41) readily reduces back to
(40) with $a=3$. For $\Phi=0$, (41) becomes a quite trivial extension of the
original KdV equation (22) by $u\rightarrow U$.

The component form of (41) reads

$\displaystyle v_{t}=(-v_{xx}+3v^{2}+3\zeta_{x}\zeta+6uv+6\xi\zeta)_{x},$
$\displaystyle u_{t}=(-u_{xx}+3u^{2}+3\zeta\xi)_{x},$
$\displaystyle\xi_{t}=(-\xi_{xx}+6u\xi+3v\xi+3u_{x}\zeta)_{x},$
$\displaystyle\zeta_{t}=(-\zeta_{xx}+6u\zeta+3v\zeta)_{x},$ (42)

with $U=u+\theta\xi$ and $\Phi=\zeta+\theta v$, where $u$ and $v$ are boson
components and $\xi$ and $\zeta$ are fermion components.

The Lax pair of (41) can be written as

$\displaystyle\Psi_{xx}=\Phi{\cal D}\Psi+(U+\lambda)\Psi,$
$\displaystyle\Psi_{t}=-4\Psi_{xxx}+6U\Psi_{x}+3U_{x}\Psi+6\Phi{\cal
D}\Psi_{x}+3\Phi_{x}{\cal D}\Psi.$ (43)

###  IV.2 Ren-symmetric integrable KdV systems

Analogous to (40), the general ren-symmetric KdV equation is expressible in
the form

$\Phi_{t}+\Phi_{xxx}+\sum_{i=0}^{\left[\beta_{1}\right]}a_{i}\left({\cal{R}}^{i}\Phi\right)\left({\cal{R}}^{\alpha+\beta-i}\Phi\right)=0,\
\beta=0,\ 1,\ 2,\ \ldots,\ \alpha-1,\ \beta_{1}\equiv\frac{\alpha+\beta}{2},$
(44)

where $\left[\beta_{1}\right]$ is the integer part of $\beta_{1}$, $a_{i},\
i=0,\ 1,\ 2,\ \ldots,\ [\beta_{1}]$, are arbitrary bosonic constants, $\beta$
is the degree of the ren-field $\Phi\equiv\Phi(x,\ t,\ \theta)$.

As in the super-symmetric ($\alpha=2$) case, one may find some possible
integrable cases by fixing the constants $a_{i}$ of the ren-symmetric KdV
equation (44). For instance, the Lax integrable systems,

$\Phi_{jt}+\Phi_{jxxx}-3{\cal{R}}^{\alpha-j}\left({\cal{R}}^{j}\Phi_{j}\right)^{2}+\rho_{j}=0,\
{\cal{R}}^{j}\rho_{j}=0,\ j=0,\ 1,\ \ldots,\ \alpha-1,\ $ (45)

with $\rho_{j}=0$ possessing Lax pair of

$\displaystyle\Psi_{xx}-({\cal{R}}^{j}\Phi_{j}+\lambda)\Psi=0,$
$\displaystyle\Psi_{t}+4\Psi_{xxx}-3\left({\cal{R}}^{j}\Phi_{jx}\right)\Psi-6\left({\cal{R}}^{j}\Phi_{j}\right)\Psi_{x}=0,$
(46)

are just the special cases of (44). The degrees of $\Phi_{j}$ and $\rho_{j}$
in (45) are $j$.

For $\alpha=3$, the ren-symmetric KdV equation (44) becomes

$\displaystyle\Phi_{0t}+\Phi_{0xxx}+a\Phi_{0}\Phi_{0x}+b({\cal{R}}\Phi_{0})({\cal{R}}^{2}\Phi_{0})=0,\
{\cal{R}}=\partial_{\theta}-q\theta^{2}\partial_{x},$ (47)
$\displaystyle\Phi_{1t}+\Phi_{1xxx}+a\Phi_{1}({\cal{R}}\Phi_{1x})+b({\cal{R}}\Phi_{1})\Phi_{1x}+c({\cal{R}}^{2}\Phi_{1})^{2}=0,$
(48)
$\displaystyle\Phi_{2t}+\Phi_{2xxx}+a\Phi_{2}({\cal{R}}^{2}\Phi_{2x})+b({\cal{R}}\Phi_{2})({\cal{R}}\Phi_{2x})+c({\cal{R}}^{2}\Phi_{2})(\Phi_{2x})=0,$
(49)

where $a,\ b$ and $c$ are arbitrary constants and $\Phi_{0},\ \Phi_{1}$ and
$\Phi_{2}$ are the ren-fields with degrees, $0,\ 1$ and $2$, respectively.

The special integrable case (45) for $\\{\alpha=3,j=0\\}$ is related to (47)
with $b=0$ up to a re-scaling. (48) with $\\{a=0,\ c=b\\}$ is equivalent to
the integrable case (45) for $\\{\alpha=3,j=1\\}$. Taking $\\{a=b=0\\}$ in
(49) leads to the equivalent special integrable ren-symmetric KdV equation
(45) with $\\{\alpha=3,j=2\\}$.

Incorporating the explicit forms for

$\Phi_{0}=u+\theta\zeta+\theta^{2}\xi,\ \Phi_{1}=\xi+\theta
u+\theta^{2}\zeta,\ \Phi_{2}=\zeta+\theta\xi+\theta^{2}u$

and the consistent commutation relations

$\xi\theta=q\theta\xi,\ \zeta\theta=q^{2}\theta\zeta,\
\xi\zeta=q^{2}\zeta\xi,\ \zeta\zeta_{x}=q\zeta_{x}\zeta,\
\zeta\zeta_{xx}=q\zeta_{xx}\zeta$ (50)

leads to the coupled component forms of (47), (48) and (49)

$\displaystyle u_{t}+u_{xxx}+auu_{x}-q^{2}b\zeta\cdot\xi=0,$
$\displaystyle\zeta_{t}+\zeta_{xxx}+a(\zeta u)_{x}+b\zeta u_{x}+bq\xi^{2}=0,$
$\displaystyle\xi_{t}+\xi_{xxx}+a(u\xi)_{x}+(a-b)\zeta_{x}\cdot\zeta=0,$ (51)
$\displaystyle\xi_{t}+\xi_{xxx}+a\xi u_{x}+bu\xi_{x}+cq\zeta^{2}=0,$
$\displaystyle
u_{t}+u_{xxx}+(a+b)uu_{x}-aq^{2}\zeta_{x}\cdot\xi-(2cq^{2}+b)\xi_{x}\cdot\zeta=0,$
$\displaystyle\zeta_{t}+\zeta_{xxx}+(a-bq-cq^{2})\zeta
u_{x}+(b-aq^{2})u\zeta_{x}-aq\xi_{xx}\cdot\xi+q(c-b)\xi_{x}^{2}=0,$ (52)

and

$\displaystyle\zeta_{t}+\zeta_{xxx}-aq^{2}\zeta
u_{x}+b\xi\cdot\xi_{x}-cq^{2}u\zeta_{x}=0,$ $\displaystyle
u_{t}+u_{xxx}-q(aq+cq-b)uu_{x}+(aq-b)\xi\cdot\zeta_{xx}-(bq^{2}+cq-c)\xi_{x}\cdot\zeta_{x}-aq\xi_{xx}\cdot\zeta=0,$
$\displaystyle\xi_{t}+\xi_{xxx}-(aq^{2}+b)\xi
u_{x}-(b+c)q^{2}u\xi_{x}+a\zeta_{xx}\cdot\zeta+c\zeta_{x}^{2}=0,$ (53)

respectively. $u$ in (51)–(53) is a bosonic field and $\xi$ and $\zeta$ are
ren-fields with degrees $1$ and $2$, respectively.

The known special integrable ren-symmetric KdV systems of (51) and (52) read

$\displaystyle u_{t}+u_{xxx}-6uu_{x}=0,$
$\displaystyle\xi_{t}+\xi_{xxx}-6(u\xi)_{x}-6\zeta_{x}\cdot\zeta=0,$
$\displaystyle\zeta_{t}+\zeta_{xxx}-6(\zeta u)_{x}=0,$ (54)

and

$\displaystyle u_{t}+u_{xxx}-6uu_{x}-6q(1-q)\xi_{x}\cdot\zeta=0,$
$\displaystyle\xi_{t}+\xi_{xxx}-6u\xi_{x}-6q\zeta^{2}=0,$
$\displaystyle\zeta_{t}+\zeta_{xxx}-6(\zeta u)_{x}=0,$ (55)

respectively. The special integrable ren-symmetric KdV system of (53) is
equivalent to that of (52) by the transformation $\zeta_{x}\rightarrow\zeta$.

The choice $\alpha=4,\ q=\mbox{\rm i},\
{\cal{R}}=\partial_{\theta}-\frac{1+q}{2}\theta^{3}\partial_{x}$ leads ren-
symmetric KdV equation (44) straightforwardly to

$\displaystyle\Phi_{0t}+\Phi_{0xxx}+a\Phi_{0}\Phi_{0x}+b({\cal{R}}\Phi_{0})({\cal{R}}^{3}\Phi_{0})+c({\cal{R}}^{2}\Phi_{0})^{2}=0,$
(56)
$\displaystyle\Phi_{1t}+\Phi_{1xxx}+a\Phi_{1}{\cal{R}}\Phi_{1x}+b({\cal{R}}\Phi_{1})\Phi_{1x}+c({\cal{R}}^{2}\Phi_{1})({\cal{R}}^{3}\Phi_{1})=0,$
(57)
$\displaystyle\Phi_{2t}+\Phi_{2xxx}+a\Phi_{2}{\cal{R}}^{2}\Phi_{2x}+b({\cal{R}}\Phi_{2})({\cal{R}}\Phi_{2x})+c({\cal{R}}^{2}\Phi_{2})(\Phi_{2x})+d({\cal{R}}^{3}\Phi_{2})^{2}=0,$
(58)
$\displaystyle\Phi_{3t}+\Phi_{3xxx}+a\Phi_{3}{\cal{R}}^{3}\Phi_{3x}+b({\cal{R}}\Phi_{3})({\cal{R}}^{2}\Phi_{3x})+c({\cal{R}}^{2}\Phi_{3})({\cal{R}}\Phi_{3x})+d\Phi_{3x}{\cal{R}}^{3}\Phi_{3}=0,$
(59)

where $\Phi_{0},\ \Phi_{1},\ \Phi_{2}$ and $\Phi_{3}$ are the ren-fields with
the degrees, $0,\ 1,\ 2$ and $3$, respectively.

##  V Summary and discussions

In retrospect, the usual Grassmann number and the super-symmetric derivative
have been straightforwardly extended to more general forms, the ren-number and
the ren-symmetric derivatives, to be applicable to describe physically
important quasi-particles, anyons. Applying the ren-numbers and ren-symmetric
derivatives to integrable theory, we have extended the super-integrable and
super-symmetric integrable systems to ren-integrable and ren-symmetric
integrable systems.

It is interesting that the ren-integrable KdV system (34) possesses the
completely same form for arbitrary $\alpha$ even for the boson case
($\alpha=\infty$) and fermion case ($\alpha=2$). The only difference is that
the degrees of the ren-fields $\xi_{1}$ and $\xi_{2}$ should be complementary,
say, $\beta$ and $\alpha-\beta$, such that $\xi_{1}\xi_{2}$ becomes a boson.

The ren-integrable system (34) can be further extended to

$\displaystyle
u_{t}=\left[3u^{2}-u_{xx}+\langle\phi|\phi\rangle+12\sum_{\alpha=2}^{\infty}\sum_{\beta_{\alpha}=0}^{\alpha-1}(\langle\xi_{\\{\alpha,\beta_{\alpha}\\}}|\zeta_{\\{\alpha,\alpha-\beta_{\alpha}\\},x}\rangle-\langle\xi_{\\{\alpha,\beta_{\alpha}\\},x}|\zeta_{\\{\alpha,\alpha-\beta_{\alpha}\\}}\rangle)\right]_{x},$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})|\phi\rangle=0,$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})|\xi_{\\{\alpha,\beta_{\alpha}\\}}\rangle=0,$
$\displaystyle(\partial_{t}+4\partial_{x}^{3}-6u\partial_{x}-3u_{x})|\zeta_{\\{\alpha,\alpha-\beta_{\alpha}\\}}\rangle=0,$
(60)

where $|\phi\rangle$ is a boson vector field,
$|\xi_{\\{\alpha,\beta_{\alpha}\\}}\rangle$ is a $\beta_{\alpha}$ order ren-
vector field and $|\zeta_{\\{\alpha,\alpha-\beta_{\alpha}\\}}\rangle$ is an
$\alpha-\beta_{\alpha}$ order ren-vector field. The general ren-integrable KdV
type system (60) is still a Lax integrable model.

Although, indeed the number of papers produced so far on construction of
solutions is incredibly large, it is necessary to develop some novel methods,
one of which may be the so-called bosonization method GA1 , to construct
special solutions of the ren-integrable KdV system (34) (or more generally
(60)) and the ren-symmetric KdV system (44).

Ren-numbers may also be used to find other types of integrable models such as
the dark equations and integrable couplings darkP ; dP . The concept of dark
equations is first introduced by Kupershmidt in Ref. Dark  where many types of
dark KdV systems are given. The modified dark KdV equations are studied in
blueRef. dark1 . The bosonization procedure of the super-symmetric systems
have offered some new types of dark integrable systems GA1 . The bosonization
of ren-symmetric integrable models may yield further dark integrable
equations. In fact, applying the bosonization assumptions

$\displaystyle\xi=\eta p,\ \zeta=\eta^{2}q$ (61)

with the $\\{x,\ t\\}$-independent ren-number $\eta$ and the $\\{x,\
t\\}$-dependent boson fields $p$ and $q$ on the integrable systems (54) and/ot
(55) yields a same standard dark equation system,

$\displaystyle u_{t}+u_{xxx}-6uu_{x}=0,$ $\displaystyle
p_{t}+p_{xxx}-6(up)_{x}=0,$ $\displaystyle q_{t}+q_{xxx}-6(uq)_{x}=0,$ (62)

because of $\eta^{3}=0$. From (62), we know that the ren-integrable systems
(54) and (55) possess a special types of exact solutions with $u$ being an
arbitrary solutions of the usual KdV equation and $\xi$ and $\zeta$ being
given by (61) while $p$ and $q$ are arbitrary symmetries of the usual KdV
equation.

The dark systems can also be considered as some special type of integrable
couplings Ma ; Ma1 . The more about the ren-integrable, ren-symmetric
integrable and dark integrable systems should be further studied later.

###### Acknowledgements.

The work was sponsored by the National Natural Science Foundations of China
(Nos.12235007, 11975131). The author is indebt to thank Profs. Q. P. Liu, B.
F. Feng, X. B. Hu, R. X. Yao, M. Jia and Drs. K. Tian, X. Z. Hao and D. D.
Zhang for their helpful discussions.

## References

## References

  * (1) E. S. Abers and B. W. Lee, Gauge theories, Phys. Rept. 9 (1973) 1.
  * (2) Y. Q. Li, J. C. Chen, Y. Chen and S. Y. Lou, Darboux transformations via lie pointsymmetries: KdV equation, Chin. Phys. Lett. 31 (2014) 010201.
  * (3) P. A. Clarkson and M. D. Kruskal, New similarity reductions of the Boussinesg equation, J. Math. Phys. 30 (1989) 2201.
  * (4) P. A. Clarkson and E. L. Mansfield, Algorithms for the nonclassical method of symmetry reductions, SIAM. Appl. Math. 54 (1994) 1693.
  * (5) S. Y. Lou, Similarity solutions of the Kadomtseu-Petviashvili equation, J. Phys. A 23 (1990) L649.
  * (6) S. Y. Lou, H. Y. Ruan, D. F. Chen and W. Z.Chen, Similarity reductions of the KP equation by a direct method, J. Phys. A 24 (1991) 1455.
  * (7) C. W. Cao and X. G. Geng, Neumann and Bargmann systems associated with the coupled Kdsoliton hierarchy, J. Phys. A 23 (1991) 4117.
  * (8) Y. Cheng and Y. S Li, The constraint of the Kadomtseu-Petviashvili equation and its special solutions, Phys. Lett. A 157 (1991) 22.
  * (9) B. Konopelchenko, J. Sidorenko and W. Strampp, (1+1)-dimensional integrable systems as symmetry constraints of (2+1)-dimensional systems, Phys. Lett. A 157 (1991) 1723.
  * (10) S. Y. Lou and X. B. Hu, Infinitely many Lax pairs and symmetry constraints of the KP equation, J. Math. Phys. 38 (1997) 6401.
  * (11) S. Y. Lou and R. X. Yao, Primary branch solutions of frst order autonomous scalar partialdifferential equations via Lie symmetry approach, J. Nonlin. Math. Phys. 24 (2017) 379.
  * (12) D. Z. Freedman and J. H. Schwarz, N=4 supergravity theory with local SU(2)$\times$SU(2) invariance, Nucl. Phys. B 137 (1978) 333.
  * (13) C. Angelantonj, S. Ferrara and M. Trigiante, New D = 4 gauged supergravities from N=4 orientifolds with fluxes, JHEP 10 (2003) 015.
  * (14) F. Cooper, A. Khare and U. Sukhatme, Supersymmetry and quantum mechanics, Phys. Rep. 251 (1995) 267.
  * (15) J. Polchinski, String theory, Vol. 2: Superstring theory and beyond, Cambridge University Press, 2007.
  * (16) Yu. I. Manin and A. O. Radul, A Supersymmetric extension of the Kadomtseu-Petuiashvili Hierarchy, Comm. Math. Phys. 98(1985) 65.
  * (17) P. Mathieu, Supersymmetric extension of the Korteweg-de Vries equation, J. Math. Phys. 29 (1988) 2499.
  * (18) K. Tian and Q. Liu, A supersymmetric Sawada-Kotera equation, Phys. Lett. A 373 (2009) 1807.
  * (19) P. Labelle and P. Mathieu, N = 2 superconformal algebra and integrable O(2) fermionic extensions of the Korteweg-de Vries equation, Phys. Lett. B 215 (1988) 718
  * (20) X. N. Gao and S. Y. Lou, Bosonization of supersymmetric KdV equation, Phys. Lett. B 707 (2012) 209.
  * (21) M. X. Zhang, Q. P. Liu, Y. L. Shen and K. Wu, Bilinear approach to N=2 supersymmetric KdV equations, Sci. China Ser. A, 52 (2009) 1973.
  * (22) D. Krakowski and A. Regev, The polynomial identities of the Grassmann algebra, Trans. Amer. Math Soc. 181 (1973) 429.
  * (23) X. N. Gao, S. Y. Lou and X. Y. Tang, Bosonization, singularity analysis, nonlocal symmetry reductions and exact solutions of supersymmetric KdV equation, JHEP, 5 (2013) 029.
  * (24) F. Wilczek, Quantum mechanics of fractional-spin particles. Phys. Rev. Lett. 49 (1982) 957.
  * (25) F. Wilczek, Fractional statistics and anyon superconductivity (World Scientific, 1990).
  * (26) C. Nayak, S. H. Simon, A. Stern, M. Freedman and S. D. Sarma, Non-Abelian anyons andtopological quantum computation, Rev. Mod. Phys. 80 (2008) 1083.
  * (27) Goodle Quantum Collaborators, Non-Abelian braiding of graph vertices in a superconducting processor, Nature, (2023) https:/doi.org/10.1038/s41586-023-05954-4.
  * (28) G. Moore and N. Read, Nonabelions in the fractional quantum Hall effect, Nucl. Phys. B 360 (1991) 362.
  * (29) D. A. Ivanov, Non-Abelian statistics of half-quantum vortices in P-wave superconductors, Phys. Rev. Lett. 86 (2001) 268.
  * (30) R. M. Lutchyn, J. D. Sau and S. Das Sarma, Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105 (2010) 077001.
  * (31) B. A. Kupershmidt, A super Korteweg-de Vries equation: an integrable system, Phys. Lett. 102A (1984) 213.
  * (32) B. A. Kupershmidt, Super-integrable systems, Proc. Nat. Acad. Sci. USA, 81 (1984), 6562.
  * (33) B. A. Kupershmidt, Bosons and fermions interacting integrably with the Korteweg-de Vries field, J. Phys. A 17 (1984), L869.
  * (34) B. A. Kupershmidt, A review of superintegrable systems, Lectures in Appl. Math. 23 Part 1 (1986) 83.
  * (35) S. Majid and M. J. Rodr aguez-Plaza, Random walk and the heat equation on superspace and anyspace, J. Math. Phys. 35 (1994) 3753.
  * (36) F. Berezin, Introduction to superanalysis (Reidel, Dordrecht, 1987), 1st ed.
  * (37) M. Kieburg, H. Kohler, and T. Guhr, Integration of Grassmann variables over invariant functions on flat superspaces, J. Math. Phys. 50 (2009) 013528.
  * (38) V. K. Mel? nikov, On equations for wave interactions, Lett. Math. Phys. 7 (1983) 129.
  * (39) V. K. Mel? nikov, A direct method for deriving a multi-soliton solution for the problem of interaction of waves on the x, y plane, Commun. Math. Phys. 112 (1987) 639.
  * (40) V. K. Mel? nikov, Interaction of solitary waves in the system described by the Kadomtsev-Petviashvili equation with a self-consistent source, Commun. Math. Phys. 126 (1989) 201.
  * (41) Y. B. Zeng, W. X. Ma and R. L. Lin, Integration of the soliton hierarchy with self-consistent sources, J. Math. Phys. 41 (2000) 5453.
  * (42) X. B. Hu and H. Y. Wang, Construction of dKP and BKP equations with self-consistent sources, Inverse Problems 22 (2006) 1903.
  * (43) X. B. Hu and H. Y. Wang, New type of Kadomtsev-Petviashvili equation with self-consistent sources and its bilinear Bäcklund transformation, Inverse Problems 23 (2007) 1433.
  * (44) H. Y. Wang and X. B. Hu, A hybrid type of soliton equations with self-consistent sources: KP and Toda cases, J. Math. Anal. Appl. 376 (2011) 393.
  * (45) D. J. Zhang, The soliton solutions of some soliton equations with self-consistent sources, Chaos, Soliton & Fractals, 18 (2003) 31.
  * (46) H. J. Tian and D. J. Zhang, Cauchy matrix approach to integrable equations with self-consistent sources and the Yajima-Oikawa system, Appl. Math. Lett. 103 (2020) 106165.
  * (47) P. Mathieu, The Painlevé property for fermionic extensions of the Korteweg-de Vries equation, J. Math. Phys. 29 (1988) 2499.
  * (48) X. N. Gao, X. D. Yang and S. Y. Lou, Exact solutions of supersymmetric KdV-a system via bosonization approach, Commun. Theor. Phys. 58 (2012) 617.
  * (49) S. Y. Lou, Dark parameterizations, equivalent partner fields and integrable systems, Commun. Theor. Phys. 55 (2011) 743.
  * (50) L. Peng, X. D. Yang and S. Y. Lou, Integrable KP coupling and its exact solution, Commun. Theor. Phys. 58 (1012) 1.
  * (51) B. A. Kupershmidt, Dark equations J. Nonlinear Math. Phys. 8 (2001) 363.
  * (52) N. Xiong, S. Y. Lou, B. Li and Y. Chen, Classification of dark modified KdV equation, Commun. Theor. Phys. 68 (2017) 13.
  * (53) W. X. Ma and B. Fuchssteiner, The bi-Hamiltonian structure of the perturbation equations of the KdV hierarchy, Phys. Lett. A 213 (1996) 49.
  * (54) W. X. Ma, Integrable couplings of vector AKNS soliton equations, J. Math. Phys. 46 (2005) 033507.


# Bayesian autoencoders with uncertainty quantification: Towards trustworthy
anomaly detection

Bang Xiang Yong<EMAIL_ADDRESS>Alexandra Brintrup<EMAIL_ADDRESS>Institute
for Manufacturing, University of Cambridge, United Kingdom

###### Abstract

Despite numerous studies of deep autoencoders (AEs) for unsupervised anomaly
detection, AEs still lack a way to express uncertainty in their predictions,
crucial for ensuring safe and trustworthy machine learning systems in high-
stake applications. Therefore, in this work, the formulation of Bayesian
autoencoders (BAEs) is adopted to quantify the total anomaly uncertainty,
comprising epistemic and aleatoric uncertainties. To evaluate the quality of
uncertainty, we consider the task of classifying anomalies with the additional
option of rejecting predictions of high uncertainty. In addition, we use the
accuracy-rejection curve and propose the weighted average accuracy as a
performance metric. Our experiments demonstrate the effectiveness of the BAE
and total anomaly uncertainty on a set of benchmark datasets and two real
datasets for manufacturing: one for condition monitoring, the other for
quality inspection.

###### keywords:

Bayesian autoencoders, anomaly detection, uncertainty, trustworthy machine
learning

††journal: Expert Systems with Applications

##  1 Introduction

The use of unsupervised neural networks (NNs) such as deep autoencoders (AEs)
has consistently achieved improvements over traditional machine learning (ML)
models for high-dimensional anomaly detection [1, 2, 3]. Nonetheless, one
missing ingredient is the measurement of predictive uncertainty. That is, in
addition to reporting the prediction of an anomaly, can the AE tell how
uncertain it is about the prediction?

The model’s prediction and its uncertainty are two different outputs: the
prediction classifies a data point as inlier or anomaly by relying on an
anomaly score, which typically measures its similarity or distance from a
reference distribution of normality; the higher the anomaly score, the higher
the chance of an abnormal observation [4]. By contrast, the anomaly
uncertainty measures the trustworthiness of the classification; the higher the
uncertainty, the less reliable is the prediction, implying a greater chance of
a predictive error. (See Fig. 1 for an illustration).

Figure 1: An uncertainty-aware anomaly detector quantifies predictive
uncertainty in addition to predicting inliers and anomalies.

One popular method to quantify uncertainty for NNs is by adopting the Bayesian
framework, resulting in the Bayesian autoencoder (BAE) formulation. Two types
of uncertainty are of interest: epistemic uncertainty captures the uncertainty
in the BAE parameters, while aleatoric uncertainty reflects the inherent
randomness in the data [5]. However, recent works on BAEs [6, 7, 8, 9, 10]
have only investigated the use of epistemic and aleatoric uncertainties as
anomaly scores, leaving the notion of anomaly uncertainty understudied.
Consequently, we are still incapable of knowing when to (dis)trust the BAEs’
predictions.

Knowing the uncertainty of each prediction is important for ensuring
reliability, especially in high-stake applications, as we can anticipate
whether there is an error with the prediction [11]. If we filter away cases of
high uncertainty, we are left with trustworthy predictions of higher accuracy.
In addition, these uncertain cases can be referred for further diagnosis,
either by humans or equipment, thereby increasing the overall predictive
performance. Therefore, in this work, we contribute in the following ways:

  1. 1.

Formulation of BAE to quantify predictive uncertainty of anomalies. We do so
by converting negative log-likelihood (NLL) estimates from the BAE into the
total anomaly uncertainty which is decomposable into epistemic and aleatoric
uncertainties.

  2. 2.

Evaluation of uncertainty in anomaly detection. We consider the task of
classifying anomalies with the option of rejecting uncertain predictions.
While these evaluations have been explored in supervised learning [12, 13,
14], we contribute by extending to unsupervised anomaly detection.

  3. 3.

Validation of proposed BAE with real use cases that demonstrate its
effectiveness and value. Our experiments demonstrate the added benefits of
quantifying anomaly uncertainty on a set of benchmark datasets, as well as two
real datasets for manufacturing applications: one for condition monitoring and
the other for quality inspection.

##  2 Background and related works

Anomaly detection***Also known as novelty detection, outlier detection, and
more recently, out-of-distribution detection. is a binary classification task
of identifying data points that differ significantly from a reference
distribution of inliers [4]. Anomaly detection has many real use cases in
high-stake applications, e.g. cybersecurity [15], financial surveillance [16],
health and medical diagnosis [17], manufacturing [18]. Many traditional ML
models have been adapted to perform anomaly detection, e.g. isolation forest,
one-class support vector machine, K-means clustering [19]. A drawback of
traditional ML models is their low complexity and flexibility, limiting their
scalability from handling high-dimensional data prevalent in modern databases
[20].

Overcoming the limitations of traditional ML models, deep NNs are expressive
and scalable models. Many works have studied their use for anomaly detection,
consistently achieving state-of-the-art results [3, 17, 18, 15, 16]. A key
property of NNs is their flexibility of modelling; each component (e.g.
activation function, layers, optimiser) can be tuned and customised to adapt
to the data at hand. Furthermore, specialised hardware have expedited their
adoption and scalability by parallelising and optimising computations [21].

An emerging challenge with using NNs for anomaly detection is the
quantification of uncertainty [20]. Uncertainty is a core component for
promoting algorithmic transparency and hence, the advancement of trustworthy
ML [22]. While most works have focused on anomaly scores, less is known about
anomaly uncertainty. To this end, Perini et al. [23] introduced the example-
wise confidence method for quantifying anomaly uncertainty, however, their
methods have only experimented with traditional ML models, leaving their
applicability to deep AEs unclear.

The development of Bayesian neural networks (BNNs) [24] for uncertainty
quantification has led to the formulation of BAEs, a probabilistic version of
the AE. Nonetheless, extant works have used uncertainty as an anomaly score,
which is different from the notion of anomaly uncertainty. We list extant
works as follows: Yong et al. [25] quantified epistemic and aleatoric
uncertainties for detecting sensor drifts. Daxberger & Lobato [8] developed
the Bayesian variational autoencoder (BVAE) to quantify uncertainty as a
measure of deviation from the training distribution. Baur et al. [10]
developed the BAE with Monte Carlo dropout for anomaly detection in medical
application. Tran et al. [9] observed that BAE-generated data have higher
uncertainty than the reconstructed data. Chandra et al. [26] studied the BAE
with Markov Chain Monte Carlo (MCMC) sampling for dimensionality reduction and
classification tasks.

Table 1: Comparisons between related works and ours on ML models for anomaly
uncertainty.

Work | ML model | Anomaly uncertainty
---|---|---
Ours | BAE | ✓
Perini et al. [23] | Traditional | ✓
Yong et al. [25] | BAE | ✗
Daxberger & Lobato [8] | BVAE | ✗
Baur et al. [10] | BAE | ✗
Tran et al. [9] | BAE | ✗
Chandra et al. [26] | BAE | ✗
Kriegel et al. [27] | Traditional | ✗
Kwon et al. [28] | Supervised BNN | ✗

Our work stands out as the first to provide a comprehensive study on the BAE’s
predictive uncertainty of an anomaly. We extend the work of Kwon et al. [28],
which decomposed the classification uncertainty of supervised BNNs into
epistemic and aleatoric uncertainties, to consider unsupervised anomaly
detection with BAEs. Further, our work is inspired by Kriegel et al. [27], who
proposed converting an ensemble of anomaly scores into anomaly probabilities,
and combining them via averaging. To improve the combination quality, they
further proposed a customised scaling. We extend their study by (1) including
BAEs, (2) quantifying anomaly uncertainty while maintaining coherence within a
probabilistic framework, and (3) finding that their customised scaling
improved the quality of uncertainty. Comparisons with related works are
summarised in Table 1.

##  3 Methods

This section details the probabilistic formulation of BAEs for quantifying
anomaly uncertainty, followed by ways to evaluate them.

###  3.1 Bayesian autoencoder

Suppose we have a set of data
$X^{train}={\\{\textbf{x}_{1},\textbf{x}_{2},\textbf{x}_{3},...\textbf{x}_{N}\\}}$,
$\textbf{x}_{i}\in\rm I\\!R^{D}$. An AE is an NN parameterised by $\theta$,
and consists of two parts: an encoder $f_{\text{encoder}}$, for mapping input
data x to a latent embedding, and a decoder $f_{\text{decoder}}$ for mapping
the latent embedding to a reconstructed signal of the input $\hat{\textbf{x}}$
(i.e.
$\hat{\textbf{x}}={f_{\theta}}(\textbf{x})=f_{\text{decoder}}(f_{\text{encoder}}(\textbf{x}))$)
[29].

Bayes’ rule can be applied to the parameters of the AE to create a BAE,

$p(\theta|X^{train})=\frac{p(X^{train}|\theta)\ p(\theta)}{p(X^{train})}\\\ ,$
(1)

where $p(X^{train}|\theta)$ is the likelihood and $p(\theta)$ is the prior
distribution of the AE parameters. The log-likelihood for a diagonal Gaussian
distribution is,

$\log{p(\textbf{x}|\theta)}=-(\frac{1}{D}\sum^{D}_{i=1}{\frac{1}{2\sigma_{i}^{2}}}(x_{i}-\hat{x_{i}})^{2}+\frac{1}{2}\log{\sigma_{i}^{2}})$
(2)

where $\sigma_{i}^{2}$ is the variance of the Gaussian distribution. For
simplicity, we use an isotropic Gaussian likelihood with $\sigma_{i}^{2}=1$ in
this study, since the mean-squared error (MSE) function is proportional to the
NLL. For the prior distribution, we use an isotropic Gaussian prior which
effectively leads to $L_{2}$ regularisation.

Since Equation 1 is analytically intractable for a deep NN, various
approximate methods have been developed such as Stochastic Gradient Markov
Chain Monte Carlo (SGHMC) [30], Monte Carlo Dropout (MCD) [31], Bayes by
Backprop (BBB) [32], and anchored ensembling [33] to sample from the posterior
distribution. In contrast, a deterministic AE has its parameters estimated
using maximum likelihood estimation (MLE) or maximum a posteriori (MAP) when
regularisation is introduced. The variational autoencoder (VAE) [34] and BAE
are AEs formulated differently within a probabilistic framework: in the VAE,
only the latent embedding is stochastic while the $f_{\text{encoder}}$ and
$f_{\text{decoder}}$ are deterministic and the model is trained using
variational inference, while, on the other hand, the BAE (similar to BNN) has
distributions over all parameters of $f_{\text{encoder}}$ and
$f_{\text{decoder}}$. See Appendix A for descriptions of the posterior
sampling methods used in our study.

In short, the training phase of BAE entails using one of the sampling methods
to obtain a set of approximate posterior samples
$\\{\hat{\theta}_{m}\\}^{M}_{m=1}$. Then, during the prediction phase, we use
the posterior samples to compute $M$ estimates of the NLL. For brevity, we
denote the posterior NLL score
$\scalebox{0.75}[1.0]{$-$}\log{p(\textbf{x}|\hat{\theta}_{m})}$, as
$\ell{}_{\hat{\theta}_{m}}({\textbf{x})}$. Next, we look into using the NLL
scores to quantify the anomaly uncertainty.

###  3.2 Quantifying anomaly uncertainties

Overview. The proposed method for quantifying anomaly uncertainty is
illustrated in Fig. 2. Although we use the NLL as an anomaly score, it is not
a true probability of an anomalous outcome. Hence, we first convert the NLL
scores into anomaly probabilities via the cumulative distribution function
(CDF). Next, using the anomaly probabilities, we compute the epistemic and
aleatoric uncertainties. Capturing both uncertainties is crucial for forming a
holistic quantification of predictive uncertainty; losing one of them leads to
incomplete quantification and hence a lower quality of uncertainty estimate.
Summing them forms the total anomaly uncertainty. Now, we formally describe
our method for obtaining the total anomaly uncertainty with the BAE.

Figure 2: Workflow for quantifying the probability of an anomalous event and
the predictive uncertainty using an ensemble of $M$-samples from the BAE
posterior.

(i) Distribution of anomaly scores.

Let $\phi_{m}$ be the vector of parameters of the distribution of NLL scores
from the $m$-th BAE posterior sample,

$\ell{}_{\hat{\theta}_{m}}({\textbf{x})}\thicksim\text{Q}(\phi_{m})$ (3)

where Q is an arbitrary distribution. The optimal choice of distribution Q is
dataset dependent, and in this study, we have experimented with setting Q as
either a Gaussian, exponential or uniform distribution. For simplicity, we
estimate $\phi_{m}$ by applying maximum likelihood estimation on the anomaly
scores computed over the training data, i.e.
$\hat{\phi}_{m}=\operatorname*{arg\,max}_{\phi_{m}}\prod_{i=1}^{N}p(\ell{}_{\hat{\theta}_{m}}({\textbf{x}_{i}^{train})}\,|\,{\phi_{m}})$.

(ii) Quantification of anomaly probability. We define $y\in\\{0,1\\}$ as the
outcome of an anomaly for a new data point $\textbf{x}^{*}$, which has a
Bernoulli distribution,

$y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}\thicksim\text{Ber}(p_{m}^{anomaly}(\textbf{x}^{*}))$
(4)

where the anomaly probability $p_{m}^{anomaly}(\textbf{x}^{*})$ is modelled
using the CDF,

$p_{m}^{anomaly}(\textbf{x}^{*})=\text{CDF}_{\hat{\phi}_{m}}(\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})})$
(5)

Using the CDF to quantify anomaly probability is not uncommon; applying a min-
max scaler [35] on the anomaly scores to rescale into $[0,1]$ is the same as
fitting to a uniform distribution and using its CDF.

Alternatively, the empirical cumulative distribution function (ECDF) can be
used to estimate $p_{m}^{anomaly}(\textbf{x}^{*})$. The ECDF essentially
captures the fraction of anomaly scores of the training data with lower values
than the observed test anomaly score. Moreover, the ECDF converges to the true
CDF as the number of examples increases according to Glivenko-Cantelli theorem
[36]. As such, the ECDF benefits from an increasing number of training data,
and we can avoid having to specify a distribution Q for the anomaly scores by
leveraging the ECDF.

$\text{ECDF}(\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})})=\frac{1}{N}\sum_{i=1}^{N}1_{\ell{}_{\hat{\theta}_{m}}({\textbf{x}_{i}^{train})}\,\leq\,{\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})}}}$
(6)

In addition to using the CDF or ECDF when converting anomaly score to
probability, we may choose to apply a customised scaling method proposed by
Kriegel et al. [27],

$p_{m}^{Z_{anomaly}}(\textbf{x}^{*})=\text{max}\bigg{\\{}0,\frac{\text{CDF}_{\hat{\phi}_{m}}(\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})})-\text{CDF}_{\hat{\phi}_{m}}(\operatorname*{\mathbb{E}}[\ell{}_{\hat{\theta}_{m}}({{\textbf{x}^{train}})}])}{1-\text{CDF}_{\hat{\phi}_{m}}(\operatorname*{\mathbb{E}}[\ell{}_{\hat{\theta}_{m}}({{\textbf{x}^{train}})}])}\bigg{\\}}$
(7)

reflecting the beliefs that the anomaly probability is zero when the observed
anomaly score is less than or equal to the average anomaly score of the
training set (i.e. $p_{m}^{Z_{anomaly}}(\textbf{x}^{*})=0$ for
$\text{CDF}_{\hat{\phi}_{m}}(\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})})\leq{\text{CDF}_{\hat{\phi}_{m}}(\operatorname*{\mathbb{E}}[\ell{}_{\hat{\theta}_{m}}({{\textbf{x}^{train}})}])}$),
and that the anomaly probability increases as the anomaly score gets higher
than the average anomaly score.

By definition, the mean and variance of the Bernoulli distribution in Eq. 4
are:

$\displaystyle\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]=p_{m}^{anomaly}(\textbf{x}^{*})$
(8)
$\displaystyle\text{Var}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]=p_{m}^{anomaly}(\textbf{x}^{*})(1-p_{m}^{anomaly}(\textbf{x}^{*}))$
(9)

See Fig. 3 for an illustrative example of Eq. 8 and Eq. 9.

Figure 3: Conversion of $\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})}$ from a
single sample of BAE into
$\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]$
(upper panel) and
$\text{Var}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]$ (lower
panel). Distribution Q is chosen to be either the Gaussian or uniform
distribution with an option of applying the customised scaling method. The
customised scaling shifts the graph to the right by treating
$\operatorname*{\mathbb{E}}[\ell{}_{\hat{\theta}_{m}}({{\textbf{x}^{train}})}]$
as a reference.

Next, to predict an anomalous event, we summarise the anomaly probabilities
from all BAE posterior samples via expectation,

$\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train}]\approx{}\operatorname*{\mathbb{E}}{}_{{\theta},{\phi}}[\,\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]\,]=\frac{1}{M}\sum_{m=1}^{M}p_{m}^{anomaly}(\textbf{x}^{*})$
(10)

A reasonable threshold can be set for the mean anomaly probability at any
point where $\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train}]\geq 0.5$ to
get a hard prediction of an anomalous event. Moreover, as a calibrated anomaly
score, its value ranges in [0,1] and facilitates better interpretation of an
anomalous event than the raw anomaly score, which takes on any real number
[27].

(iii) Decomposition of anomaly uncertainty. Based on the law of total variance
[37], we decompose the total anomaly uncertainty, $U^{total}$ into its
epistemic and aleatoric components,

$\begin{split}\text{Var}[y|\textbf{x*},X^{train}]\approx{\underbrace{\operatorname*{\mathbb{E}}{}_{{\theta},{\phi}}[\text{Var}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]]}_{\text{$U^{aleatoric}$}}}\,+\\\
\underbrace{\text{Var}_{{\theta},{\phi}}[\operatorname*{\mathbb{E}}{}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]]}_{\text{$U^{epistemic}$}}\end{split}$
(11)

Substituting with Eq. 8 and Eq. 9, the $U^{aleatoric}$ and $U^{epistemic}$ can
be computed as

$\displaystyle\operatorname*{\mathbb{E}}{}_{{\theta},{\phi}}[\text{Var}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]]=\operatorname*{\mathbb{E}}{}_{m}[\\{p_{m}^{anomaly}(\textbf{x}^{*})(1-p_{m}^{anomaly}(\textbf{x}^{*}))\\}_{m=1}^{M}]$
(12)
$\displaystyle\text{Var}_{{\theta},{\phi}}[\operatorname*{\mathbb{E}}{}[y|\textbf{x*},X^{train},\hat{\theta}_{m},\hat{\phi}_{m}]]=\text{Var}_{m}[\\{p_{m}^{anomaly}(\textbf{x}^{*})\\}_{m=1}^{M}]$
(13)

Our method propagates the model uncertainty from $M$ samples of the BAE
posterior to yield a holistic quantification of anomaly uncertainty. In
contrast, when we have a deterministic AE (i.e. $M=1$), the $U^{epistemic}$
term is zero and the remaining component is the $U^{aleatoric}$, signifying
the inability of the deterministic AE to capture epistemic uncertainty.

Note on interpretability. Reading the $U^{total}$ can be unintuitive since the
variance of a Bernoulli distribution ranges in [0,0.25]. Therefore, to enhance
interpretability, we multiply $U^{total}$ $\times 4$ to range between [0,1] in
this study.

###  3.3 Alternative ways to quantify anomaly uncertainties

We investigate two other methods for quantifying anomaly uncertainty. Firstly,
by applying the method proposed by Perini et al. [23] on the BAE under the
assumption that the training set is not polluted with anomalies, we arrive at
the following uncertainty estimate,

$U^{exceed}=\begin{cases}(\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train}])^{N}&,\,\operatorname*{\mathbb{E}}{}_{\theta}[\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})}]<\text{max}\\{\operatorname*{\mathbb{E}}{}_{\theta}[\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{train}_{n})}]\\}_{n=1}^{N}\\\
1-(\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train}])^{N}&,\,\text{otherwise}\end{cases}$
(14)

where $N$ is the number of training examples. Like our method, the
$U^{exceed}$ relies on quantifying the $p_{m}^{anomaly}(\textbf{x}^{*})$. The
difference, however, is in the conversion from anomaly probability to anomaly
uncertainty, which we posit a problem with the $U^{exceed}$ is in the regime
of high $N$, as this leads to extreme values of uncertainty estimate. In
effect, it is prone to over- and underconfidence when the training examples
are abundant, limiting its scalability to larger datasets. In addition, the
$U^{exceed}$ does not account for epistemic uncertainty, leading to a poorer
uncertainty estimate.

The second alternative way to measure anomaly uncertainty is by taking the
variance across the $M$ samples of NLL scores from the BAE or VAE,

$\text{Var}_{\theta}[\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})}]=\frac{1}{M}\sum_{m=1}^{M}(\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})}-\operatorname*{\mathbb{E}}{}_{\theta}[\ell{}_{\hat{\theta}_{m}}({\textbf{x}^{*})}])^{2}$
(15)

We shall refer to this method as $\mathrm{Var_{\theta}[NLL]}$. Note that
$\mathrm{Var_{\theta}[NLL]}$ is not applicable with the deterministic AE since
the number of sampled predictions is 1 ($\mathrm{Var_{\theta}[NLL]}$ =0). See
Fig. 4 for an illustration of the various methods for quantifying anomaly
uncertainty.

Figure 4: Visualisation of (a)
$\operatorname*{\mathbb{E}}[y|\textbf{x*},X^{train}]$ as anomaly scores, and
measures of uncertainties using (b) $U^{epistemic}$, (c) $U^{aleatoric}$, (d)
$U^{total}$, (e) $U^{exceed}$, and (f) $\mathrm{Var_{\theta}[NLL]}$ using a
BAE-Ensemble on a toy dataset with 2 features. Green dots are training
examples while darker regions indicate higher values. $U^{epistemic}$ and
$U^{aleatoric}$ complement each other, forming reasonable regions of
uncertainty with the $U^{total}$. By contrast, $U^{exceed}$ is overconfident
in regions surrounding the training examples.

###  3.4 Evaluating predictive uncertainty in anomaly detection

We evaluate the quality of uncertainty in the task of classifying anomalies
with a reject option. Specifically, in addition to classifying a data point as
either inlier or anomaly, the BAE has the option of withholding its prediction
by not assigning any label. The rationale is that predictions with high
uncertainty are likely to be erroneous; hence rather than accepting all
predictions in a high-stake application, it is safer to accept only
predictions of low uncertainty, which are more accurate [11].

In this setup, we use the accuracy-rejection curve (ARC) proposed by Nadeem et
al. [12], which plots the accuracy of the retained predictions against the
rejection rate evaluated at multiple anomaly uncertainty thresholds. In
effect, the ARC illustrates the trade-off between accuracy and rejection rate.
As a measure of accuracy, while any metrics for binary classification can be
used, we prefer the geometric mean of the sensitivity and specificity ($GSS$)
for evaluating classifiers on an imbalanced dataset [38]. Then, to summarise
the model’s performance, we compute the weighted average accuracy

$W^{GSS}=\frac{\sum_{i}^{T}(100-r_{i})GSS_{i}}{\sum_{i}^{T}(100-r_{i})}$ (16)

where $T$ is the number of uncertainty threshold settings and $r$ is the
rejection rate in percentage. For the anomaly uncertainty to be skillful, the
$W^{GSS}$ should be higher than the baseline $GSS$, $Base^{GSS}$, where $r=0$,
(i.e. not using the anomaly uncertainty). Hence, an important metric is the
difference between the $W^{GSS}$ and the baseline,

$Gain^{GSS}=W^{GSS}-Base^{GSS}$ (17)

which measures the average improvement in accuracy when uncertainty is
included. It is possible for the $Gain^{GSS}$ to be negative in value,
indicating a drop in performance due to the poor utility of uncertainty in
rejecting erroneous predictions.

##  4 Experiment setup

The full code for reproducing the results will be released upon acceptance of
this paper.

###  4.1 Datasets

We use a collection of benchmark datasets, Outlier Detection Datasets (ODDS)
[39], which are common in extant studies [40, 41]. In addition, we use two
publicly available datasets for industrial applications: ZeMA [42] for
condition monitoring and STRATH [43] for quality inspection. In the ZeMA
dataset, there are 4 tasks; each task uses sensor measurements as inputs for
detecting anomalies in hydraulic subsystems, namely, (i) the cooler, (ii) the
valve, (iii) the pump, and (iv) the accumulator. While, in the STRATH dataset,
the tasks are to detect defective parts manufactured from a radial forging
facility. Different sensors are used for each task: (i) a position sensor
(L-ACTpos), (ii) a speed sensor (A-ACTspd), (iii) a servo (Feedback-SPA), and
(iv) all sensors from tasks (i-iii). See [44, 45] for detailed descriptions of
the facility setup and data collection of ZeMA and STRATH datasets.

For all datasets, we split the inliers into train-test sets with a ratio of
70:30 and include all anomalies in the test set. We use the min-max scaler
[35] to transform the input features into [0,1]. The resulting dimensions of
the datasets is tabulated in Table 2. Additional details of preprocessing are
included in Appendix B.

Table 2: Number of examples in train and test sets, and number of features for
each task in ODDS, ZeMA and STRATH datasets.

Datasets & tasks | Train (inliers) | Test (inliers) | Test (anomalies) | Features
---|---|---|---|---
ODDS |  |  |  |
(i) Cardio | 1324 | 331 | 176 | 21
(ii) Lympho | 113 | 29 | 6 | 18
(iii) Optdigits | 4052 | 1014 | 150 | 64
(iv) Pendigits | 5371 | 1343 | 156 | 16
(v) Thyroid | 2943 | 736 | 93 | 6
(vi) Ionosphere | 180 | 45 | 126 | 33
(vii) Pima | 400 | 100 | 268 | 8
(viii) Vowels | 1124 | 282 | 50 | 12
ZeMA |  |  |  |
(i) Cooler | 518 | 223 | 242 | 60
(ii) Valve | 787 | 338 | 30 | 60
(iii) Pump | 854 | 367 | 22 | 60
(iv) Accumulator | 419 | 180 | 272 | 60
STRATH |  |  |  |
(i) L-ACTpos | 51 | 23 | 7 | 562
(ii) A-ACTspd | 51 | 23 | 7 | 562
(iii) Feedback-SPA | 51 | 23 | 7 | 562
(iv) All sensors | 51 | 23 | 7 | 1686

###  4.2 Model hyperparameters

We fit the deterministic AE, VAE, and BAE models to the train set using the
Adam optimiser [46] for 100 epochs with a fixed learning rate of 0.001 for
STRATH. For ODDS and ZeMA, we do not use a fixed learning rate but instead
employ an automatic learning rate finder [47]. The number of posterior samples
is set to $M=100$ for the VAE, BAE-MCD and BAE-BBB, while $M=10$ for the BAE-
Ensemble. The prior scaling term (weight decay) for all models is
$1\times{10^{-10}}$. We specify the architecture of the encoder in Table 3(b).
U-Net skip connections [48] are implemented in all models except for those in
STRATH tasks (i) and (ii).

Table 3: Encoder architecture for ODDS, ZeMA and STRATH datasets. We set the
latent dimensions to be half of the flattened input dimensions. The decoder is
a reflection of the encoder, in which the Conv1D layers are replaced by
Conv1D-Transpose layers. We apply layer normalisation after each Conv1D and
linear layer for ODDS and STRATH. The leaky ReLu [49] is used as the
activation function with a slope of 0.01 while the sigmoid function is used at
the decoder’s final layer.

(a) ODDS

Layer | Output channels/nodes | Kernel | Strides
---|---|---|---
Linear |  Input dimensions $\times 4$ | - | -
Linear |  Input dimensions $\times 4$ | - | -
Linear | Latent dimensions | - | -

(b) ZeMA and STRATH

Layer | Output channels/nodes | Kernel | Strides
---|---|---|---
Conv1D | 10 | 8 | 2
Conv1D | 20 | 2 | 2
Reshape | - | - | -
Linear | 1000 | - | -
Linear | Latent dimensions | - | -

We repeat the experiment 10 times with different random seeds. For each run,
we vary the choice of anomaly score distribution Q to be either the Gaussian,
exponential, or uniform distribution. Due to the large combinations of
methods, we report results where Q is implicitly chosen (unless stated) to
maximise the $W^{GSS}$ criterion.

##  5 Results and discussion

Incorrect predictions have higher uncertainty. From Fig. 5, we observe that
inaccurate predictions are usually accompanied by higher $U^{total}$ scores,
demonstrating the use of $U^{total}$ to anticipate both Type I and Type II
predictive errors [50], in the absence of true labels. Nonetheless, this is
not always the case; for example, in Fig. 5(a)(iii) the VAE assigns similar or
higher $U^{total}$ scores to correct predictions. Likewise, the deterministic
AE does the same for ZeMA task (iv). In such cases, the uncertainty quality is
poor, highlighting the ongoing research challenge of obtaining reliable
uncertainty estimates. Hence, a careful assessment of uncertainty quality is
required before deploying in production.

(a) ODDS

(b) ZeMA

(c) STRATH

Figure 5: Box plots of $U^{total}$ for correct predictions (TP+TN) and
erroneous predictions (FP+FN) using various models. (TP = true positives, TN =
true negatives, FP = true positives, FN = true negatives).

Rejecting uncertain predictions may yield higher predictive accuracy.
Strikingly, the detection performance can rise to > 90% $GSS$ after rejecting
40% of uncertain predictions for most tasks (Fig. 6). Without rejection,
however, for instance, in Fig. 6(c)(iii), most models scored much lower $GSS$
$(\approx 70\%)$. Note that not all models yield positive performance gains as
we reject predictions of high uncertainty, shown by a constant or
deteriorating trend, e.g. BAE-MCD on Fig. 6(a)(ii) and Fig. 6(b)(iv)).

(a) ODDS

(b) ZeMA

(c) STRATH

Figure 6: ARCs for comparisons of deterministic AE, VAE and BAEs with
posterior approximated under various techniques. Mean and standard error of
GSS are evaluated on (a) ODDS, (b) ZeMA and (c) STRATH datasets over 10
experiment runs. $U^{total}$ is used as the rejection criterion with the Q
distribution chosen based on the best $W^{GSS}$ score.

Should we use the $U^{exceed}$ or $\mathrm{Var_{\theta}[NLL]}$ method to
quantify anomaly uncertainty? Earlier, we mentioned that $U^{exceed}$ (Eq. 14)
does not scale well with the number of training examples. We notice it has a
poor mean score of 0.1% $Gain^{GSS}$ on the ZeMA dataset (Table 4(c)) while
having a better score on the STRATH dataset, which has fewer training examples
than ZeMA. Further inspection confirms the performance of $U^{exceed}$ is poor
and does not scale well with number of training examples Fig. 7. This lack of
scalability makes the $U^{exceed}$ less practical when there is an abundance
of training data, or when data is incrementally collected. On the other hand,
we do not observe such behaviours with the $U^{total}$ as it performs well
even as training examples increase in number.

Figure 7: Mean and standard error of $Gain^{GSS}$ against the number of
training examples in ODDS, ZeMA and STRATH datasets. Statistics are calculated
over 10 experiment runs and 5 types of AEs. Dotted line at $Gain^{GSS}=0$ is
the baseline.

All in all, $U^{total}$ outperforms $U^{exceed}$ and
$\mathrm{Var_{\theta}[NLL]}$. Despite its implementation simplicity, the
$\mathrm{Var_{\theta}[NLL]}$ is unreliable as a measure of anomaly uncertainty
as we find that $\mathrm{Var_{\theta}[NLL]}$ yields the lowest scores for both
ZeMA and STRATH datasets in Table 4(c). The $U^{total}$ also performs better
than $U^{epistemic}$ and $U^{aleatoric}$, proving the benefit of combining
both types of uncertainties compared to using either one of them.

Table 4: Mean $W^{GSS}$ ($Gain^{GSS}$) is evaluated evaluated on (a) ODDS, (b)
ZeMA and (c) STRATH datasets. Each dataset consists of multiple tasks and 10
experiment runs. Uncertainty estimation method with the highest mean $W^{GSS}$
score is bolded. Values are shown in percentage.

(a) ODDS

Model | $U^{epistemic}$ | $U^{aleatoric}$ | $U^{total}$ | $U^{exceed}$ | $\mathrm{Var_{\theta}[NLL]}$
---|---|---|---|---|---
Deterministic AE | - | 84.6(+4.3) | 84.6(+4.3) | 80.2(0.0) | -
VAE | 81.3(+1.8) | 83.8(+4.1) | 84.3(+4.6) | 79.4(-0.3) | 78.4(+0.1)
BAE-MCD | 81.4(+1.7) | 83.8(+3.8) | 84.2(+4.2) | 79.9(0.0) | 77.9(-1.3)
BAE-BBB | 86.7(+4.0) | 87.8(+4.9) | 89.2(+6.3) | 83.2(0.0) | 80.5(-0.7)
BAE-Ensemble | 85.2(+1.5) | 86.4(+3.0) | 88.0(+4.6) | 84.0(+0.3) | 81.2(-0.8)
Mean | 83.6(+2.2) | 85.3(+4.0) | 86.1(+4.8) | 81.3(0.0) | 79.5(-0.7)

(b) ZeMA

Model | $U^{epistemic}$ | $U^{aleatoric}$ | $U^{total}$ | $U^{exceed}$ | $\mathrm{Var_{\theta}[NLL]}$
---|---|---|---|---|---
Deterministic AE | - | 88.1(+2.9) | 88.1(+2.9) | 85.4(+0.4) | -
VAE | 89.4(+1.2) | 92.6(+5.0) | 92.9(+5.4) | 88.2(+0.2) | 89.4(+1.4)
BAE-MCD | 82.8(-0.1) | 86.1(+3.8) | 85.8(+3.6) | 83.2(+0.3) | 84.4(+1.8)
BAE-BBB | 89.8(+2.4) | 91.9(+3.5) | 92.6(+4.1) | 88.0(-0.4) | 86.5(-0.1)
BAE-Ensemble | 94.9(+1.3) | 95.1(+1.6) | 96.7(+3.2) | 93.5(0.0) | 94.0(+1.1)
Mean | 89.2(+1.2) | 90.8(+3.4) | 91.2(+3.8) | 87.7(+0.1) | 88.6(+1.0)

(c) STRATH

Model | $U^{epistemic}$ | $U^{aleatoric}$ | $U^{total}$ | $U^{exceed}$ | $\mathrm{Var_{\theta}[NLL]}$
---|---|---|---|---|---
Deterministic AE | - | 83.9(+6.3) | 83.9(+6.3) | 81.4(+2.5) | -
VAE | 82.6(+3.7) | 84.1(+6.2) | 84.2(+6.7) | 81.2(+4.1) | 81.0(+1.6)
BAE-MCD | 83.7(+5.7) | 84.2(+6.2) | 84.5(+6.5) | 80.2(+4.2) | 79.0(+0.3)
BAE-BBB | 82.6(+5.5) | 82.4(+5.8) | 83.2(+6.7) | 80.4(+6.4) | 78.3(+0.1)
BAE-Ensemble | 81.9(+2.5) | 84.9(+6.7) | 84.7(+6.5) | 81.8(+4.9) | 80.0(0.0)
Mean | 82.7(+4.4) | 83.9(+6.2) | 84.1(+6.5) | 81.0(+4.4) | 79.6(+0.5)

Deterministic AE vs stochastic AEs (VAE and BAE). We discuss differences
between AE models when using the $U^{total}$ to estimate the predictive
uncertainty. The deterministic AE, which quantifies only the aleatoric
component, stands as a baseline against the stochastic AEs. For most models,
the stochastic AEs have a higher $W^{GSS}$ than the deterministic AE, showing
the benefits of accounting for model uncertainty in a probabilistic framework.
However, this improvement does not hold for some models; one reason could be
the lower quality of posterior sampling [51, 33]. On the other hand, the BAE-
Ensemble outperforms the baseline AE and the VAE, scoring the highest
$W^{GSS}$ on ZeMA and STRATH, and second highest $W^{GSS}$ on ODDS.

Besides the performance gain, the BAE-Ensemble is relatively easy to implement
compared to other stochastic AEs. It does not require any changes to the
architecture or the optimiser of a deterministic AE, allowing a seamless
transition from a deterministic AE to a stochastic AE. While the BAE does
require $M$ times more computations than the deterministic AE, the
computations are highly parallelisable since we can compute each ensemble
member independently. Hence, the BAE can leverage the development of dedicated
hardware for NN computations, e.g. edge computing devices for industrial
environments [21].

Effect of anomaly probability conversion. Quantifying $U^{total}$ depends on
the method of anomaly probability conversion (e.g. using either the CDF or
ECDF, and whether to apply customised scaling); the switch from one method of
conversion to the other strongly influences the quality of uncertainty.
Notably, we find that customised scaling does improve the quality of
uncertainty (see Table 5) as we notice large gains in the percentage of
positive experiment runs after applying the customised scaling.

Table 5: Percentage of experiment runs where $Base^{GSS}$ >0 when using
$U^{total}$ while varying the method of anomaly probability conversion.
Percentage greater than 80% are in bold.

Conversion method | % (Positives/Total runs)
---|---
| ODDS | ZeMA | STRATH
ECDF | 18.2% (73/400) | 41.5% (83/200) | 48.5% (97/200)
ECDF + Scale | 83.0% (332/400) | 94.5% (189/200) | 89.5% (179/200)
Gaussian | 37.8% (151/400) | 47.5% (95/200) | 50.5% (101/200)
Gaussian + Scale | 80.5% (322/400) | 69.0% (138/200) | 90.0% (180/200)
Exponen. | 16.2% (65/400) | 48.5% (97/200) | 41.0% (82/200)
Exponen. + Scale | 82.5% (330/400) | 89.0% (178/200) | 44.0% (88/200)
Uniform | 45.8% (183/400) | 22.0% (44/200) | 50.5% (101/200)
Uniform + Scale | 50.0% (200/400) | 25.5% (51/200) | 61.5% (123/200)

In cases where $U^{total}$ did not perform as well, we suggest one reason is
the improper choice of distribution of NLL scores for the dataset. For
instance, although the exponential distribution worked reliably well for ODDS
and ZeMA, it did not work well for the STRATH dataset. Instead of committing
to a specific distribution, we can use the ECDF, a non-parametric method,
which works reliably well for all datasets. We report additional results in
Appendix C.

##  6 Conclusion

In this work, we have formulated the BAEs to estimate the predictive
uncertainty for anomaly detection. Our method computes the $U^{epistemic}$ and
$U^{aleatoric}$, and adding the two uncertainties forms the $U^{total}$.

Our experiments have evaluated the quality of uncertainty in the regime where
the model is given the option to reject predictions of high uncertainty, and
have demonstrated improvements using the proposed methods on multiple publicly
available datasets. The best performing BAE consistently outperform the
deterministic AE, highlighting the benefit of the Bayesian formulation.

Several directions for future work are possible. There is a need to further
understand the conditions in which the anomaly uncertainty may fail, for
instance, in adversarial attacks. Comprehensive case studies are necessary to
understand the practical challenges in deploying BAEs. Future work should also
explore the use of anomaly uncertainty for incremental and active learning.

## Acknowledgments

The work reported here was supported by the European Metrology Programme for
Innovation and Research (EMPIR) under the project Metrology for the Factory of
the Future (MET4FOF), project number 17IND12 and part-sponsored by Research
England’s Connecting Capability Fund award CCF18-7157: Promoting the Internet
of Things via Collaboration between HEIs and Industry (Pitch-In).

## References

  * [1] R. Chalapathy, S. Chawla, Deep learning for anomaly detection: A survey, arXiv preprint arXiv:1901.03407 (2019).
  * [2] M. Munir, M. A. Chattha, A. Dengel, S. Ahmed, A comparative analysis of traditional and deep learning-based anomaly detection methods for streaming data, in: 2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA), IEEE, 2019, pp. 561–566.
  * [3] G. Pang, C. Shen, L. Cao, A. V. D. Hengel, Deep learning for anomaly detection: A review, ACM Computing Surveys 54 (2) (mar 2021).  doi:10.1145/3439950.
URL https://doi.org/10.1145/3439950

  * [4] V. Chandola, A. Banerjee, V. Kumar, Anomaly detection: A survey, ACM computing surveys (CSUR) 41 (3) (2009) 1–58.
  * [5] A. D. Kiureghian, O. Ditlevsen, Aleatory or epistemic? does it matter?, Structural Safety 31 (2) (2009) 105–112, risk Acceptance and Risk Communication.
  * [6] A. Legrand, H. Trannois, A. Cournier, Use of uncertainty with autoencoder neural networks for anomaly detection, in: 2019 IEEE Second International Conference on Artificial Intelligence and Knowledge Engineering (AIKE), IEEE, 2019, pp. 32–35.
  * [7] B. X. Yong, T. Pearce, A. Brintrup, Bayesian autoencoders: Analysing and fixing the bernoulli likelihood for out-of-distribution detection, in: ICML Workshop on Uncertainty and Robustness in Deep Learning, 2020.
  * [8] E. Daxberger, J. M. Hernández-Lobato, Bayesian variational autoencoders for unsupervised out-of-distribution detection, arXiv preprint arXiv:1912.05651 (2019).
  * [9] B.-H. Tran, S. Rossi, D. Milios, P. Michiardi, E. V. Bonilla, M. Filippone, Model selection for bayesian autoencoders, arXiv preprint arXiv:2106.06245 (2021).
  * [10] C. Baur, B. Wiestler, S. Albarqouni, N. Navab, Bayesian skip-autoencoders for unsupervised hyperintense anomaly detection in high resolution brain mri, in: 2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI), IEEE, 2020, pp. 1905–1909.
  * [11] J. Gao, J. Yao, Y. Shao, Towards reliable learning for high stakes applications, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, 2019, pp. 3614–3621.
  * [12] M. S. A. Nadeem, J.-D. Zucker, B. Hanczar, Accuracy-rejection curves (arcs) for comparing classification methods with a reject option, in: Machine Learning in Systems Biology, PMLR, 2009, pp. 65–81.
  * [13] K. Hendrickx, L. Perini, D. Van der Plas, W. Meert, J. Davis, Machine learning with a reject option: A survey, arXiv preprint arXiv:2107.11277 (2021).
  * [14] C. Leibig, V. Allken, M. S. Ayhan, P. Berens, S. Wahl, Leveraging uncertainty information from deep neural networks for disease detection, Scientific reports 7 (1) (2017) 1–14.
  * [15] D. Kwon, H. Kim, J. Kim, S. C. Suh, I. Kim, K. J. Kim, A survey of deep learning-based network anomaly detection, Cluster Computing 22 (1) (2019) 949–961.
  * [16] M. Ahmed, A. N. Mahmood, M. R. Islam, A survey of anomaly detection techniques in financial domain, Future Generation Computer Systems 55 (2016) 278–288.  doi:https://doi.org/10.1016/j.future.2015.01.001.
URL https://www.sciencedirect.com/science/article/pii/S0167739X15000023

  * [17] T. Fernando, H. Gammulle, S. Denman, S. Sridharan, C. Fookes, Deep learning for medical anomaly detection – a survey, ACM Comput. Surv. 54 (7) (jul 2021).  doi:10.1145/3464423.
URL https://doi.org/10.1145/3464423

  * [18] P. Kamat, R. Sugandhi, Anomaly detection for predictive maintenance in industry 4.0-a survey, in: E3S Web of Conferences, Vol. 170, EDP Sciences, 2020, p. 02007\.
  * [19] S. Agrawal, J. Agrawal, Survey on anomaly detection using data mining techniques, Procedia Computer Science 60 (2015) 708–713.  doi:https://doi.org/10.1016/j.procs.2015.08.220.
URL https://www.sciencedirect.com/science/article/pii/S1877050915023479

  * [20] S. Thudumu, P. Branch, J. Jin, J. J. Singh, A comprehensive survey of anomaly detection techniques for high dimensional big data, Journal of Big Data 7 (1) (2020) 1–30.
  * [21] E. Li, L. Zeng, Z. Zhou, X. Chen, Edge ai: On-demand accelerating deep neural network inference via edge computing, IEEE Transactions on Wireless Communications 19 (1) (2019) 447–457.
  * [22] U. Bhatt, J. Antorán, Y. Zhang, Q. V. Liao, P. Sattigeri, R. Fogliato, G. Melançon, R. Krishnan, J. Stanley, O. Tickoo, et al., Uncertainty as a form of transparency: Measuring, communicating, and using uncertainty, in: Proceedings of the 2021 AAAI/ACM Conference on AI, Ethics, and Society, 2021, pp. 401–413.
  * [23] L. Perini, V. Vercruyssen, J. Davis, Quantifying the confidence of anomaly detectors in their example-wise predictions, in: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, Springer, 2020, pp. 227–243.
  * [24] C. Bishop, Pattern recognition and machine learning (2007).
  * [25] B. X. Yong, Y. Fathy, A. Brintrup, Bayesian autoencoders for drift detection in industrial environments, in: 2020 IEEE International Workshop on Metrology for Industry 4.0 & IoT, IEEE, 2020, pp. 627–631.
  * [26] R. Chandra, M. Jain, M. Maharana, P. N. Krivitsky, Revisiting bayesian autoencoders with mcmc, arXiv preprint arXiv:2104.05915 (2021).
  * [27] H.-P. Kriegel, P. Kroger, E. Schubert, A. Zimek, Interpreting and unifying outlier scores, in: Proceedings of the 2011 SIAM International Conference on Data Mining, SIAM, 2011, pp. 13–24.
  * [28] Y. Kwon, J.-H. Won, B. J. Kim, M. C. Paik, Uncertainty quantification using bayesian neural networks in classification: Application to biomedical image segmentation, Computational Statistics & Data Analysis 142 (2020) 106816.
  * [29] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016.
  * [30] T. Chen, E. Fox, C. Guestrin, Stochastic gradient hamiltonian monte carlo, in: International Conference on machine learning, 2014, pp. 1683–1691.
  * [31] Y. Gal, Z. Ghahramani, Dropout as a bayesian approximation: Representing model uncertainty in deep learning, in: International Conference on Machine Learning, 2016, pp. 1050–1059.
  * [32] C. Blundell, J. Cornebise, K. Kavukcuoglu, D. Wierstra, Weight uncertainty in neural network, in: International Conference on Machine Learning, PMLR, 2015, pp. 1613–1622.
  * [33] T. Pearce, F. Leibfried, A. Brintrup, Uncertainty in neural networks: Approximately bayesian ensembling, in: International conference on artificial intelligence and statistics, PMLR, 2020, pp. 234–244.
  * [34] D. P. Kingma, M. Welling, Auto-encoding variational bayes, arXiv preprint arXiv:1312.6114 (2013).
  * [35] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12 (2011) 2825–2830.
  * [36] H. G. Tucker, A generalization of the glivenko-cantelli theorem, The Annals of Mathematical Statistics 30 (3) (1959) 828–830.
  * [37] N. A. Weiss, A course in probability, Pearson Addison-Wesley, 2005.
  * [38] L. I. Kuncheva, Á. Arnaiz-González, J.-F. Díez-Pastor, I. A. Gunn, Instance selection improves geometric mean accuracy: a study on imbalanced data classification, Progress in Artificial Intelligence 8 (2) (2019) 215–228.
  * [39] S. Rayana, ODDS library (2016).
URL http://odds.cs.stonybrook.edu

  * [40] G. O. Campos, A. Zimek, J. Sander, R. J. Campello, B. Micenková, E. Schubert, I. Assent, M. E. Houle, On the evaluation of unsupervised outlier detection: measures, datasets, and an empirical study, Data mining and knowledge discovery 30 (4) (2016) 891–927.
  * [41] J. Chen, S. Sathe, C. Aggarwal, D. Turaga, Outlier detection with autoencoder ensembles, in: Proceedings of the 2017 SIAM international conference on data mining, SIAM, 2017, pp. 90–98.
  * [42] T. Schneider, S. Klein, M. Bastuck, Condition monitoring of hydraulic systems Data Set at ZeMA (Apr. 2018).  doi:10.5281/zenodo.1323611.
URL https://doi.org/10.5281/zenodo.1323611

  * [43] C. Tachtatzis, G. Gourlay, I. Andonovic, O. Panni, Sensor data set radial forging at afrc testbed v2 (Sep. 2019).  doi:10.5281/zenodo.3405265.
URL https://doi.org/10.5281/zenodo.3405265

  * [44] N. Helwig, E. Pignanelli, A. Schütze, Condition monitoring of a complex hydraulic system using multivariate statistics, in: 2015 IEEE International Instrumentation and Measurement Technology Conference (I2MTC) Proceedings, 2015, pp. 210–215.  doi:10.1109/I2MTC.2015.7151267.
  * [45] Y. Luo, P. Harris, Uncertainty in data analysis for strath testbed, in: 2020 IEEE International Workshop on Metrology for Industry 4.0 & IoT, IEEE, 2020, pp. 95–100.
  * [46] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [47] L. N. Smith, Cyclical learning rates for training neural networks, in: 2017 IEEE Winter Conference on Applications of Computer Vision (WACV), IEEE, 2017, pp. 464–472.
  * [48] O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation, in: N. Navab, J. Hornegger, W. M. Wells, A. F. Frangi (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, Springer International Publishing, Cham, 2015, pp. 234–241.
  * [49] A. L. Maas, A. Y. Hannun, A. Y. Ng, Rectifier nonlinearities improve neural network acoustic models, in: ICML Workshop on Deep Learning for Audio, Speech and Language Processing, 2013.
  * [50] A. Banerjee, U. Chitnis, S. Jadhav, J. Bhawalkar, S. Chaudhury, Hypothesis testing, type i and type ii errors, Industrial psychiatry journal 18 (2) (2009) 127.
  * [51] J. Yao, W. Pan, S. Ghosh, F. Doshi-Velez, Quality of uncertainty quantification for bayesian neural network inference, arXiv preprint arXiv:1906.09686 (2019).
  * [52] S. Kullback, R. A. Leibler, On information and sufficiency, The annals of mathematical statistics 22 (1) (1951) 79–86.
  * [53] J. W. Tukey, et al., Exploratory data analysis, Vol. 2, Reading, Mass., 1977.
  * [54] V. Bewick, L. Cheek, J. Ball, Statistics review 13: receiver operating characteristic curves, Critical care 8 (6) (2004) 1–5.

##  Appendix A BAE posterior sampling methods

We describe the following methods to sample from the BAE: BBB [32], MCD [31]
and anchored ensembling [33].

###  A.1 Variational inference

BBB and MCD are examples of variational inference. The idea of variational
inference is to approximate the posterior distribution over parameters by
introducing a variational distribution.

$q_{\psi}{}(\theta)\approx p(\theta|X)$ (18)

where $\psi$ is the vector of parameters of the variational distribution, $q$.
The objective of training is thus, to minimise the Kullback-Leibner (KL)
divergence [52], which is a measure of similarity between two distributions:
the variational distribution and the true posterior. This yields an optimised
variational distribution to be sampled from during prediction.

$\psi^{\mathrm{opt}}=\operatorname*{arg\,min}_{\psi}\mathrm{KL}[q_{\psi}{}(\theta)\,||\,p(\theta|X)]$
(19)

The KL divergence can be approximated with $M$ samples of AE parameters
$\\{\theta_{m}\\}_{m=1}^{M}$ from the variational distribution $q$. Note that
minimising the KL divergence is the same as maximising the log evidence lower
bound (ELBO).

$\mathrm{KL}[q_{\psi}{}(\theta)\,||\,p(\theta|X)]\approx\sum_{m=1}^{M}-\log{p(X|\theta_{m})}-\lambda\,\mathrm{KL}[q_{\psi}{}(\theta_{m})\,||\,p(\theta_{m})]$
(20)

where the weight decay, $\lambda$ scales the KL divergence of the prior,
effectively controlling the strength of regularisation. Note that the first
term of the sum is the log-likelihood of a sample from the variational
distribution.

Now, we discuss popular models of the variational distribution. In BBB, we use
a diagonal Gaussian distribution as the variational distribution, while the
prior is a mixture of 2 diagonal Gaussian distributions,

$\displaystyle\psi=(\mu,\sigma^{2})$
$\displaystyle\epsilon\sim\mathrm{N}(0,1)$
$\displaystyle\theta=\mu+\sigma\epsilon$

where $\epsilon$ is the noise term introduced for implementing the
reparameterisation trick [32], necessary for backpropagating the gradients
during optimisation. The KL prior loss is

$\begin{split}\mathrm{KL}[q_{\psi}{}(\theta_{m})\,||\,p(\theta_{m})]^{\mathrm{BBB}}=\pi(\frac{||\theta_{m}||^{2}}{2\tau_{1}^{2}}-\log{\tau_{1}})\,+(1-\pi)(\frac{||\theta_{m}||^{2}}{2\tau_{2}^{2}}-\log{\tau_{2})}\,+\\\
\frac{||\theta_{m}-\mu||^{2}}{2\sigma^{2}}-\log{\sigma}\end{split}$ (21)

where $\pi$ , $\tau_{1}$ and $\tau_{2}$ are parameters of the prior, which we
fix in our experiments as 0.5, 1.0 and 0.1, respectively.

On the other hand, the MCD, which is also type of variational inference, is
implemented by adding a dropout layer after each linear or convolutional
layer. The corresponding variational distribution is a mixture of Gaussian and
Bernoulli distributions.

$\displaystyle\psi=W$
$\displaystyle\epsilon\sim\mathrm{Ber}(1-p^{\mathrm{dropout}})$
$\displaystyle\theta=W\cdot\epsilon$

where $p^{\mathrm{dropout}}$ is the dropout probability and $W$ refers to the
AE parameters when they are not dropped out. Further, the
$p^{\mathrm{dropout}}$ is a hyperparameter and is not updated during training.
We set $p^{\mathrm{dropout}}=0.01$ in our experiments. Lastly, the KL prior
term can be approximated as

$\mathrm{KL}[q_{\psi}{}(\theta_{m})\,||\,p(\theta_{m})]^{\mathrm{MCD}}\approx||\theta_{m}||^{2}$
(22)

###  A.2 Anchored ensembling

In anchored ensembling, posteriors are approximated by Bayesian inference
under the family of methods called randomised maximum a posteriori (MAP)
sampling, where model parameters are regularised by values drawn from a
distribution (so-called anchor distribution), which can be set equal to the
prior.

Assume our ensemble consists of $M$ independent autoencoders and each j-th
autoencoder contains a set of parameters, $\theta_{m}$ where
$m\in{\\{1,2,3...M\\}}$. In anchored ensembling for approximating the
posterior distribution, the ‘anchored weights‘ for each autoencoder are unique
and sampled during initialisation from a prior distribution
$\theta_{m}^{\mathrm{anc}}\sim
N(\mu_{m}^{\mathrm{anc}},\sigma_{m}^{\mathrm{anc}^{2}})$ and remain fixed
throughout the training procedure.

The autoencoders are trained by minimising the loss function, which is the
negative sum of log-likelihood (based on i.i.d assumption) and log-prior where
both are assumed to be Gaussian. For each member of the ensemble, the loss to
be optimised is

$\mathcal{L}(\theta_{m},X)=-\log{p(X|\theta_{m})}+\lambda||\theta_{m}-\theta_{m}^{\mathrm{anc}}||^{2}$
(23)

where $\lambda$ is a hyperparameter for scaling the regulariser term arising
from the prior.

##  Appendix B Data preprocessing

ODDS. For each dataset in the ODDS collection, we do not apply any
preprocessing steps aside from min-max scaling and data splitting.

For ZeMA and STRATH, the inputs to the AE have the dimensions of $(B\times
L\times K\,)$, where $B$ is the batch size, $L$ is the sequence length, and
$K$ is the number of sensors.

ZeMA. For each datapoint in the batch, The labels of each hydraulic
subsystem’s condition are available in the dataset, and we label the best
working conditions as inliers while the remaining states as anomalies. We have
chosen the following sensor-subsystem pairs in our analysis: temperature
sensor (TS4) for the cooler, valve and accumulator, and pressure sensor (PS6)
for the pump. We downsample the pressure sensor to 1Hz; no resampling is
applied on the temperature data, hence $L=60$ and $K=1$.

STRATH. The radial forging process is divided into the heating and forging
phases. We segment the sensors data to consider only the forging phase. We
downsample the data by 10 fold, reducing the sequence length to $L=562$ for
all tasks. For tasks (i-iii), $K=1$ and $K=3$ for task (iv). For each forged
part, measurements of its geometric dimensions are available as quality
indicators; we focus our analysis on the 38 diameter@200 dimension. To flag
parts as anomalies, we first obtain the absolute difference between the
measured and nominal dimensions, and subsequently, we apply the Tukey’s fences
method [53]. The remaining parts are labelled as inliers. The data
preprocessing steps are outlined in Fig. 8.

Figure 8: Data preprocessing pipeline for STRATH data set

As mentioned, we apply min-max scaling for all datasets. For ODDS, we apply
the min-max scaling to each feature independently. By contrast, for ZeMA and
STRATH, we obtain the min-max values from the train set for each sensor
independently, instead of each feature, to retain the shape of signal. Note
that we prevent train-test bias when by fitting the scaler to the train set
only instead of the entire dataset.

##  Appendix C Additional results

Results of using the area under the receiver operating characteristic curve
(AUROC) [54] as a measure of classifier performance for ARCs are reported in
Fig. 9 and Table 6(c). In addition, we report the ARCs to compare the choices
of distribution Q (see Fig. 10).

(a) ODDS

(b) ZeMA

(c) STRATH

Figure 9: ARCs for comparisons of deterministic AE, VAE and BAEs with
posterior approximated under various techniques. Mean and standard error of
AUROC are evaluated on (a) ODDS, (b) ZeMA and (c) STRATH datasets over 10
experiment runs. $U^{total}$ is used as the rejection criterion with the Q
distribution chosen based on the best $W^{AUROC}$ score. Table 6: Mean
$W^{AUROC}$ ($Gain^{AUROC}$) is evaluated on (a) ODDS, (b) ZeMA and (c) STRATH
datasets. Each dataset consists of multiple tasks and 10 experiment runs.
Uncertainty estimation method with the highest mean $W^{AUROC}$ score is
bolded. Values are shown in percentage.

(a) ODDS

Model | $U^{epistemic}$ | $U^{aleatoric}$ | $U^{total}$ | $U^{exceed}$ | $\mathrm{Var_{\theta}[NLL]}$
---|---|---|---|---|---
Deterministic AE | - | 91.6(+1.8) | 91.6(+1.8) | 90.1(+0.2) | -
VAE | 89.0(+1.2) | 89.8(+2.0) | 90.0(+2.2) | 88.0(+0.2) | 86.4(-1.3)
BAE-MCD | 91.5(+0.8) | 92.6(+1.9) | 92.6(+2.0) | 90.5(-0.2) | 88.3(-2.4)
BAE-BBB | 92.8(+1.4) | 93.5(+2.1) | 93.9(+2.5) | 91.6(+0.3) | 89.7(-1.7)
BAE-Ensemble | 93.2(+0.9) | 93.8(+1.5) | 94.3(+2.0) | 92.4(+0.2) | 90.2(-2.1)
Mean | 91.6(+1.1) | 92.3(+1.9) | 92.5(+2.1) | 90.5(+0.1) | 88.6(-1.9)

(b) ZeMA

Model | $U^{epistemic}$ | $U^{aleatoric}$ | $U^{total}$ | $U^{exceed}$ | $\mathrm{Var_{\theta}[NLL]}$
---|---|---|---|---|---
Deterministic AE | - | 92.7(+2.6) | 92.7(+2.6) | 90.8(+0.6) | -
VAE | 94.7(+1.3) | 95.7(+2.3) | 95.8(+2.4) | 94.0(+0.5) | 94.9(+1.4)
BAE-MCD | 91.7(+1.8) | 92.8(+3.0) | 92.9(+3.0) | 90.4(+0.6) | 92.2(+2.4)
BAE-BBB | 96.0(+1.4) | 96.3(+1.6) | 96.5(+1.8) | 95.3(+0.6) | 95.1(+0.4)
BAE-Ensemble | 97.5(+1.5) | 97.5(+1.6) | 97.8(+1.8) | 96.2(+0.3) | 97.0(+1.0)
Mean | 95.0(+1.5) | 95.0(+2.2) | 95.1(+2.3) | 93.3(+0.5) | 94.8(+1.3)

(c) STRATH

Model | $U^{epistemic}$ | $U^{aleatoric}$ | $U^{total}$ | $U^{exceed}$ | $\mathrm{Var_{\theta}[NLL]}$
---|---|---|---|---|---
Deterministic AE | - | 89.4(+3.1) | 89.4(+3.1) | 88.5(+2.2) | -
VAE | 89.4(+2.9) | 89.6(+3.1) | 89.6(+3.1) | 88.7(+2.2) | 88.2(+1.7)
BAE-MCD | 89.3(+3.0) | 89.5(+3.2) | 89.4(+3.1) | 88.2(+1.9) | 86.9(+0.6)
BAE-BBB | 88.9(+3.5) | 88.6(+3.2) | 88.4(+3.1) | 87.6(+2.2) | 86.8(+1.5)
BAE-Ensemble | 89.1(+2.3) | 89.8(+3.0) | 89.8(+3.0) | 89.0(+2.2) | 87.7(+0.9)
Mean | 89.2(+2.9) | 89.4(+3.1) | 89.3(+3.1) | 88.4(+2.1) | 87.4(+1.2)

(a) ODDS

(b) ZeMA

(c) STRATH

Figure 10: ARCs for comparisons of anomaly probability conversions. Mean and
standard error of GSS are evaluated on (a) ODDS, (b), ZeMA and (c) STRATH
datasets over 10 experiment runs. Results are shown for the BAE-Ensemble model
using $U^{total}$ as the rejection criterion.