# Cosmology in the mimetic higher-curvature $f(R,R_{\mu\nu}R^{\mu\nu})$
gravity

Adam Z. Kaczmarek1<EMAIL_ADDRESS>Dominik Szczȩśniak2
1Department of Physics, Czȩstochowa University of Technology, 19 Armii
Krajowej Ave., 42200 Czȩstochowa, Poland  2 Department of Theoretical Physics,
Faculty of Science and Technology, Jan Długosz University in Czȩstochowa,
13/15 Armii Krajowej Ave., 42200 Czȩstochowa, Poland

###### Abstract

In the framework of the mimetic approach, we study the
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity with the Lagrange multiplier constraint
and the scalar potential. We introduce field equations for the discussed
theory and overview their properties. By using the general reconstruction
scheme we obtain the power law cosmology model for the
$f(R,R_{\mu\nu}R^{\mu\nu})=R+d(R_{\mu\nu}R^{\mu\nu})^{p}$ case as well as the
model that describes symmetric bounce. Moreover, we reconstruct model,
unifying both matter dominated and accelerated phases, where ordinary matter
is neglected. Using inverted reconstruction scheme we recover specific
$f(R,R_{\mu\nu}R^{\mu\nu})$ function which give rise to the de-Sitter
evolution. Finally, by employing the perfect fluid approach, we demonstrate
that this model can realize inflation consistent with the bounds coming from
the BICEP2/Keck array and the Planck data. Thus, it is suggested that the
introduced extension of the mimetic regime may describe any given cosmological
model.

##  I Introduction

In the 1990s, the field of cosmology encountered serious breakthrough with the
discovery of accelerated expansion of the Universe. This led to the current
evolutionary status that is determined by the three contributions to the total
density parameter ($\Omega$): the ordinary matter ($\Omega_{m}\approx 5\%$),
the dark matter ($\Omega_{DM}\approx 27\%$) and the dark energy
($\Omega_{DE}\approx 68\%$), where the latter contribution is responsible for
the accelerated expansion phenomena [1]. Importantly, the nature of the dark
sector in the energetic budget is also the reason for intense experimental and
theoretical studies that exposed shortcomings of the general relativity (GR)
[2, 3]. In this respect, one of the most promising proposals to deal with the
mentioned issues is based on the extensions and modifications of the
Einstein’s framework of GR [4]. Note that the proposals of this type are quite
successful not only in the context of the aforementioned cosmic expansion but
also in other issues of the modern cosmology [5]. Moreover, in terms of the
quantum gravity, the standard GR cannot be renormalized but higher-curvature
terms in the gravitational action may tackle this issue [6, 7, 8].

In the literature the GR has been extended in many ways, mainly by including
additional curvature terms, scalar fields or coupling with matter. For a well-
written reviews on the modified gravities we refer reader to [9, 10, 11, 4].
Among many interesting approaches, one of the most popular proposals is the
so-called $f(R)$ gravity introduced by Buchdachl in 70’s [12, 13, 14], where
standard GR action is generalized to the arbitrary function of the Ricci
scalar ($R\rightarrow f(R)$). Further curvature invariants were included in
the $f(\mathcal{G})$ gravity, where arbitrary function of the Gauss-Bonnet
invariant is added to the Lagrangian [15]. In a wide class of the higher-
curvature models, the $f(R,\mathcal{G})$, $f(R,R_{\mu\nu}R^{\mu\nu})$ or
theories with the higher-order of derivatives have also been introduced and
studied extensively [16, 17, 18, 19, 20, 18]. The authors of [21] considered
combinations of the Ricci scalar, the Ricci squared ($R_{\mu\nu}R^{\mu\nu}$)
and the Riemann squared ($R_{\mu\nu\alpha\beta}R^{\mu\nu\alpha\beta}$) terms
and shown that such models can provide alternative to the dark energy.
Moreover, theories of this type were recently studied in the context of
divergences and renormalization group frameworks and play important role in
the perturbative approach to the quantum gravity [22, 8, 23]. Note that the
general relativity can be further extended by the scalar-geometry or scalar-
matter couplings. In this case gravitational action includes terms with
contributions not only from geometry but also various functions of the scalar
field, matter Lagrangian or trace of the energy-momentum tensor.
Representative example of the scalar-geometry coupling is the Jordan-Brans-
Dicke theory [24, 25], further generalized to the $f(\phi,R)$ and
$f(\phi,\mathcal{G})$ gravity [9]. For the matter-geometry couplings, the
$f(R,T)$ and $f(R,T,R_{\mu\nu}T^{\mu\nu})$ models ($T$ denotes trace of the
energy-momentum tensor $T_{\mu\nu}$) with further extensions are good examples
studied widely in the literature [26, 27, 28, 29, 30, 31, 32].

On the other hand, the dark matter is even more mysterious since very few
information is known except the fact that it modifies galactic rotation, it is
not baryonic, has no electric and color charges and do not interact with the
light [33, 34]. Apart from the particle candidates for the dark matter, such
as supersymmetry neutralino, an interesting candidate has been provided by
Chamseddine and Mukhanov. In their so-called mimetic modification of the GR,
the leading idea was to isolate gravity’s conformal degree of freedom by
parametrizing physical metric ($g_{\mu\nu}$) in terms of the scalar field
($\phi$) and the auxiliary metric ($\hat{g}_{\mu\nu}$) [35]. As a consequence,
the conformal degree of freedom turns out to be dynamical, even in the absence
of the standard matter fluids, and may mimic behavior of the cold dark matter.
It is worth to remark here, that the mimetic theory has been recently
extensively studied in various contexts [36, 37, 38, 39, 40, 41, 42, 43, 44].
Moreover, the mimetic approach has been combined with other modifications of
gravity, such as $f(R)$, $f(\mathcal{G})$ or $f(R,T)$ theories with emphasis
on the cosmological applications [45, 46, 47, 48, 49, 50, 51, 52]. The
unimodular extensions of the mimetic approach were also introduced and studied
in recent articles [53, 54]. For introduction and review of the mimetic theory
see [55].

Besides the riddle of the dark Universe, the modified gravity is also
important in terms of the early stages of the Universe. In this context, the
theory of inflation is by now most possible scenario that describes early
Universe history [56, 57]. This framework, developed in the 80’s, provides
explanation of important issues in the Big Bang cosmology such as the horizon,
flatness and monopole problems [58]. According to the inflation paradigms,
early stages of the Universe exhibited exponentially accelerated expansion,
which after deceleration led to the standard cosmological eras [59]. However,
no direct proof that this early epoch occurred is given up to date. One should
note that inflation is well described by the slow-roll models which usually
contain single scalar field (inflaton) that drives inflation. Moreover, the
form of the scalar field potential can be related to the spectral index as
well as to the scalar and tensor perturbations as generated during the
inflationary epoch. Thus, scalar potential may be consistent with
observations. Recent observational data coming from Planck and BICEP2/Keck
array heavily constrained scalar-to-tensor ratio and spectral index of
primordial curvature perturbations corresponding to the inflation [60]. One
should note, that the bounce paradigm may also be viable characteristic of the
early Universe and an interesting alternative [61, 62, 63, 64]. From the
modified gravity perspective, inflation can be successfully described and the
unification of early-time and late-time acceleration eras is also possible
[10]. In the literature, inflation in the modified gravity is still heavily
studied and debated topic. In particular, for the $f(R)$ gravities a
successful model of the Starobinsky’s inflation exists, which is consistent
with the local and large scale constraints of gravity [65]. Alternatively,
other modifications of the GR were studied in the context of inflation as well
as the unification of accelerated phases of the Universe [5, 66]. Note that
inflation for the mimetic gravity with its extensions, was also recently
investigated [47, 67, 68, 19, 43].

In the context of the above, herein we combine mimetic idea with the general
higher-curvature $f(R,R_{\mu\nu}R^{\mu\nu})$ gravity. We intent to extend
previously considered mimetic framework with more general
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity, with respect to the cosmological
reconstruction of various models. After brief discussion of the field
equations, by using the reconstruction method, we obtain the power law
evolution model, the symmetric bounce as well as the transition between
radiation and matter dominated phases. Then, in reference the current
experimental data, we discuss inflationary model of the mimetic
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity. By using the reconstruction methods, we
obtain Lagrange multiplier and potential responsible for inflationary
evolution described by the scale factor $H(N)=\alpha e^{\beta N}+\gamma$.
Finally by employing the perfect fluid approach, we show that it is possible
to achieve inflationary evolution in the discussed theory, in agreement with
the current data.

##  II Theoretical Model

The main idea behind the mimetic gravity is that the $g_{\mu\nu}$ metric (the
fundamental variable in gravity) may be expressed by the new degrees of
freedom [19]. By doing so, such degrees of freedom can admit wider class of
solutions than the standard GR. In particular, the conformal degree of freedom
can be isolated by expressing the physical metric in terms of the scalar field
$\phi$ and the auxiliary metric $\hat{g}_{\mu\nu}$ as follows:

$\displaystyle
g_{\mu\nu}=-\hat{g}^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi\hat{g}_{\mu\nu}.$
(1)

Note that the above parametrization posses Weyl symmetry under the conformal
transformation, i.e. $\hat{g}_{\alpha\beta}=e^{\omega(x)}g_{\alpha\beta}$. It
is also worth to remark, that instead of using the physical metric
$g_{\mu\nu}$ in variation of the gravitational action, one can use variation
with respect to the metric $\hat{g}_{\mu\nu}$ and the scalar field $\phi$. As
a result, from Eq. (1) the following constraint can be obtained:

$\displaystyle
g^{\mu\nu}(\phi,\hat{g}_{\mu\nu})\partial_{\mu}\phi\partial_{\nu}\phi=-1.$ (2)

In order to impose given constraint at the action level we use the Lagrange
multiplier formalism [69, 70]. In the mimetic regime we need to introduce the
Lagrange multiplier ($\lambda$) that corresponds to the mimetic constraint
[19]. Hence, for the mimetic $f(R,R^{\mu\nu}R_{\mu\nu})$ gravity we have:

$\displaystyle S=\int
d^{4}x\Big{[}\sqrt{-g}\big{(}f(R,R_{\mu\nu}R^{\mu\nu})-V(\phi)+\lambda(g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+1)\big{)}\Big{]},$
(3)

where $f(R,R^{\mu\nu}R_{\mu\nu})$ is the analytical function of the Ricci
scalar $R$ and contraction of the Ricci tensors (the Ricci squared term)
$R^{\mu\nu}R_{\mu\nu}$. Moreover, $V(\phi)$ denotes scalar (mimetic)
potential. We note that the equations obtained from the variation with respect
to the physical metric $g_{\mu\nu}$ with imposed mimetic constraint are fully
equivalent to the equations that one can derive by using action written in
terms of the auxiliary metric $\hat{g}_{\mu\nu}$. We emphasize that standard
matter fields are not part of our considerations since one of the core
characteristics of the mimetic theory is to mimic matter content [19]. In what
follows, the variation of the action (3) with respect to the components of the
metric tensor $g^{\mu\nu}$, gives the following field equation:

$\displaystyle f_{R}R_{\mu\nu}$
$\displaystyle-\nabla_{\mu}\nabla_{\nu}f_{R}+g_{\mu\nu}\Box
f_{R}+\Box(R_{\mu\nu}f_{Y})+g_{\mu\nu}\nabla^{\alpha}\nabla^{\beta}(R_{\alpha\beta}f_{Y})-2\nabla_{\alpha}\nabla_{\beta}(R^{\alpha}_{(\mu}\delta^{\beta}_{\nu)}f_{Y})+2R^{\alpha}_{\mu}R_{\alpha\nu}f_{Y}$
$\displaystyle-\frac{1}{2}g_{\mu\nu}(f-V(\phi)+\lambda(g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi+1))+\lambda\partial_{\mu}\phi\partial_{\nu}\phi=0,$
(4)

where $Y=R_{\mu\nu}R^{\mu\nu}$, $f_{Y}=\frac{\partial f(R,Y)}{\partial R}$ and
$f_{Y}=\frac{\partial f(R,Y)}{\partial Y}$, whereas the indices in brackets
denote symmetrization ($A_{(\mu\nu)}=\frac{1}{2}(A_{\mu\nu}+A_{\nu\mu})$).
Next, by varying the action with respect to the scalar field $\phi$ the
following scalar equation can be obtained:

$\displaystyle 2\nabla^{\mu}(\lambda\partial_{\mu}\phi)+V^{\prime}(\phi)=0,$
(5)

where prime stands for the differentiation with respect to the scalar field
($V^{\prime}(\phi)=\frac{dV(\phi)}{d\phi}$). It is important to remark, that
the variation of the action with respect to $\lambda$ returns mimetic
constraint.

Further in this work we assume that the background geometry is described by
the FLRW metric, with the following line element:

$\displaystyle ds^{2}=-dt^{2}+a^{2}(t)dx_{i}dx^{i},$ (6)

where $i=1,2,3$ and $a(t)$ is the scale factor. For simplicity, we also assume
that the scalar field depends only on time, i.e. $\phi=\phi(t)$. By using the
FLRW metric (6), the corresponding $(t,t)$ component of the mimetic
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity equations (4) is:

$\displaystyle-3(\dot{H}+H^{2})f_{R}$
$\displaystyle+3H\partial_{t}f_{R}+\frac{1}{2}(f+\lambda(\dot{\phi}^{2}+1)-V(\phi))+6(2H\ddot{H}-2\dot{H}^{2}+3H^{2}\dot{H}-3H^{4})f_{Y}$
$\displaystyle+6(2H\dot{H}+3H^{3})\partial_{t}f_{Y}=0,$ (7)

while the $(i,j)$ component reads:

$\displaystyle(\dot{H}+3H^{2})f_{R}-2H\partial_{t}f_{R}-\partial_{tt}f_{R}-\frac{1}{2}(f-\lambda(\dot{\phi}^{2}-1)-V(\phi))-2(2\dddot{H}+12H\ddot{H}-9H^{4}+9H^{2}\dot{H}+6\dot{H}^{2})f_{Y}$
$\displaystyle-2(18H\dot{H}+6H^{3}+4\ddot{H})\partial_{t}f_{Y}-2(2\dot{H}+3H^{2})\partial_{tt}f_{Y}=0.$
(8)

From Eq. (5) the corresponding scalar equation is given as:

$\displaystyle
6H\lambda\dot{\phi}+2(\dot{\lambda}\dot{\phi}+\lambda\ddot{\phi})-V^{\prime}(\phi)=0.$
(9)

Note, that dot denotes differentiation with respect to the cosmic time $t$ and
$H=\dot{a}/a$ is the Hubble rate. For brevity, we note that the Ricci scalar
and the Ricci tensor squared in the FLRW spacetime are:

$\displaystyle
R=6(\dot{H}+2H^{2}),\;\;\;\;Y=R_{\mu\nu}R^{\mu\nu}=12\big{(}\dot{H}^{2}+3H^{2}(\dot{H}+H^{2})\big{)}.$
(10)

Finally, the mimetic constraint for the metric given by Eq. (6) is:

$\displaystyle\dot{\phi}^{2}=1,$ (11)

and yields:

$\displaystyle\phi=t+C,$ (12)

with $C$ being integration constant. Keeping in mind that the scalar field
appears only with derivatives, shift symmetry allows for setting constant
$C=0$ without losing generality. By using $\phi=t$ the field and scalar
equations (7-9) reduce to:

$\displaystyle-3(\dot{H}+H^{2})f_{R}$
$\displaystyle+3H\partial_{t}f_{R}+\frac{1}{2}(f-V(\phi))+6(2H\ddot{H}-2\dot{H}^{2}+3H^{2}\dot{H}-3H^{4})f_{Y}$
$\displaystyle+6(2H\ddot{H}+3H^{3})\partial_{t}f_{Y}+\lambda=0,$ (13)
$\displaystyle(\dot{H}+3H^{2})f_{R}-2H\partial_{t}f_{R}-\partial_{tt}f_{R}-\frac{1}{2}(f-V(\phi))-2(2\dddot{H}+12H\ddot{H}-9H^{4}+9H^{2}\dot{H}+6\dot{H}^{2})f_{Y}$
$\displaystyle-2(18H\dot{H}+6H^{3}+4\ddot{H})\partial_{t}f_{Y}-2(2\dot{H}+3H^{2})\partial_{tt}f_{Y}=0,$
(14)

and

$\displaystyle 6H\lambda+2\dot{\lambda}-V^{\prime}(\phi)=0.$ (15)

Hence, Eq. (14) can be used to determine potential $V(\phi)$ that generates
evolution of the Hubble parameter for the specific $f(R,Y)$ gravity:

$\displaystyle V(\phi=t)=$
$\displaystyle-2(\dot{H}+3H^{2})f_{R}+4H\partial_{t}f_{R}+2\partial_{tt}f_{R}+f+4\big{[}(2\dddot{H}+12H\ddot{H}-9H^{4}+9H^{2}\dot{H}+8\dot{H}^{2})f_{Y}$
$\displaystyle+(18H\dot{H}+6H^{3}+4\ddot{H})\partial_{t}f_{Y}+(2\dot{H}+3H^{2})\partial_{tt}f_{Y}\big{]}.$
(16)

By combining Eqs. (13) and (14), the analytic form of the Lagrange multiplier
$\lambda(t)$ is provided:

$\displaystyle\lambda(t)=2\dot{H}f_{R}-H\partial_{t}f_{R}+\partial_{tt}f_{R}+4(3H\ddot{H}+\dddot{H}+\dot{H}^{2})f_{Y}-2(3H^{3}-12H\dot{H}-4\ddot{H})\partial_{t}f_{Y}+2(2\dot{H}+3H^{2})\partial_{tt}f_{Y}.$
(17)

This equation allows to obtain the desired evolution in various ways. Given
the Hubble parameter $H(t)$ for the specified $f(R,Y)$ function, one can find
corresponding multiplier $\lambda$ and mimetic potential $V(\phi)$. Also,
given the form of the mimetic potential and the Lagrange multiplier $\lambda$,
one can solve above equations to find proper $f(R,R_{\mu\nu}R^{\mu\nu})$
gravity model that realizes cosmological scenario of interest.

##  III Reconstruction of the FLRW cosmologies

In the following section, we use reconstruction methods in order to obtain
various cosmological scenarios in the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$
gravity. The procedure of cosmological reconstruction is widely used in the
context of modified theories of gravity and helps with the extraction of
physical implications in the theory of interest [67, 71, 72, 73, 74, 75, 76,
32]. In general, obtaining solutions for the GR’s extensions may be
troublesome, since even small modifications can drastically increase
difficulty of equations. The reconstruction scheme used in the extensions of
the GR relies on the fact that arbitrary function (or potential) is used in
the definition of the modified gravity. Once the Hubble factor $H$ is
specified, the field equations are solved to obtain the model of interest.
This technique is a good way to survey modified theories of gravity, since
satisfying well established cosmological models in gravity’s extensions is
desired. Moreover, this technique yields difference between the GR and its
modifications, namely specific class (or classes) of the extended gravity may
describe any given spacetime [10]. In this section, once the function $f(R,Y)$
is given, we obtain the Lagrange multiplier $\lambda$ and potential $V$
satisfying specified evolution. In general, by using this technique one can
achieve model describing any given Hubble rate such as inflation or
transitions between phases for the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$
approach.

In the present study, we assume specific form of the
$f(R,R_{\mu\nu}R^{\mu\nu})$ function, where Ricci squared term enters with
arbitrary power $p$, i. e.:

$\displaystyle f(R,Y)=R+dY^{p}.$ (18)

Thus, the rearranged equations for this choice of the $f(R,Y)$ gravity are:

$\displaystyle
V(t)=-2(\dot{H}+3H^{2})f_{R}+4H\partial_{t}f_{R}+2\partial_{tt}f_{R}+f+4dp\big{[}(2\dddot{H}+12H\ddot{H}-9H^{4}+9H^{2}\dot{H}+8\dot{H}^{2})Y^{p-1}$
$\displaystyle+(18H\dot{H}+6H^{3}+4\ddot{H})(p-1)Y^{p-2}\dot{Y}+(2\dot{H}+3H^{2})((p-1)(p-2)Y^{p-3}\dot{Y}^{2}+(p-1)Y^{p-2}\ddot{Y})\big{]},$
(19)

and

$\displaystyle\lambda(t)=2\dot{H}f_{R}-H\partial_{t}f_{R}+\partial_{tt}f_{R}+4dp(3H\ddot{H}+\dddot{H}+\dot{H}^{2})Y^{p-1}-2dp(3H^{3}-12H\dot{H}-4\ddot{H})(p-1)Y^{p-2}\dot{Y}$
$\displaystyle+2dp(2\dot{H}+3H^{2})\big{(}(p-1)(p-2)Y^{p-3}\dot{Y}^{2}+(p-1)Y^{p-2}\ddot{Y}\big{)}.$
(20)

As an example, let us consider the power-law scale factor and the
corresponding Hubble parameter:

$\displaystyle a(t)=a_{0}t^{n},\;\;\;H(t)=\frac{n}{t}.$ (21)

The solutions of this type are very useful in cosmology, especially in
describing various phases in the history of the Universe [59]. We note that
for $0<n<1$, the decelerated universe occurs where for dust dominated Universe
$n=\frac{2}{3}$ and for the radiation era $n=\frac{1}{2}$. The corresponding
Ricci scalar and Ricci squared are:

$\displaystyle
R=-6\frac{n}{t^{2}}(1-2n),\;\;\;\;Y=12n^{2}\frac{(3n^{2}-3n+1)}{t^{4}}.$ (22)

By using this expressions in the field equations, one can reconstruct the
Lagrange multiplier and the mimetic potential for the
$f(R,R_{\mu\nu}R^{\mu\nu}$ gravity. The general expressions obtained for the
Lagrange multiplier and potential are given in the Appendix. For the dust
$n=2/3$ and the inverse $p=-1$ power of $Y$, one gets:

$\displaystyle\lambda(t)=d\frac{45}{8}t^{4}-\frac{4}{3t^{2}}$ (23)

and

$\displaystyle V(t)=d\frac{135}{16}t^{4}-\frac{8}{3t^{2}}.$ (24)

It is worth to mention, that the Ricci squared contribution leads to the
higher-order terms in the expressions for $V$ and $\lambda$. In the absence of
this term ($d=0$), the potential takes form $\frac{8}{3}t^{-2}$ that is
consistent with the result obtained in the standard mimetic gravity, i.e.
$V(t)=C_{0}t^{-2}$ for constant $C_{0}=8/3$ [19].

As an another example we consider the exponential symmetric bounce cosmology,
associated with the scale factor and the Hubble rate:

$\displaystyle a(t)=e^{\alpha t^{2}},\;\;\;\;H(t)=2\alpha t.$ (25)

One of the main characteristics of the bounce cosmologies is the absence of
the initial singularity. These theories constitute also interesting
alternatives to the well-established inflationary paradigm [77, 62]. For these
scenarios, the Universe contraction decreases the effective radius of the
Universe to the minimal size. Then, accelerated expansion occurs. In this case
the Ricci scalar and Ricci tensor squared are:

$\displaystyle
R=-6(2\alpha+8t^{2}\alpha^{2}),\;\;\;\;Y=48\alpha^{2}(1+6t^{2}\alpha+12t^{4}\alpha^{2}).$
(26)

The obtained Lagrange multiplier and mimetic potential take form:

$\displaystyle\lambda(t)=4\alpha-\frac{d(3456\alpha^{5}t^{10}+11376\alpha^{4}t^{8}+8208\alpha^{3}t^{6}+1788\alpha^{2}t^{4}+48\alpha
t^{2}-11)}{144\alpha^{2}(12\alpha^{2}t^{4}+6\alpha t^{2}+1)^{4}},$ (27)
$\displaystyle
V(t)=-6(2\alpha+8t^{2}\alpha^{2})-2(2\alpha+12t^{2}\alpha^{2})+\frac{d(1152\alpha^{6}t^{12}+2976\alpha^{5}t^{10}+816\alpha^{4}t^{8}-472\alpha^{3}t^{6}-124\alpha^{2}t^{4}+14\alpha
t^{2}+3)}{16\alpha^{2}(12\alpha^{2}t^{4}+6\alpha t^{2}+1)^{4}},$ (28)

when again we use $p=-1$ and the general form of the potential and the
Lagrange multiplier is listed in the Appendix.

In the following part, we provide reconstruction of the model that describes
transition between matter-dominated and accelerated phases of the Universe’s
history. In this scenario, the Hubble rate is given by [29]:

$\displaystyle H=g_{0}+\frac{g_{1}}{t}.$ (29)

It is important to remark, that for small $t$, $H\approx\frac{g_{1}}{t}$ and
the Universe is filled with the perfect fluid with the EoS (equation of state)
parameter $w=-1+2/3g_{1}$. On the other hand, for the large $t$, the Hubble
rate approaches $H\rightarrow g_{0}$ and the Universe looks like the de-Sitter
one [10, 29]. We note that in this scenario, there is no real matter and the
contribution of the mimetic $f(R,Y)$ gravity plays the role of the standard
matter content.

In the small $t$ limit, the Lagrange multiplier and the mimetic potential are
equal to:

$\displaystyle\lambda(t)\approxeq\frac{12dg_{1}{}^{2}t^{6}-43dg_{1}t^{6}+30dt^{6}-648g_{1}{}^{8}+1296g_{1}{}^{7}-1080g_{1}{}^{6}+432g_{1}{}^{5}-72g_{1}{}^{4}}{36g_{1}{}^{3}(3g_{1}{}^{2}-3g_{1}+1){}^{2}t^{2}},$
(30) $\displaystyle
V(t)\approxeq\frac{6dg_{1}{}^{3}t^{6}-19dg_{1}{}^{2}t^{6}-15dg_{1}t^{6}+20dt^{6}+648g_{1}{}^{9}-1728g_{1}{}^{8}+1944g_{1}{}^{7}-1152g_{1}{}^{6}+360g_{1}{}^{5}-48g_{1}{}^{4}}{12g_{1}{}^{3}(3g_{1}{}^{2}-3g_{1}+1){}^{2}t^{2}},$
(31)

while in the large $t$ limit one gets:

$\displaystyle\lambda(t)\approxeq 0,\;\;\;V(t)\approxeq
6g_{0}^{2}+\frac{d}{16g_{0}^{4}}$ (32)

where again we have chosen $p=-1$. Note also that the full form of the
obtained potential and the Lagrange multiplier can be found in the Appendix.
Thus, in principle, the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$ gravity is able to
unify phases dominated by matter with the transition to the late time
accelerated evolution even in the absence of real matter for the appropriate
mimetic potential and Lagrange multiplier.

To this end, we reconstruct the particular class of
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity that satisfy mimetic field equations, once
particular mimetic potential and the Lagrange multiplier are assumed. This is
inversion of the previously used reconstruction. In this approach, as an
example, we consider the de-Sitter space as specified by the constant Hubble
parameter:

$\displaystyle H=H_{0}.$ (33)

We remark that this spacetime usually serves as a description of the
accelerated expansion phase and the early Universe inflation [59]. By assuming
the constant potential $V$ and the multiplier $\lambda$ as:

$\displaystyle V(t)=\alpha H_{0},\;\;\;\;\lambda(t)=2\beta H_{0},$ (34)

we obtain the following partial differential equation (PDE):

$\displaystyle
3H_{0}^{2}f_{R}+\frac{1}{2}f+6H_{0}^{2}+\frac{(\beta-\alpha)}{2}H_{0}-18H_{0}^{2}f_{Y}=0.$
(35)

This PDE has the following solution:

$\displaystyle
f(R,Y)=12H_{0}^{2}+(\alpha-\beta)H_{0}+e^{-\frac{R}{6H_{0}^{2}}}F(6H_{0}^{2}R+Y),$
(36)

where $F$ is an arbitrary function.

In conclusion, the reconstruction technique can work in both ways: one can
assume specific functional form of the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$ and
cosmological evolution to obtain the mimetic potential and the Lagrange
multiplier and _vice versa_. In general, by using this procedure one can
obtain any cosmological scenario for the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$
gravity.

##  IV Inflation

Inflation is believed to be one of the fundamental building blocks of the
modern cosmology that may solve some of its problems [58, 65]. This is to say,
the main goal of the inflationary theory is to explain primordial fluctuations
i.e. the density variations occurring in very early stages of the cosmic
evolution. Since inflation is the most plausible scenario for the early
Universe, the modified theories of gravity should be able to properly describe
this phase of the cosmic history. Moreover, the deviations from GR can be
interpreted as a quantum-induced corrections or motivated by the ultra-violet
(UV) behavior of the quantum gravity, playing an important role in the
inflationary phase of the early Universe [78]. Thus, viable description of
this stage in any modified gravity is desired. Therefore, in the following
section we analyse inflation in the context of the mimetic
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity to explore ability to obtain model
comparable with the recent BICEP2/Keck observations [60, 78]. Note, that the
cosmological reconstruction can also be successfully applied to the discussed
scenario. In this section we obtain the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$
gravity that describes inflationary model which is consistent with the recent
observational data. Again, the functional form of Eq. (18) is assumed.

While discussing the inflationary cosmology it is worth to use the number of
e-foldings $N$ (intervals for which the scale factor grows by the $N$ factors
of $e$) instead of the cosmic time $t$. The relation between the scale factor
and the e-folding number $N$ is following, $e^{N}=\frac{a}{a_{0}}$, where
$a_{0}$ is the initial value of the scale factor in the initial time instance.
For brevity, we list transformation rules for the time derivatives with
respect to the e-foldings number:

$\displaystyle\frac{d}{dt}=H(N)\frac{d}{dN},\;\;\;\;\frac{d^{2}}{dt^{2}}=H^{2}(N)\frac{d^{2}}{dN^{2}}+H(N)H^{\prime}(N)\frac{d}{dN},$
$\displaystyle\frac{d^{3}}{dt^{3}}=3H^{2}(N)H^{\prime}(N)\frac{d^{2}}{dN^{2}}+H^{2}(N)H^{\prime\prime}(N)\frac{d}{dN}+H(N)H^{\prime
2}(N)\frac{d}{dN}+H^{3}(N)\frac{d^{3}}{dN^{3}}.$ (37)

We note that prime symbols (′) correspond to the derivatives
$\frac{dH(N)}{dN}$.

In order to study the slow-roll indices we use perfect fluid approach
developed in [79]. In this approach extra terms in the gravitational action
(1) can be regarded as the perfect fluid. This formalism allows to obtain the
spectral indices independently from the model. According to this approach, the
slow-roll parameters are given as a functions of the Hubble rate $H$:

$\displaystyle\epsilon$
$\displaystyle=\frac{-4H(N)}{H^{\prime}(N)}\Big{(}\frac{H^{\prime
2}(N)+6H^{\prime}(N)H(N)+H^{\prime\prime}(N)H(N)}{H^{\prime}(N)H(N)+3H^{2}(N)}\Big{)}^{2},$
$\displaystyle\eta$
$\displaystyle=-\frac{\Big{(}9\frac{H^{\prime}(N)}{H(N)}+3\frac{H^{\prime\prime}(N)}{H(N)}+\frac{1}{2}\Big{(}\frac{H^{\prime}(N)}{H(N)}\Big{)}^{2}-\frac{1}{2}\Big{(}\frac{H^{\prime\prime}(N)}{H(N)}\Big{)}^{2}+\frac{H^{\prime\prime}(N)}{H^{\prime}(N)}+\frac{H^{\prime\prime\prime}(N)}{H^{\prime}(N)}\Big{)}}{2\big{(}3+\frac{H^{\prime}(N)}{H(N)}\big{)}}.$
(38)

By using the above, the spectral index of the curvature perturbations $n_{s}$
and the scalar-to-tensor ratio $r$ in terms of the slow-roll parameters can be
provided:

$\displaystyle n_{s}=-6\epsilon+2\eta+1,\;\;\;\;r=16\epsilon.$ (39)

We remark that the spectral index describes variation of the density
fluctuations with respect to the scale. On the other hand, the scalar-tensor
ratio relates spectras of the scalar and tensor perturbations. By using Eqs.
(16), (17) and (37), the field equations in terms of the e-foldings number $N$
are:

$\displaystyle\lambda(N)=H(N)^{2}f_{R}^{\prime\prime}(N)-H(N)^{3}f_{R}^{\prime}(N)H^{\prime}(N)+H(N)f_{R}^{\prime}(N)H^{\prime}(N)+2f_{R}(N)H(N)H^{\prime}(N)+4H(N)^{3}f_{Y}^{\prime\prime}H^{\prime}(N)$
$\displaystyle+12H(N)^{2}f_{Y}^{\prime}H^{\prime}(N)^{2}+H(N)^{3}f_{Y}^{\prime}\Big{(}8H^{\prime\prime}(N)+30H^{\prime}(N)\Big{)}-6H(N)^{4}\Big{(}f_{Y}^{\prime}(N)-f_{Y}^{\prime\prime}(N)\Big{)}+4f_{Y}(N)H(N)H^{\prime}(N)^{3}$
$\displaystyle+4f_{Y}H(N)^{3}\Big{(}H^{(3)}(N)+3H^{\prime\prime}(N)\Big{)}+16f_{Y}H(N)^{2}H^{\prime}(N)\Big{(}H^{\prime\prime}(N)+H^{\prime}(N)\Big{)},$
(40) $\displaystyle
V(N)=f+2H(N)^{2}f_{R}^{\prime\prime}(N)+2H(N)f_{R}^{\prime}(N)H^{\prime}(N)+4H(N)^{2}f_{R}^{\prime}(N)-2f_{R}(N)H(N)\Big{(}H(N)(H^{\prime\prime}(N)+3)+H^{\prime}(N)^{2}\Big{)}$
$\displaystyle+8H(N)^{3}f_{Y}^{\prime\prime}(N)H^{\prime}(N)+12H(N)^{4}f_{Y}^{\prime\prime}(N)+16H(N)^{3}f_{Y}^{\prime}(N)H^{\prime\prime}(N)+84H(N)^{3}f_{Y}^{\prime}(N)H^{\prime}(N)+24H(N)^{2}f_{Y}^{\prime}(N)$
$\displaystyle\times
H^{\prime}(N)^{2}+24H(N)^{4}f_{Y}^{\prime}(N)+4f_{Y}(N)H(N)\Big{(}2H^{\prime}(N)^{3}+2H(N)^{2}(H^{(3)}(N)+6H^{\prime\prime}(N)\Big{)}+4H(N)H^{\prime}(N)\Big{(}2H^{\prime\prime}(N)$
$\displaystyle+3H^{\prime}(N)\Big{)}-9H(N)^{3}),$ (41)

where we have used relations from Eqs. (37). Thus, one can obtain the Lagrange
multiplier and potential corresponding to the scale factor $H(N)$, effectively
reconstructing any given inflationary model. We note, that inverse
reconstruction method is also possible once $\lambda(N)$ and $V(N)$ are
specified.

As an example of the inflationary reconstruction, we consider a inflationary
model described by the Hubble factor [47]:

$\displaystyle H(N)=\alpha e^{\beta N}+\gamma.$ (42)

The corresponding mimetic potential and the Lagrange multiplier are given by:

$\displaystyle\lambda(N)$ $\displaystyle=\frac{1}{{36(\gamma+\alpha e^{\beta
N})^{3}}}\Bigg{[}\alpha\beta e^{\beta N}(72(\gamma+\alpha e^{\beta
N})^{4}-\frac{d}{{(3\gamma^{2}+\alpha^{2}(\beta^{2}+3\beta+3)e^{2\beta
n}+3\alpha(\beta+2)\gamma e^{\beta n})^{4}}}$
$\displaystyle\times\Big{[}-18(\beta^{2}+3\beta-6)\gamma^{6}+\alpha^{6}(\beta^{2}+3\beta+3)^{2}(30\beta^{2}+43\beta+12)e^{6\beta
n}+\alpha^{5}(30\beta^{6}+394\beta^{5}+1806\beta^{4}+4095\beta^{3}\beta$
$\displaystyle+4941\beta^{2}+2961+648)\gamma e^{5\beta
n}+3\alpha^{4}(3\beta^{6}+69\beta^{5}+494\beta^{4}+1575\beta^{3}+2502\beta^{2}+1920\beta+540)\gamma^{2}e^{4\beta
n}+3\alpha^{3}$
$\displaystyle\times(12\beta^{5}+142\beta^{4}+725\beta^{3}+1698\beta^{2}+1830\beta+720)\gamma^{3}e^{3\beta
n}+3\alpha^{2}(10\beta^{4}+90\beta^{3}+432\beta^{2}+825\beta+540)\gamma^{4}e^{2\beta
n}$ $\displaystyle-9\alpha(2\beta^{3}+3\beta^{2}-37\beta-72)\gamma^{5}e^{\beta
n})\Big{]}\Bigg{]},$ (43) $\displaystyle V(N)$
$\displaystyle=-\frac{1}{{36(\gamma+\alpha e^{\beta
N})^{3}}}\Bigg{[}-\frac{3d(\gamma+\alpha e^{\beta
N})}{3\gamma^{2}+\alpha^{2}(\beta^{2}+3\beta+3)e^{2\beta
N}+3\alpha(\beta+2)\gamma e^{\beta N}}-216\alpha\beta e^{\beta
N}(\gamma+\alpha e^{\beta N})^{4}$
$\displaystyle-\frac{1}{{(3\gamma^{2}+\alpha^{2}(\beta^{2}+3\beta+3)e^{2\beta
N}+3\alpha(\beta+2)\gamma e^{\beta N})^{3}}}\Big{[}4\alpha\beta de^{\beta
N}(3\gamma^{2}+\alpha^{2}(4\beta^{2}+9\beta+3)e^{2\beta N}$
$\displaystyle+\alpha(2\beta^{2}+9\beta+6)\gamma e^{\beta
N})(3(\beta+4)\gamma^{2}+4\alpha^{2}(\beta^{2}+3\beta+3)e^{2\beta
N}+\alpha(2\beta^{2}+15\beta+24)\gamma e^{\beta N})\Big{]}$
$\displaystyle+\frac{d(-9\gamma^{3}+3\alpha^{3}(4\beta^{3}+8\beta^{2}-3)e^{3\beta
N}+3\alpha^{2}(4\beta^{3}+12\beta^{2}-9)\gamma e^{2\beta
N}+\alpha(2\beta^{3}+12\beta^{2}-27)\gamma^{2}e^{\beta
N})}{(3\gamma^{2}+\alpha^{2}(\beta^{2}+3\beta+3)e^{2\beta
N}+3\alpha(\beta+2)\gamma e^{\beta N})^{2}}$
$\displaystyle+\frac{1}{(3\gamma^{2}+\alpha^{2}(\beta^{2}+3\beta+3)e^{2\beta
N}+3\alpha(\beta+2)\gamma e^{\beta N})^{4}}\big{(}2\alpha\beta^{2}de^{\beta
N}(3\gamma+\alpha(2\beta+3)e^{\beta N})(-9(\beta+4)\gamma^{5}$
$\displaystyle+28\alpha^{5}(\beta^{2}+3\beta+3)^{2}e^{5\beta
N}+4\alpha^{4}(7\beta^{4}+72\beta^{3}+255\beta^{2}+396\beta+243)\gamma
e^{4\beta N}$
$\displaystyle+\alpha^{3}(8\beta^{4}+141\beta^{3}+786\beta^{2}+1719\beta+1368)\gamma^{2}e^{3\beta
N}+3\alpha^{2}(7\beta^{3}+64\beta^{2}+231\beta+264)\gamma^{3}e^{2\beta N}$
$\displaystyle+3\alpha(2\beta^{2}+15\beta+36)\gamma^{4}e^{\beta
N})\big{)}+72\alpha^{2}\beta^{2}e^{2\beta N}(\gamma+\alpha e^{\beta
N})^{4}+72\alpha\beta^{2}e^{\beta N}(\gamma+\alpha e^{\beta N})^{5}\Bigg{]}.$
(44)

Further, the slow-roll parameters from the perfect fluid approach are:

$\displaystyle\epsilon=-\frac{\alpha\beta e^{\beta
n}((\beta+6)\gamma+2\alpha(\beta+3)e^{\beta n})^{2}}{4(\gamma+\alpha e^{\beta
n})(3\gamma+\alpha(\beta+3)e^{\beta n})^{2}}$ (45)

and

$\displaystyle\eta=-\frac{\beta((\beta+6)\gamma^{2}+8\alpha^{2}(\beta+3)e^{2\beta
n}+2\alpha(4\beta+15)\gamma e^{\beta n})}{4(\gamma+\alpha e^{\beta
n})(3\gamma+\alpha(\beta+3)e^{\beta n})}.$ (46)

Finally the corresponding spectral index and the scalar-tensor ratios are
listed below:

$\displaystyle n_{s}$ $\displaystyle=\frac{1}{2(\gamma+\alpha e^{\beta
n})(3\gamma+\alpha(\beta+3)e^{\beta
n})^{2}}\Big{[}-3(\beta^{2}+6\beta-6)\gamma^{3}+2\alpha^{3}(\beta+3)^{2}(2\beta+1)e^{3\beta
n}$ $\displaystyle+2\alpha^{2}(2\beta^{3}+16\beta^{2}+39\beta+27)\gamma
e^{2\beta n}+\alpha(2\beta^{3}+3\beta^{2}+12\beta+54)\gamma^{2}e^{\beta
n}\Big{]},$ (47) $\displaystyle r=-\frac{4\alpha\beta e^{\beta
n}((\beta+6)\gamma+2\alpha(\beta+3)e^{\beta n})^{2}}{(\gamma+\alpha e^{\beta
n})(3\gamma+\alpha(\beta+3)e^{\beta n})^{2}}.$ (48)

The inflationary parameters depend on the e-foldings number $N$ and the three
free parameters $\alpha,\beta,\gamma$, which are associated with the Hubble
rate (42). We note that this set can be further simplified by introducing
$\Upsilon=\alpha/\gamma$. In Fig. (1) we present the behavior of the
observational indices as a function of the $\beta$ parameter. For exemplary
purposes, we have chosen three values of ratio $\Upsilon$. Gray regions
indicate parameter range that is compatible with the BICEP2/Keck array data
[60]. Based on the obtained results we observe that the spectral index $n_{s}$
diminishes with increase of the $\beta$, while scalar to tensor ratio $r$
exhibits increasing behavior. This inflationary model was extensively studied
by Zhong and collaborators [47] in mimetic $f(\mathcal{G})$. We refer reader
to their work for the detailed survey of the ranges of parameters that are
compatible with the recent observational data provided by BICEP2/Keck array
[60].

The main goal of this part of our study was to reconstruct proper inflationary
model in the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$ theory. For more detailed
discussion of other inflationary scenarios in the mimetic gravity, we refer
reader to [47]. We conclude that extension of the GR presented here is
comparable with other mimetic gravities and may consitiute interesting
alternative [47, 53].

Figure 1: Behavior of the observational indices as a functions of parameter
$\beta$ for specific ratios $\Upsilon=\alpha/\gamma$ ($N=60$). Gray regions
indicate compatibility of the $n_{s}$ and $r$ with recent Planck and
BICEP2/Keck array data.

##  V Summary

In this work, we presented mimetic extension of the
$f(R,R_{\mu\nu}R^{\mu\nu})$ gravity using the Lagrange multiplier formalism
with the mimetic scalar potential added to the Lagrangian, where mimetic field
isolates conformal degree of freedom. Introduced extension was implemented by
using the well-known reconstruction technique in order to obtain models of
interest. Assuming functional form as $f(R,Y)=R+dY^{p}$, we obtained the
Lagrange multiplier $\lambda$ and the mimetic potential $U(\phi)$ that satisfy
the power law evolution of the Universe. Another reconstructed model described
symmetric bounce cosmology. Moreover, we have reconstructed model that satisfy
transition between matter dominated and accelerated phase in the history of
the Universe, showing that this extension of the mimetic gravity is capable to
unify various stages of cosmic history. Additionally, using the inverse
reconstruction method and choosing appropriate forms of $\lambda$ and $V$ we
have obtained $f(R,Y)$ model describing the de-Sitter model with constant
Hubble factor $H_{0}$. Our work is closed with discussion of the inflationary
cosmology in a given regime. By using the Hubble parameter $H(N)=\alpha
e^{\beta N}+\gamma$ we get Lagrange multiplier and the mimetic potential which
successfully describes inflationary model. The reconstructed inflationary
model is phenomenologically viable for the wide range of the parameters when
confronted with BICEP2/Keck data [60, 47].

In conclusion, the mimetic extension of the $f(R,T,R_{\mu\nu}T^{\mu\nu})$ can
describe any given early or late-time cosmological model in particularly clear
way. In general, the procedure presented here can be extended to other higher-
curvature modifications of the general relativity or theories that include
coupling with the trace of the energy-momentum tensor such as the
$f(R,T,R_{\mu\nu}T^{\mu\nu})$ or $f(\mathcal{G},T)$ approaches [18, 72]. We
also remark that the matter fields were not used in our considerations, since
main feature of the mimetic gravity is to mimic cosmological behavior driven
by the matter fields [19, 54]. In the approach presented here, the mimetic
condition is supported by the higher curvature terms of the
$f(R,R_{\mu\nu}R^{\mu\nu})$ theory, such that the geometry of the considered
theory incarnates matter. Moreover, the mimetic condition can play an
important role in further considerations of the various higher order
extensions of the general relativity and can be extended by inclusion of the
unimodular concept in the considerations [54]. We also wish to point out that
by choosing an appropriate functions and parameters one can obtain analytical
results in the mimetic $f(R,R_{\mu\nu}R^{\mu\nu})$ approach, which generally
is hard task in the context of the modified theories of gravity [5]. Since
there are many other proposals of the modified gravities that satisfy
realistic models with sufficient degree of accuracy, the local constraints
should be employed, such as the Post-Newtonian or Solar-System tests [80].
Future studies should be devoted to this topic, by addressing not only mimetic
extension presented here, but also other proposals such as mimetic versions of
the $f(R,T)$, $f(R)$ or $f(\mathcal{G})$ gravity.

##  Appendix A Appendix: Full form of the Lagrange multipliers and mimetic
potentials

### The power law evolution

$\displaystyle\lambda(t)$
$\displaystyle=\frac{1}{3n(3n^{2}-3n+1)t^{2}}\Bigg{[}d(-2^{2p+1})3^{p}p(8p^{2}-6p+1)t^{2}(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}+dn12^{p}p(24p^{2}-18p+1)t^{2}$
(49)
$\displaystyle\times(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}+6n^{2}(d12^{p}(p-1)pt^{2}(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}-1)-18n^{4}+18n^{3}\Bigg{]}$
$\displaystyle V(t)=$
$\displaystyle\frac{1}{(3n(3n^{2}-3n+1)t^{2})}\Bigg{[}-d3^{p}4^{p+1}p(8p^{2}-6p+1)t^{2}(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}$
(50)
$\displaystyle+3n^{2}(d12^{p}(8p^{2}-8p+3)t^{2}(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}-8)$
$\displaystyle-
dn3^{p+1}4^{p}(16p^{3}-4p^{2}-4p+1)t^{2}(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}+162n^{5}-234n^{4}$
$\displaystyle+9n^{3}(d12^{p}(p-1)t^{2}(\frac{n^{2}(3n^{2}-3n+1)}{t^{4}})^{p}+14)\Bigg{]}$

### The symmetric bounce

$\displaystyle\lambda(t)=2\Big{[}2\alpha+8\alpha^{2}dp(48\alpha^{2}+576\alpha^{4}t^{4}+288\alpha^{3}t^{2})^{p-1}\frac{-6\alpha
d(p-1)pt^{2}(\alpha t^{2}-2)(4\alpha
t^{2}+1)(48\alpha^{2}+576\alpha^{4}t^{4}+288\alpha^{3}t^{2})^{p}}{(12\alpha^{2}t^{4}+6\alpha
t^{2}+1)^{2}}$ (51) $\displaystyle+\frac{d48^{p}(p-1)p(3\alpha
t^{2}+1)(\alpha^{2}(12\alpha^{2}t^{4}+6\alpha
t^{2}+1))^{p}(48\alpha^{3}(4p-5)t^{6}+12\alpha^{2}(8p-9)t^{4}+6\alpha(2p-1)t^{2}+1)}{(12\alpha^{2}t^{4}+6\alpha
t^{2}+1)^{3}})\Big{]},$ $\displaystyle V(t)$
$\displaystyle=d48^{p}(\alpha^{2}(12\alpha^{2}t^{4}+6\alpha
t^{2}+1))^{p}+\alpha^{5}d2^{4p+3}3^{p+1}(p-1)pt^{2}(2\alpha t^{2}+3)(4\alpha
t^{2}+1)(\alpha^{2}(12\alpha^{2}t^{4}+6\alpha t^{2}+1))^{p-2}$
$\displaystyle-576\alpha^{4}dpt^{4}(48\alpha^{2}+576\alpha^{4}t^{4}+288\alpha^{3}t^{2})^{p-1}-12\alpha(4\alpha
t^{2}+1)-4\alpha(6\alpha t^{2}+1)$ (52)
$\displaystyle+\frac{d3^{p}4^{2p+1}(p-1)p(3\alpha
t^{2}+1)(48\alpha^{3}(4p-5)t^{6}+12\alpha^{2}(8p-9)t^{4}+6\alpha(2p-1)t^{2}+1)(\alpha^{2}(12\alpha^{2}t^{4}+6\alpha
t^{2}+1))^{p}}{(12\alpha^{2}t^{4}+6\alpha t^{2}+1)^{3}}.$

### Transition between matter-dominated and accelerated phases

For the small $t$ limit ($g_{0}\rightarrow 0$):

$\displaystyle\lambda(t)$
$\displaystyle=\frac{1}{3g_{1}\left(3g_{1}{}^{2}-3g_{1}+1\right)t^{2}}\Bigg{[}d\left(-2^{2p+1}\right)3^{p}p\left(8p^{2}-6p+1\right)t^{2}\left(\frac{g_{1}{}^{2}\left(3g_{1}{}^{2}-3g_{1}+1\right)}{t^{4}}\right){}^{p}+dg_{1}12^{p}p\left(24p^{2}-18p+1\right)$
(53) $\displaystyle
t^{2}\left(\frac{g_{1}{}^{2}\left(3g_{1}{}^{2}-3g_{1}+1\right)}{t^{4}}\right){}^{p}+6g_{1}{}^{2}\left(d12^{p}(p-1)pt^{2}\left(\frac{g_{1}{}^{2}\left(3g_{1}{}^{2}-3g_{1}+1\right)}{t^{4}}\right){}^{p}-1\right)-18g_{1}{}^{4}+18g_{1}{}^{3}\Bigg{]},$
(54) $\displaystyle V(t)$
$\displaystyle=-\frac{1}{3g_{1}(3g_{1}{}^{2}-3g_{1}+1)t^{2}}\Bigg{[}d3^{p}4^{p+1}p(8p^{2}-6p+1)t^{2}\Big{(}\frac{g_{1}{}^{2}(3g_{1}{}^{2}-3g_{1}+1)}{t^{4}}\Big{)}{}^{p}+3g_{1}{}^{2}\Big{(}d12^{p}(8p^{2}-8p+3)t^{2}$
$\displaystyle\times\Big{(}\frac{g_{1}{}^{2}(3g_{1}{}^{2}-3g_{1}+1)}{t^{4}}\Big{)}{}^{p}-2\Big{)}-dg_{1}3^{p+1}4^{p}(16p^{3}-4p^{2}-4p+1)t^{2}\left(\frac{g_{1}{}^{2}(3g_{1}{}^{2}-3g_{1}+1)}{t^{4}}\right){}^{p}$
$\displaystyle+\Big{(}9g_{1}{}^{3}(d12^{p}(p-1)t^{2}\Big{(}\frac{g_{1}{}^{2}(3g_{1}{}^{2}-3g_{1}+1)}{t^{4}}\Big{)}{}^{p}+4\Big{)}+54g_{1}{}^{5}-72g_{1}{}^{4}\Bigg{]}$
(55)

In the large $t$ limit ($g_{1}\rightarrow 0$):

$\displaystyle\lambda(t)=0,\;\;\;v(t)=6g_{0}{}^{2}-d36^{p}(p-1)\left(g_{0}{}^{4}\right){}^{p}.$
(56)

## References

  * Riess et al. [1998] A. G. Riess, A. V. Filippenko, et al., The Astronomical Journal 116, 1009 (1998).
  * Huterer and Shafer [2017] D. Huterer and D. L. Shafer, Reports on Progress in Physics 81, 016901 (2017).
  * Copeland et al. [2006] E. J. Copeland, M. Sami, and S. Tsujikawa, International Journal of Modern Physics D 15, 1753 (2006).
  * Joyce et al. [2016] A. Joyce, L. Lombriser, and F. Schmidt, Annual Review of Nuclear and Particle Science 66, 95 (2016).
  * Nojiri et al. [2017] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Physics Reports 692, 1 (2017).
  * Julve and Tonin [1978] J. Julve and M. Tonin, Il Nuovo Cimento B 46, 137 (1978).
  * Kiefer [2006] C. Kiefer, Annalen der Physik 15, 129 (2006).
  * Hamber [2009] H. W. Hamber, _Quantum gravitation: The Feynman path integral approach_ (Springer, Berlin, 2009), ISBN 978-3-540-85293-3.
  * Capozziello and De Laurentis [2011] S. Capozziello and M. De Laurentis, Physics Reports 509, 167 (2011).
  * Nojiri and Odintsov [2011] S. Nojiri and S. D. Odintsov, Physics Reports 505, 59 (2011).
  * Clifton et al. [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Physics Reports 513, 1 (2012).
  * Buchdahl [1970] H. A. Buchdahl, Monthly Notices of the Royal Astronomical Society 150, 1 (1970).
  * Sotiriou and Faraoni [2010] T. P. Sotiriou and V. Faraoni, Reviews of Modern Physics 82, 451 (2010).
  * Böhmer et al. [2008] C. G. Böhmer, T. Harko, and F. S. N. Lobo, Astroparticle Physics 29, 386 (2008).
  * Nojiri and Odintsov [2005] S. Nojiri and S. D. Odintsov, Physics Letters B 631, 1 (2005).
  * Elizalde et al. [2010] E. Elizalde, R. Myrzakulov, V. V. Obukhov, and D. Sáez-Gómez, Classical and Quantum Gravity 27, 095007 (2010).
  * de la Cruz-Dombriz and Sáez-Gómez [2012] Á. de la Cruz-Dombriz and D. Sáez-Gómez, Classical and Quantum Gravity 29, 245014 (2012).
  * Sharif and Ikram [2016] M. Sharif and A. Ikram, The European Physical Journal C 76, 640 (2016).
  * Lambiase et al. [2015] G. Lambiase, M. Sakellariadou, A. Stabile, and A. Stabile, Journal of Cosmology and Astroparticle Physics 2015, 003 (2015).
  * Easson [2004] D. A. Easson, International Journal of Modern Physics A 19, 5343 (2004).
  * Carroll et al. [2005] S. M. Carroll, A. De Felice, V. Duvvuri, D. A. Easson, M. Trodden, and M. S. Turner, Physical Review D 71, 063513 (2005).
  * Avramidy and Barvinsky [1985] I. Avramidy and A. Barvinsky, Physics Letters B 159, 269 (1985).
  * Ohta et al. [2018] N. Ohta, R. Percacci, and A. D. Pereira, Physical Review D 97, 104039 (2018).
  * Jordan [1955] P. Jordan, _Schwerkraft und Weltall_ (Friedrich Vieweg und Sohn, Braunschweig, 1955).
  * Brans and Dicke [1961] C. Brans and R. H. Dicke, Physical Review 124, 925 (1961).
  * Haghani et al. [2013] Z. Haghani, T. Harko, F. S. N. Lobo, H. R. Sepangi, and S. Shahidi, Physical Review D 88, 044023 (2013).
  * Odintsov and Sáez-Gómez [2013] S. D. Odintsov and D. Sáez-Gómez, Physics Letters B 725, 437 (2013).
  * Harko et al. [2011] T. Harko, F. S. N. Lobo, S. Nojiri, and S. D. Odintsov, Physical Review D 84, 024020 (2011).
  * Houndjo [2012] M. J. S. Houndjo, International Journal of Modern Physics D 21, 1250003 (2012).
  * Harko and Lobo [2014] T. Harko and F. S. N. Lobo, Galaxies 2, 410 (2014).
  * Sharif and Zubair [2013] M. Sharif and M. Zubair, Journal of Cosmology and Astroparticle Physics 2013, 042 (2013).
  * Kaczmarek and Szczȩśniak [2020] A. Z. Kaczmarek and D. Szczȩśniak, Scientific Reports 10, 18076 (2020).
  * Oikonomou et al. [2007] V. Oikonomou, J. Vergados, and C. Moustakidis, Nuclear Physics B 773, 19 (2007).
  * Shafi et al. [2015] Q. Shafi, Şükrü Hanif Tanyıldızı, and C. S. Ün, Nuclear Physics B 900, 400 (2015).
  * Chamseddine and Mukhanov [2013] A. H. Chamseddine and V. Mukhanov, Journal of High Energy Physics 2013 (2013), ISSN 1029-8479.
  * Dutta et al. [2018] J. Dutta, W. Khyllep, E. N. Saridakis, N. Tamanini, and S. Vagnozzi, Journal of Cosmology and Astroparticle Physics 2018, 041 (2018).
  * Gorji et al. [2020] M. A. Gorji, A. Allahyari, M. Khodadi, and H. Firouzjahi, Physical Review D 101, 124060 (2020).
  * Golovnev [2014] A. Golovnev, Physics Letters B 728, 39 (2014).
  * Malaeb [2015] O. Malaeb, Physical Review D 91, 103526 (2015).
  * Ganz et al. [2019] A. Ganz, P. Karmakar, S. Matarrese, and D. Sorokin, Physical Review D 99, 064009 (2019).
  * Bezerra and Miranda [2019] E. Bezerra and O. D. Miranda, The European Physical Journal C 79, 310 (2019).
  * Chamseddine et al. [2019] A. H. Chamseddine, V. Mukhanov, and T. B. Russ, The European Physical Journal C 79, 558 (2019).
  * Mansoori et al. [2021] S. A. H. Mansoori, A. Talebian, and H. Firouzjahi, Journal of High Energy Physics 2021, 183 (2021).
  * Cárdenas et al. [2021] V. H. Cárdenas, M. Cruz, S. Lepe, and P. Salgado, Physics of the Dark Universe 31, 100775 (2021).
  * Odintsov and Oikonomou [2015] S. Odintsov and V. Oikonomou, Annals of Physics 363, 503 (2015).
  * Oikonomou [2015] V. K. Oikonomou, Physical Review D 92, 124027 (2015).
  * Zhong and Sáez-Chillón Gómez [2018] Y. Zhong and D. Sáez-Chillón Gómez, Symmetry 10, 170 (2018).
  * Leon and Saridakis [2015] G. Leon and E. N. Saridakis, Journal of Cosmology and Astroparticle Physics 2015, 031 (2015).
  * Haghani et al. [2018] Z. Haghani, M. Shiravand, and S. Shahidi, International Journal of Modern Physics D 27, 1850049 (2018).
  * Gorji et al. [2018] M. A. Gorji, S. A. H. Mansoori, and H. Firouzjahi, Journal of Cosmology and Astroparticle Physics 2018, 020 (2018).
  * Baffou et al. [2017] E. H. Baffou, M. J. S. Houndjo, M. Hamani-Daouda, and F. G. Alvarenga, The European Physical Journal C 77, 708 (2017).
  * Bhattacharjee et al. [2020] S. Bhattacharjee, J. R. L. Santos, P. H. R. S. Moraes, and P. K. Sahoo, The European Physical Journal Plus 135, 576 (2020).
  * Odintsov and Oikonomou [2016] S. D. Odintsov and V. K. Oikonomou, Astrophysics and Space Science 361, 236 (2016).
  * Nojiri et al. [2016] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Classical and Quantum Gravity 33, 125017 (2016).
  * Sebastiani et al. [2017] L. Sebastiani, S. Vagnozzi, and R. Myrzakulov, Advances in High Energy Physics 2017, 3156915 (2017).
  * Linde [2007] A. Linde, _Inflationary Cosmology_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007), pp. 1–54, ISBN 978-3-540-74353-8.
  * Lyth and Riotto [1999] D. H. Lyth and A. Riotto, Physics Reports 314, 1 (1999).
  * Guth [1981] A. H. Guth, Physical Review D 23, 347 (1981).
  * Hobson et al. [2006] M. Hobson, G. Efstathiou, and A. Lasenby, _General Relativity: An Introduction for Physicists_ (Cambridge University Press, 2006), ISBN 9781139447546\.
  * Ade et al. [2016] P. A. R. Ade, Z. Ahmed, R. W. Aikin, et al. (Keck Array and BICEP2 Collaborations), Physical Review Letters 116, 031302 (2016).
  * Bars et al. [2014] I. Bars, P. Steinhardt, and N. Turok, Physical Review D 89, 061302 (2014).
  * Brandenberger et al. [1993] R. Brandenberger, V. Mukhanov, and A. Sornborger, Physical Review D 48, 1629 (1993).
  * Novello and Bergliaffa [2008] M. Novello and S. P. Bergliaffa, Physics Reports 463, 127 (2008).
  * Cai and Wilson-Ewing [2014] Y.-F. Cai and E. Wilson-Ewing, Journal of Cosmology and Astroparticle Physics 2014, 026 (2014).
  * Starobinsky [1980] A. Starobinsky, Physics Letters B 91, 99 (1980).
  * Cognola et al. [2008] G. Cognola, E. Elizalde, S. Nojiri, S. D. Odintsov, L. Sebastiani, and S. Zerbini, Physical Review D 77, 046009 (2008).
  * Nojiri and Odintsov [2014] S. Nojiri and S. D. Odintsov, Modern Physics Letters A 29, 1450211 (2014).
  * Odintsov and Oikonomou [2018] S. D. Odintsov and V. K. Oikonomou, Physical Review D 98, 044039 (2018).
  * Capozziello et al. [2010] S. Capozziello, J. Matsumoto, S. Nojiri, and S. D. Odintsov, Physics Letters B 693, 198 (2010).
  * Makarenko [2016] A. N. Makarenko, International Journal of Geometric Methods in Modern Physics 13, 1630006 (2016).
  * Nojiri and Odintsov [2006] S. Nojiri and S. D. Odintsov, Physical Review D 74, 086005 (2006).
  * Zubair and Kousar [2016] M. Zubair and F. Kousar, The European Physical Journal C 76, 254 (2016).
  * Bamba et al. [2012] K. Bamba, R. Myrzakulov, S. Nojiri, and S. D. Odintsov, Physical Review D 85, 104036 (2012).
  * Hossienkhani et al. [2014] H. Hossienkhani, A. Najafi, and N. Azimi, Astrophysics and Space Science 353, 311 (2014).
  * Rudra [2016] P. Rudra, Communications in Theoretical Physics 66, 149 (2016).
  * Carloni et al. [2012] S. Carloni, R. Goswami, and P. K. S. Dunsby, Classical and Quantum Gravity 29, 135012 (2012).
  * Brandenberger and Peter [2017] R. Brandenberger and P. Peter, Foundations of Physics 47, 797 (2017).
  * Bamba and Odintsov [2015] K. Bamba and S. D. Odintsov, Symmetry 7, 220 (2015).
  * Bamba et al. [2014] K. Bamba, S. Nojiri, S. D. Odintsov, and D. Sáez-Gómez, Physical Review D 90, 124061 (2014).
  * Vagnozzi [2017] S. Vagnozzi, Classical and Quantum Gravity 34, 185006 (2017).


# Anomaly Detection with Selective Dictionary Learning

Denis C. ILIE-ABLACHIM
Faculty of Automatic Control and Computers
University Politehnica of Bucharest
<EMAIL_ADDRESS>Bogdan DUMITRESCU
This work was supported by a grant of the Ministry of Research, Innovation and
Digitization, CNCS - UEFISCDI, project number PN-III-P4-PCE-2021-0154, within
PNCDI III. It has also been funded by the Ministry of Investments and European
Projects through the Human Capital Sectoral Operational Program 2014-2020,
Contract no. 62461/03.06.2022, SMIS code 153735. Faculty of Automatic Control
and Computers
University Politehnica of Bucharest
<EMAIL_ADDRESS>
###### Abstract

In this paper we present new methods of anomaly detection based on Dictionary
Learning (DL) and Kernel Dictionary Learning (KDL). The main contribution
consists in the adaption of known DL and KDL algorithms in the form of
unsupervised methods, used for outlier detection. We propose a reduced kernel
version (RKDL), which is useful for problems with large data sets, due to the
large kernel matrix. We also improve the DL and RKDL methods by the use of a
random selection of signals, which aims to eliminate the outliers from the
training procedure. All our algorithms are introduced in an anomaly detection
toolbox and are compared to standard benchmark results.

##  I Introduction

Dictionary Learning (DL) is a representation learning method which aims to
find a sparse representation for a set of signals $\bm{Y}$, represented as a
matrix with $N$ columns (signals) of size $m$. The representation is achieved
by computing a dictionary $\bm{D}$ of size $m\times n$ and a sparse
representation $\bm{X}$ of size $n\times N$ such that a good approximation
$\bm{Y}\approx\bm{DX}$ is obtained. Most applications with dictionary learning
are in problems with image denoising, inpainting, signal reconstruction,
clustering or classification.

In this paper we present novel methods for unsupervised learning, in
particular outlier detection, using DL. The main idea is based on finding a
suited dictionary, which is capable of well representing most signals in a
dataset, while the outlier signals representation should obtain large errors.
Considering the number of outliers significantly lower than the rest of the
signals, we expect the dictionary optimization to generally follow the
directions of the normal signals. Our developments cover both the standard and
the nonlinear (kernel) DL.

Anomaly detection (outlier detection) is the identification of a subset of
signals that have a different representation in relation to the rest of the
data. There are several successful anomaly detection methods, such as
Isolation Forest (IForest) [1], Minimum Covariance Determinant (MCD) [2, 3],
One-class SVM detector (OCSVM) or Principal Component Analysis (PCA) Outlier
Detector [4].

There are also several successful sparse coding algorithms used for anomaly
detection. An idea was presented in [5, 6]. These methods consider the data
representation as a joint sparse linear combination of training data. By
following this technique, the authors try to achieve a direct correlation
between all the available signals. Naturally, non-correlated signals are
considered as being anomalies. Another example is given in [7], where the
anomalies are identified in terms of deviation from a trained model. This
method tries to achieve good sparse representation for jointly distributed
signals, while the other independent signals should be isolated. An overview
of DL can be found in [8].

The paper is organized as follows. Section II introduces a natural way of
solving outlier detection problems using DL algorithms. Section III formulates
a new DL algorithm, called Selective Dictionary Learning, which aims to
improve the anomaly detection algorithm by randomly selecting signals for the
training procedure in order to discourage dictionary adaptation to outliers.
In Section IV we present a reduced kernel version of the DL problem and its
Selective form. Section V contains the experimental results, obtained by
running tests on multivariate data and comparing the results with those of
methods available in a Python toolkit for outlier detection.

##  II Anomaly Detection via Dictionary Learning

The DL problem is formulated as following

$\begin{array}[]{ll}\displaystyle\min_{\bm{D},\bm{X}}&\|\bm{Y}-\bm{D}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{x}_{\ell}\right\|_{0}\leq s,\ell=1:N\\\
&\left\|\bm{d}_{j}\right\|=1,j=1:n,\end{array}$ (1)

where $\left\|\cdot\right\|_{0}$ represents the $0$-pseudo-norm and $s$ is the
sparsity level.

The standard dictionary learning problem can be solved by using simple
strategies. In order to overcome the nonconvexity and the huge dimension of
the problem, the optimization procedure is organized in two steps. This method
is also known as DL by Alternate Optimization. In this way, the problem is
divided in two subproblems: sparse coding and dictionary update. By
alternating these two stages for a given number of iterations, the method can
obtain good local solutions. An iteration consists of computing the sparse
representation $\bm{X}$, while the dictionary $\bm{D}$ is fixed, and then
successively updating the dictionary columns, named atoms, while the sparse
representation is fixed. For sparse coding we use Orthogonal Matching Pursuit
(OMP) [9]. For the dictionary update we use the Approximate version of the
K-SVD algorithm (AK-SVD) [10, 11], which optimizes the atoms and their
representations successively.

A simple strategy for anomaly detection is to compute the representation error

$\bm{E}=\bm{Y}-\bm{D}\bm{X}$ (2)

and identify the signals that obtain bad representations. The score of signal
$i$ is simply the norm $\|\bm{e}_{i}\|$ of the $i$-th column of $\bm{E}$. The
largest the error norm, the more likely that signal is an outlier. The
underlying assumption is that signals that are alike can be better represented
by the dictionary designed when solving (1). However, the dictionary size $n$
and the sparsity level $s$ must be taken smaller than usual, otherwise the
representation may be uniformly good for all signals and even outliers can be
well represented. A small dictionary favors good representations for signals
that are similar, tuning the atoms for this purpose; a bad representation of
the outliers has little effect on the objective of (1), since they are few.
This trade-off is naturally obtained during the optimization.

Of course, since sparse representation is linear, similarity and dissimilarity
can be thought in terms of direction. Normal signals belong to a small number
of low dimensional subspaces and the outliers lie on very different subspaces.
This is a model that is appropriate for some anomaly detection problems but
not suited for others.

##  III Selective Dictionary Learning

In the standard DL algorithm, during the training procedure, both stages could
be affected by the presence of outliers in the training dataset. The problem
of anomaly detection can be solved more easily if we could train the
dictionary only on normal data. By neglecting the outliers from the training
dataset, we expect to obtain higher representation errors for anomalies.
However, this is not possible, since we do not know which signals are normal
and which are outliers.

To describe our strategy for eliminating most of the outliers from the
training process, we introduce two new parameters in the DL algorithm:
$train\\_perc$ (training percent) and $train\\_drop\\_prec$ (training dropout
percent). The first one represents the percent of data that are used during
the sparse coding stage. At each iteration, we first apply a random sampling
on the training data and only $train\\_perc$% of the signals are used for
sparse coding. In the dictionary update stage, we further drop off
$train\\_drop\\_perc$% of the signals, namely those having the worst
representations (largest representation errors). Although the first random
selection can eliminate both normal signals and outliers from a training
iteration, the representation of normal signals is less likely to suffer,
since there are still signals in the current training set that are similar to
them. On the contrary, outliers are more likely to lack good proxies and so
their representation will worsen. The second random selection, that of signals
with bad representation, aims to directly remove outliers from the training
process. The dictionary will be updated to better represent the signals that
already have good representations. Hence, again, the outliers representation
will worsen, but the representation of normal signals not present in the
current selection will not be significantly altered.

The DL problem can be formulated by the use of a zero extended permutation
matrix $\bm{P}$ that is modified at each stage and has the role of randomly
selecting the signals:

$\begin{array}[]{ll}\displaystyle\min_{\bm{D},\bm{X}}&\|\bm{Y}\bm{P}-\bm{D}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{x}_{\ell}\right\|_{0}\leq s,\ell=1:N\\\
&\left\|\bm{d}_{j}\right\|=1,j=1:n.\end{array}$ (3)

##  IV Reduced Kernel Dictionary Learning

Linear spaces can usually hinder good representations. In order to overcome
this problem, the standard DL can easily be extended to a nonlinear space.
This method is called Kernel Dictionary Learning (KDL) and was introduced in
[12] and [13]. By this, we reproject each signal $\bm{y}$ to a nonlinear space
$\phi(\bm{y})$, where $\phi(\cdot)$ is a nonlinear function. The dictionary
$\bm{D}$ is also extended to $\phi(\bm{Y})\bm{A}$, where $\bm{A}$ is a matrix
with unknown coefficients, taking the role of dictionary. The KDL problem is
formulated as

$\begin{array}[]{ll}\displaystyle\min_{\bm{A},\bm{X}}&\|\varphi(\bm{Y})-\varphi(\bm{Y})\bm{A}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{x}_{\ell}\right\|_{0}\leq s,\ell=1:N\\\
&\left\|\varphi(\bm{Y})\bm{a}_{j}\right\|=1,j=1:n.\end{array}$ (4)

The KDL problem can be solved similarly to the DL problem (1) if Mercer
kernels are used, which allows the substitution of a scalar product of feature
vectors $\varphi(\bm{x})^{\top}\varphi(\bm{y})$ with a kernel function
$k(\bm{x},\bm{y})$. However, the problem becomes difficult when using large
datasets, due to the large kernel matrix
$\varphi(\bm{Y})^{\top}\varphi(\bm{Y})$ that results. The size of the kernel
matrix scales linearly with the volume of the data, which leads to a large
volume of memory. Thus this strategy might not be tractable for problems with
large datasets.

In order to overcome this limitation we extend the dictionary $\bm{D}$ to a
smaller nonlinear space by $\varphi(\bm{\bar{Y}})\bm{A}$, where $\bm{\bar{Y}}$
represents a small batch of signals from the original dataset. Permuting the
signals such that $\bm{Y}=[\bm{\bar{Y}}\ \bm{\tilde{Y}}]$, we can write

$\varphi(\bm{\bar{Y}})=\underbrace{[\varphi(\bm{\bar{Y}})\quad\varphi(\bm{\tilde{Y}})}_{\bm{\varphi(Y)}}]\underbrace{\left[\begin{array}[]{l}\bm{I}\\\
\bm{0}\end{array}\right]}_{\bm{P}}.$ (5)

The KDL problem becomes

$\begin{array}[]{ll}\displaystyle\min_{\bm{A},\bm{X}}&\|\varphi(\bm{Y})-\varphi(\bm{\bar{Y}})\bm{A}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{x}_{\ell}\right\|_{0}\leq s,\ell=1:N\\\
&\left\|\varphi(\bm{\bar{Y}})\bm{a}_{j}\right\|=1,j=1:n.\end{array}$ (6)

From (5) and (6) we obtain a new optimization problem

$\begin{array}[]{ll}\displaystyle\min_{\bm{A},\bm{X}}&\|\varphi(\bm{Y})(\bm{I}-\bm{PAX})\|_{F}^{2}.\end{array}$
(7)

We denote

$\bm{E}=\bm{I}-\bm{PAX}$ (8)

the representation error and

$\bm{F}=\left[\bm{I}-\bm{P}\sum_{i\neq
j}\bm{a}_{i}\bm{x}_{i}^{T}\right]_{I_{j}}$ (9)

the representation error without the contribution of the current atom
$\bm{a}_{j}$; by $I_{j}$ we denote the set of signal indices to whose
representation $\bm{a}_{j}$ contributes. In order to solve the optimization
problem (6), we update the current atom while the other atoms and the
representation are fixed. Removing the index $j$ for a lighter notation, the
atom update problem becomes

$\begin{array}[]{ll}\displaystyle\min_{\bm{a}}&\left\|\varphi(\bm{Y})\left(\bm{F}-\bm{Pax}^{\top}\right)\right\|_{F}^{2}.\end{array}$
(10)

Using the trace form of the squared Frobenius norm, the objective function
becomes

$\begin{array}[]{l}\operatorname{Tr}\left[\left(\bm{F}^{\top}-\bm{xa}^{\top}\bm{P}^{\top}\right)\varphi^{\top}(\bm{Y})\varphi(\bm{Y})\left(\bm{F}-\bm{Pax}^{\top}\right)\right]=\\\
=\operatorname{Tr}\left[\bm{F}^{\top}\bm{KF}\right]-2\bm{x}^{\top}\bm{F}^{\top}\bm{KPa}+\|\bm{x}\|^{2}\bm{a}^{\top}\bm{P}^{\top}\bm{KPa}.\end{array}$
(11)

We compute the partial derivative of the objective function with respect to
the current atom

$\frac{\partial(\cdot)}{\partial\bm{a}}=2\|\bm{x}\|^{2}\underbrace{\bm{P}^{\top}\bm{KP}}_{\bar{\bm{K}}}\bm{a}-2\underbrace{\bm{P}^{\top}\bm{K}}_{\hat{\bm{K}}^{\top}}\bm{Fx}$
(12)

and so the optimal atom is

$\bm{a}=\left(\|\bm{x}\|^{2}\bar{\bm{K}}\right)^{-1}\hat{\bm{K}}^{\top}\bm{Fx}.$
(13)

The atom is normalized after each update; note that the normalizing factor is
$\left(\bm{a}^{\top}\bar{\bm{K}}\bm{a}\right)^{\frac{1}{2}}$ in order to
obtain $\left\|\bm{a}_{j}\right\|=1$, as required by the original DL problem.

We call Reduced Kernel Dictionary Learning using a Sampled kernel (RKDL-S) the
method solving problem (6) and summarize its update step in Algorithm 1. The
optimal representation from step 6 is computed by setting to zero the partial
derivative of (11) with respect to $\bm{x}$. The sparse representation step,
not listed here, is made using the Kernel OMP algorithm [13].

Data: reduced kernel matrix $\bar{\bm{K}}\in\mathbb{R}^{p\times p}$

partial kernel matrix $\hat{\bm{K}}\in\mathbb{R}^{N\times p}$

current dictionary $\bm{A}\in\mathbb{R}^{N\times n}$

representation matrix $\bm{X}\in\mathbb{R}^{n\times N}$

Result: updated dictionary $\bm{A}$

1 Compute error $\bm{E}=\bm{I}-\bm{P}\bm{A}\bm{X}$

2 for  _$j=1$ to $n$_ do

3 Modify error:
$\bm{F}=\bm{E}_{\mathcal{I}_{j}}+\bm{P}\bm{a}_{j}\bm{X}_{j,\mathcal{I}_{j}}$

4 Update atom:
$\bm{a}_{j}=\left(\|\bm{x}\|_{2}^{2}\bar{\bm{K}}\right)^{-1}\hat{\bm{K}}^{\top}\bm{F}\bm{X}_{j,\mathcal{I}_{j}}$

5 Normalize atom:
$\bm{a}_{j}\leftarrow\bm{a}_{j}/\left(\bm{a}_{j}^{\top}\bm{\bar{K}}\bm{a}_{j}\right)^{\frac{1}{2}}$

6 Update representation:
$\bm{X}_{j,\mathcal{I}_{j}}^{\top}\leftarrow\bm{F}^{\top}\bm{\hat{K}}\bm{a}_{j}$

7 Recompute error:
$\bm{E}_{\mathcal{I}_{j}}=\bm{F}-\bm{P}\bm{a}_{j}\bm{X}_{j,\mathcal{I}_{j}}$

8

Algorithm 1 RKDL-S

RKDL-S achieves good results, but nevertheless in the training process there
are chances to use abnormal signals, by the use of random sampling extraction.
This fact can lead to a decrease in accuracy and performance. A better
strategy that could overcome this problem would be to use a trained dictionary
instead of $\bm{\bar{Y}}$ signals. This can be achieved by using a dictionary,
denoted $\bm{\bar{D}}$, obtained from the linear cases in the previous
sections. The corresponding optimization problem is

$\begin{array}[]{ll}\displaystyle\min_{\bm{A},\bm{X}}&\|\varphi(\bm{Y})-\varphi(\bm{\bar{D}})\bm{A}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{x}_{\ell}\right\|_{0}\leq s,\ell=1:N\\\
&\left\|\varphi(\bm{\bar{D}})\bm{a}_{j}\right\|=1,j=1:n.\end{array}$ (14)

We name it RKDL-D, the last letter indicating the use of dictionary instead of
sampled signals. In order to update the current atom, we rewrite the new
optimization problem as follows

$\min_{\bm{a}_{j}}\left\|\varphi(\bm{Y})-\varphi(\bar{\bm{D}})\sum_{i\neq
j}\bm{a}_{i}\bm{x}_{i}^{\top}-\varphi(\bar{\bm{D}})\bm{a}_{j}\bm{x}_{j}^{\top}\right\|_{F}^{2}.$
(15)

Expressing the Frobenius norm via its trace form, (15) becomes

$\begin{array}[]{r}\displaystyle\min_{\bm{a}_{j}}\operatorname{Tr}\left[\left(\varphi^{\top}(\bm{Y})-\sum_{i\neq
j}\bm{x}_{i}\bm{a}_{i}^{\top}\varphi^{\top}(\bar{\bm{D}})-\bm{x}_{j}\bm{a}_{j}^{\top}\varphi^{\top}(\bar{\bm{D}})\right)\right.\\\
\displaystyle\left.\left(\varphi(\bm{Y})-\varphi(\bar{\bm{D}})\sum_{i\neq
j}\bm{a}_{i}\bm{x}_{i}^{\top}-\varphi(\bar{\bm{D}})\bm{a}_{j}\bm{x}_{j}^{\top}\right)\right].\end{array}$
(16)

After the substitution of scalar products with the kernel function and
negleting the terms that do not depend on $\bm{a}_{j}$ the final optimization
problem is

$\begin{array}[]{r}\displaystyle\min_{\bm{a}_{j}}\operatorname{Tr}\left[2\sum_{i\neq
j}\bm{x}_{i}\bm{a}_{i}^{\top}K(\bm{\bar{\bm{D}}},\bm{\bar{\bm{D}}})\bm{a}_{j}\bm{x}_{j}^{\top}+\bm{x}_{j}\bm{a}_{j}^{\top}\underbrace{K(\bar{\bm{D}},\bar{\bm{D}}}_{\bar{K}_{\bar{\bm{D}}}})\bm{a}_{j}\bm{x}_{j}^{\top}\right.\\\
\left.-2\underbrace{K(\bm{Y},\bar{\bm{D}})}_{\hat{K}_{\bar{\bm{D}}}}\bm{a}_{j}\bm{x}_{j}^{\top}\right].\end{array}$
(17)

Algorithm 1 can be easily modified for solving (17), following the same line
of reasoning as above. In particular, the atom update relation is

$\bm{a}_{j}=\left(\|\bm{x}\|_{2}^{2}\bar{\bm{K}}_{\bar{\bm{D}}}\right)^{-1}(\hat{\bm{K}}^{\top}_{\bar{\bm{D}}}+\bar{\bm{K}}_{\bar{\bm{D}}}R)\bm{X}_{j}$

and the representation update is

$\bm{X}_{j}^{\top}\leftarrow(\bm{\hat{K}}_{\bar{\bm{D}}}-R\bar{\bm{K}}_{\bar{\bm{D}}})\bm{a}_{j},$

where we denoted $\bar{\bm{K}}_{D}$ the reduced kernel matrix
$k(\bar{\bm{D}},\bar{\bm{D}})$, $\hat{\bm{K}}_{D}$ the partial kernel matrix
$k(\bm{Y},\bar{\bm{D}})$ and
$\bm{R}=\bm{X}^{\top}\bm{A}^{\top}-\bm{X}_{j}\bm{a}_{j}^{\top}$ the
transposition representation product with respect to the current atom
$\bm{a}_{j}$. The new method is summarized in Algorithm 2.

Data: reduced kernel matrix $\bar{\bm{K}}_{\bar{\bm{D}}}\in\mathbb{R}^{p\times
p}$

partial kernel matrix $\hat{\bm{K}}_{\bar{\bm{D}}}\in\mathbb{R}^{N\times p}$

current dictionary $\bm{A}\in\mathbb{R}^{N\times n}$

representation matrix $\bm{X}\in\mathbb{R}^{n\times N}$

Result: updated dictionary $\bm{A}$

1 Compute sum
$\displaystyle\bm{S}=\sum_{i=1}^{n}\bm{X}_{i}^{\top}\bm{a}_{i}^{\top}$

2 for  _$j=1$ to $n$_ do

3 Modify sum: $\bm{R}=\bm{S}-\bm{X}_{j}\bm{a}_{j}^{\top}$

4 Update atom:
$\bm{a}_{j}=\left(\|\bm{x}\|_{2}^{2}\bar{\bm{K}}_{\bar{\bm{D}}}\right)^{-1}(\hat{\bm{K}}_{\bar{\bm{D}}}^{\top}+\bar{\bm{K}}_{\bar{\bm{D}}}R)\bm{X}_{j}$

5 Normalize atom:
$\bm{a}_{j}\leftarrow\bm{a}_{j}/\left(\bm{a}_{j}^{\top}\bm{\bar{K}}_{\bar{\bm{D}}}\bm{a}_{j}\right)^{\frac{1}{2}}$

6 Update representation:
$\bm{X}_{j}^{\top}\leftarrow(\bm{\hat{K}}_{\bar{\bm{D}}}-R\bar{\bm{K}}_{\bar{\bm{D}}})\bm{a}_{j}$

7 Recompute error: $\bm{S}=\bm{R}+\bm{X}_{j}\bm{a}_{j}^{\top}$

8

Algorithm 2 RKDL-D

Following the same strategy presented in Section III, the RKDL methods can
easily be adapted to their Selective form. The Selective Reduced Kernel
Dictionary Learning (SRKDL) problem is solved as the previous one, by
introducing two additional steps for the randomly selection of signals, one
for the kernel OMP subproblem and the second one for the matrix coefficients
update subproblem. In both cases the random sampling selection is made
according to the entire data set (including the abnormal signals).

##  V Experiments

In this section we present the main results obtained with the proposed DL
algorithms for anomaly detection. All algorithms have been developed in Python
and have been introduced in the framework of the PyOD [14] anomaly detection
toolbox. For the evaluation, all vectors of a dataset were normalized and were
split into two sets: 60% for training and 40% for testing. Each experiment was
repeated ten times independently with random splits. In terms of performance,
we measure and compute the mean of the area under the receiver operating
characteristic (ROC) curve and the precision @ rank n score. We used 16 real-
world datasets from different domains, more precisely those gathered in ODDS
(Outlier Detection DataSets)111http://odds.cs.stonybrook.edu/ and used as
benchmark in PyOD, and 2 synthetic datasets.

All the algorithms were implemented in Python on a Desktop PC with Ubuntu
20.04 as operating system, having a processor of base frequency of 2.90 GHz
(Max Turbo Frequency 4.80 GHz) and 80GB RAM memory (although a 16 GB RAM
memory is sufficient). During the experiments, ten different rounds were run.
The execution time, receiver operating characteristic value and precision n
score were measured based on the average of all rounds. For the nonlinear
versions we used two different kernels: radial basis function kernel
$k(\bm{x},\bm{y})=\exp{(-\gamma||\bm{x}-\bm{y}||_{2}^{2})}$ and polynomial
kernel $k(\bm{x},\bm{y})=(\gamma\bm{x}^{\top}\bm{y}+\alpha)^{\beta}$. The
hyperparameters of the kernel functions were chosen according to a grid
search. Based on the average results on all the datasets, they were set as
following: $\gamma=1/m$, $\alpha=1$ and $\beta=3$, for the synthetic datasets,
while for the rest we used $\gamma=0.1/m$ for the rbf kernel and $\gamma=10/m$
for the polynomial kernel; we remind that $m$ is the size of a signal. All the
implementations are available at https://github.com/denisilie94/pyod-dl,
including the two synthetic datasets.

The first synthetic dataset was generated based on two different sparse coded
sets of signals. Using two dictionaries, $D_{i}$, the dictionary for inliers,
and $D_{o}$, the dictionary for outliers, two sets of signals were generated
having the sparsity constraint $s=4$. For the numerical experiment we set the
number of inliers $N_{i}=512$ and number of outliers $N_{o}=64$, while the
dictionary size are $n_{i}=50$ and $n_{o}=400$. The signals size was set to
$m=64$. For the outliers signals we used an overcomplete dictionary, since its
representation ability is much more diverse than in the case of the dictionary
for inliers. The second dataset consists of random samples from two normal
(Gaussian) distributions, of different mean and standard deviation. We kept
the same number of normal and abnormal signals of size 64 as in the previous
dataset. The two Gaussian distributions were generated so that the
distribution of normal signals clearly overlaps with the distribution of
abnormal signals. More exactly the inlier mean and variance are $\mu_{i}=0$
and $\sigma_{i}=0.5$, while the outliers parameters are $\mu_{o}=-0.1$ and
$\sigma_{o}=0.45$.

For the DL methods we used small dictionaries of size $n=50$, while the
sparsity constraint was $s=5$. All the dictionaries were trained using $20$
iterations using the AK-SVD method. For the SDL method the $train\\_perc=0.7$
and $train\\_drop\\_perc=0.4$. For the RKDL method, the size of the matrices
$\bm{\bar{Y}}$ from (6) and $\bm{\bar{D}}$ from (14) was set to $10\%$ of the
number of signals. The selective version of RKDL used the parameters
$train\\_perc=0.8$ and $train\\_drop\\_perc=0.3$.

The results show the good behaviour of our algorithms in detecting outliers
via sparse coding. In terms of performance, the DL methods obtain competitive
results. The main results are summarized in Tables I, III for the public PyOD
methods and in Tables II, IV for the DL methods. In all the tables we
highlight the best three results from both sets of methods (PyOD and DL) taken
together. For the synthetic datasets, we noticed that the PyOD methods do not
obtain good results. The DL methods obtain better classification results for
the dataset generated with sparse coding and the dataset with Gaussian
distribution. For ODDS datasets, the overall results are predominantly better
for PyOD methods. However, there are a few datasets where DL methods stand out
as being better. For example, for the cardio dataset, the DL methods achieve
the third place in top, while for the ionosphere and satellite datasets it
occupies the second and third place. An interesting dataset is vertebral where
the DL methods are the best, occupying all three positions of the top.

Data | Samples | Dim. | Out. Perc. | ABOD | CBLOF | FB | HBOS | IForest | KNN | LOF | MCD | OCSVM | PCA
---|---|---|---|---|---|---|---|---|---|---|---|---|---
dl_out | 576 | 64 | 11.1111 | 0.84496 | 0.52271 | 0.58823 | 0.48706 | 0.50225 | 0.56253 | 0.59419 | 0.8734 | 0.5106 | 0.49162
2gauss_out | 576 | 64 | 11.1111 | 0.00633 | 0 | 0.19077 | 0.47946 | 0.28293 | 0 | 0.26359 | 0.00015 | 0.36793 | 0.54262
arrhythmia | 452 | 274 | 14.6018 | 0.76875 | 0.78382 | 0.77807 | 0.82193 | 0.7996 | 0.7861 | 0.77866 | 0.77897 | 0.78116 | 0.7815
cardio | 1831 | 21 | 9.6122 | 0.56917 | 0.81003 | 0.58673 | 0.8351 | 0.91844 | 0.72363 | 0.57357 | 0.82715 | 0.93484 | 0.95038
glass | 214 | 9 | 4.2056 | 0.79507 | 0.84125 | 0.87261 | 0.73887 | 0.74977 | 0.85076 | 0.8644 | 0.79006 | 0.63236 | 0.6747
ionosphere | 351 | 33 | 35.8974 | 0.92476 | 0.89718 | 0.87304 | 0.56144 | 0.85411 | 0.92674 | 0.87535 | 0.95566 | 0.84192 | 0.7962
letter | 1600 | 32 | 6.25 | 0.87825 | 0.78306 | 0.86605 | 0.59268 | 0.64011 | 0.87656 | 0.85935 | 0.8074 | 0.61182 | 0.5283
lympho | 148 | 18 | 4.0541 | 0.91097 | 0.96731 | 0.97528 | 0.99569 | 0.99288 | 0.9745 | 0.97709 | 0.91125 | 0.97587 | 0.98467
mnist | 7603 | 100 | 9.2069 | 0.78153 | 0.84041 | 0.72046 | 0.57419 | 0.80673 | 0.84813 | 0.71608 | 0.86661 | 0.85289 | 0.85266
musk | 3062 | 166 | 3.1679 | 0.18444 | 1 | 0.52626 | 0.99998 | 0.99984 | 0.79857 | 0.52867 | 0.99997 | 1 | 0.99995
optdigits | 5216 | 64 | 2.8758 | 0.46674 | 0.7692 | 0.44336 | 0.87325 | 0.70608 | 0.37076 | 0.45004 | 0.3979 | 0.49972 | 0.50856
pendigits | 6870 | 16 | 2.2707 | 0.68776 | 0.89307 | 0.45953 | 0.92381 | 0.94964 | 0.74865 | 0.46975 | 0.83439 | 0.93031 | 0.93525
pima | 768 | 8 | 34.8958 | 0.67938 | 0.65781 | 0.62345 | 0.69995 | 0.67798 | 0.70781 | 0.62705 | 0.67528 | 0.6215 | 0.64811
satellite | 6435 | 36 | 31.6395 | 0.57137 | 0.74942 | 0.55717 | 0.75811 | 0.6937 | 0.68364 | 0.55727 | 0.80304 | 0.66224 | 0.59884
satimage-2 | 5803 | 36 | 1.2235 | 0.81896 | 0.99922 | 0.45701 | 0.98042 | 0.99384 | 0.9536 | 0.45774 | 0.99593 | 0.9978 | 0.98218
vertebral | 240 | 6 | 12.5 | 0.42615 | 0.43309 | 0.41658 | 0.32625 | 0.39276 | 0.38166 | 0.40811 | 0.39158 | 0.44308 | 0.40269
vowels | 1456 | 12 | 3.4341 | 0.96059 | 0.92221 | 0.94252 | 0.67267 | 0.75966 | 0.968 | 0.94096 | 0.80761 | 0.78021 | 0.60267
wbc | 378 | 30 | 5.5556 | 0.90473 | 0.92005 | 0.93254 | 0.95163 | 0.93073 | 0.93662 | 0.93488 | 0.92102 | 0.93189 | 0.91587
TABLE I: ROC Performance - PyOD methods Data | DL | SDL | RKDL-S | RKDL-D | SRKDL-S | SRKDL-D
---|---|---|---|---|---|---
|  |  | rbf | poly | rbf | poly | rbf | poly | rbf | poly
dl_out | 0.89666 | 0.85336 | 0.36259 | 0.36424 | 0.39653 | 0.37279 | 0.32207 | 0.34155 | 0.3849 | 0.36028
2gauss_out | 0.91278 | 0.90274 | 0.05688 | 0.02432 | 0.0165 | 0.01193 | 1 | 0.82143 | 0.00569 | 0.0036
arrhythmia | 0.77057 | 0.77194 | 0.68557 | 0.63535 | 0.70723 | 0.76356 | 0.72746 | 0.72333 | 0.72112 | 0.78032
cardio | 0.70023 | 0.72884 | 0.63584 | 0.92797 | 0.60367 | 0.82092 | 0.69322 | 0.93747 | 0.63524 | 0.89624
glass | 0.67484 | 0.64841 | 0.77257 | 0.65649 | 0.80463 | 0.68711 | 0.79033 | 0.68143 | 0.62142 | 0.72519
ionosphere | 0.93401 | 0.93923 | 0.87097 | 0.62494 | 0.88993 | 0.80313 | 0.85409 | 0.60376 | 0.8578 | 0.77142
letter | 0.82366 | 0.81978 | 0.72686 | 0.33016 | 0.74923 | 0.41022 | 0.72863 | 0.3328 | 0.71789 | 0.42632
lympho | 0.91635 | 0.93615 | 0.54727 | 0.711 | 0.7846 | 0.95522 | 0.60362 | 0.82363 | 0.82788 | 0.9102
mnist | 0.81029 | 0.79723 | 0.66476 | 0.80017 | 0.54424 | 0.57414 | 0.6443 | 0.71678 | 0.5917 | 0.57642
musk | 0.88574 | 0.89346 | 0.69204 | 0.75747 | 0.55197 | 0.70492 | 0.77949 | 0.9375 | 0.70581 | 0.82653
optdigits | 0.40227 | 0.4084 | 0.36615 | 0.41572 | 0.40347 | 0.56437 | 0.39185 | 0.43933 | 0.48264 | 0.61547
pendigits | 0.62745 | 0.63382 | 0.8356 | 0.91127 | 0.78032 | 0.88282 | 0.83539 | 0.92798 | 0.84152 | 0.91862
pima | 0.56365 | 0.55959 | 0.60599 | 0.64188 | 0.6131 | 0.65217 | 0.62009 | 0.6357 | 0.63752 | 0.65449
satellite | 0.65351 | 0.64655 | 0.65338 | 0.66762 | 0.64236 | 0.59866 | 0.77197 | 0.68252 | 0.77682 | 0.69671
satimage-2 | 0.59438 | 0.55493 | 0.85817 | 0.97076 | 0.90414 | 0.92801 | 0.98003 | 0.96482 | 0.98159 | 0.9646
vertebral | 0.48265 | 0.46904 | 0.39465 | 0.41046 | 0.48752 | 0.39953 | 0.46184 | 0.40483 | 0.50182 | 0.38861
vowels | 0.77689 | 0.80236 | 0.80882 | 0.519 | 0.78853 | 0.66121 | 0.86525 | 0.52215 | 0.84647 | 0.678
wbc | 0.81705 | 0.84737 | 0.67586 | 0.88776 | 0.71915 | 0.9167 | 0.81133 | 0.89594 | 0.73305 | 0.90848
TABLE II: ROC Performance - DL methods Data | Samples | Dim. | Out. Perc. | ABOD | CBLOF | FB | HBOS | IForest | KNN | LOF | MCD | OCSVM | PCA
---|---|---|---|---|---|---|---|---|---|---|---|---|---
dl_out | 576 | 64 | 11.1111 | 0.47533 | 0.11883 | 0.1958 | 0.10082 | 0.07764 | 0.15684 | 0.19525 | 0.47469 | 0.11109 | 0.0999
2gauss_out | 576 | 64 | 11.1111 | 0 | 0 | 0 | 0.10901 | 0.03221 | 0 | 0.004 | 0 | 0.02137 | 0.10618
arrhythmia | 452 | 274 | 14.6018 | 0.38076 | 0.45385 | 0.42297 | 0.51108 | 0.49992 | 0.44637 | 0.43343 | 0.39952 | 0.4614 | 0.46129
cardio | 1831 | 21 | 9.6122 | 0.23743 | 0.42966 | 0.169 | 0.44761 | 0.49186 | 0.33227 | 0.15409 | 0.42084 | 0.50112 | 0.609
glass | 214 | 9 | 4.2056 | 0.17023 | 0.07262 | 0.14762 | 0 | 0.07262 | 0.07262 | 0.14762 | 0 | 0.17262 | 0.07262
ionosphere | 351 | 33 | 35.8974 | 0.84415 | 0.77489 | 0.70558 | 0.32951 | 0.64743 | 0.86021 | 0.70634 | 0.88065 | 0.70005 | 0.57286
letter | 1600 | 32 | 6.25 | 0.38009 | 0.23969 | 0.36419 | 0.07152 | 0.08828 | 0.33117 | 0.36411 | 0.19327 | 0.15096 | 0.08747
lympho | 148 | 18 | 4.0541 | 0.44834 | 0.75167 | 0.75167 | 0.84667 | 0.87667 | 0.75167 | 0.75167 | 0.56833 | 0.75167 | 0.75167
mnist | 7603 | 100 | 9.2069 | 0.3555 | 0.40231 | 0.32986 | 0.11882 | 0.30346 | 0.42043 | 0.33429 | 0.3462 | 0.39619 | 0.38461
musk | 3062 | 166 | 3.1679 | 0.05075 | 1 | 0.22297 | 0.97832 | 0.98069 | 0.2733 | 0.16955 | 0.98889 | 1 | 0.97994
optdigits | 5216 | 64 | 2.8758 | 0.00602 | 0 | 0.02445 | 0.2194 | 0.0271 | 0 | 0.02335 | 0 | 0 | 0
pendigits | 6870 | 16 | 2.2707 | 0.08125 | 0.23974 | 0.06579 | 0.29793 | 0.35505 | 0.09844 | 0.06529 | 0.08928 | 0.32866 | 0.31865
pima | 768 | 8 | 34.8958 | 0.51929 | 0.48378 | 0.44802 | 0.54238 | 0.50233 | 0.54133 | 0.45552 | 0.49625 | 0.47035 | 0.49429
satellite | 6435 | 36 | 31.6395 | 0.39023 | 0.57978 | 0.39016 | 0.56903 | 0.55766 | 0.49945 | 0.38929 | 0.68452 | 0.53455 | 0.47844
satimage-2 | 5803 | 36 | 1.2235 | 0.21305 | 0.93759 | 0.0555 | 0.6939 | 0.8775 | 0.38087 | 0.05551 | 0.64813 | 0.93556 | 0.80408
vertebral | 240 | 6 | 12.5 | 0.06005 | 0.03381 | 0.06439 | 0.00714 | 0.05337 | 0.02381 | 0.05059 | 0 | 0.02381 | 0.02262
vowels | 1456 | 12 | 3.4341 | 0.57102 | 0.36427 | 0.3224 | 0.12974 | 0.19602 | 0.50929 | 0.35506 | 0.2186 | 0.27907 | 0.13636
wbc | 378 | 30 | 5.5556 | 0.30604 | 0.48064 | 0.51879 | 0.58166 | 0.50879 | 0.49518 | 0.51879 | 0.45771 | 0.51249 | 0.47673
TABLE III: Precision @ N Performances - PyOD methods Data | DL | SDL | RKDL-S | RKDL-D | SRKDL-S | SRKDL-D
---|---|---|---|---|---|---
|  |  | rbf | poly | rbf | poly | rbf | poly | rbf | poly
dl_out | 0.54116 | 0.53812 | 0.01536 | 0.02597 | 0.01429 | 0.02937 | 0.00385 | 0.01629 | 0.01087 | 0.01447
2gauss_out | 0.51047 | 0.49212 | 0.00333 | 0 | 0 | 0 | 1 | 0.78568 | 0 | 0
arrhythmia | 0.42828 | 0.42783 | 0.32057 | 0.31022 | 0.35562 | 0.42643 | 0.38116 | 0.37517 | 0.37499 | 0.45478
cardio | 0.30468 | 0.30867 | 0.22022 | 0.5416 | 0.19782 | 0.35779 | 0.24374 | 0.58545 | 0.1814 | 0.50126
glass | 0.1369 | 0.04762 | 0.12262 | 0.09762 | 0.11429 | 0.07262 | 0.14762 | 0.03929 | 0.125 | 0.09762
ionosphere | 0.8081 | 0.81393 | 0.73023 | 0.43417 | 0.78044 | 0.61622 | 0.70715 | 0.42508 | 0.72037 | 0.58191
letter | 0.27955 | 0.26227 | 0.1753 | 0.01572 | 0.2086 | 0.03195 | 0.15485 | 0.02421 | 0.16185 | 0.02761
lympho | 0.49833 | 0.45666 | 0.125 | 0.12833 | 0.25833 | 0.42833 | 0.22333 | 0.265 | 0.26667 | 0.465
mnist | 0.37567 | 0.36151 | 0.23029 | 0.35623 | 0.13266 | 0.12771 | 0.18196 | 0.28904 | 0.16703 | 0.13247
musk | 0.4075 | 0.37901 | 0.12873 | 0.23353 | 0.10472 | 0.17828 | 0.21692 | 0.68521 | 0.26918 | 0.0882
optdigits | 0.00963 | 0.00696 | 0.0111 | 0 | 0.02303 | 0.04896 | 0.00474 | 0 | 0.0437 | 0.1143
pendigits | 0.09973 | 0.09568 | 0.16584 | 0.28128 | 0.10907 | 0.27757 | 0.1476 | 0.27513 | 0.16232 | 0.25404
pima | 0.41686 | 0.39683 | 0.43692 | 0.4744 | 0.44025 | 0.48196 | 0.45952 | 0.48992 | 0.4757 | 0.49287
satellite | 0.46056 | 0.45341 | 0.47054 | 0.51423 | 0.45615 | 0.4217 | 0.61256 | 0.54732 | 0.62188 | 0.55516
satimage-2 | 0.07119 | 0.04584 | 0.09132 | 0.64813 | 0.33605 | 0.38399 | 0.58768 | 0.64744 | 0.56207 | 0.62129
vertebral | 0.09294 | 0.07273 | 0.06198 | 0.05927 | 0.10508 | 0.02143 | 0.05048 | 0.03214 | 0.07903 | 0.00667
vowels | 0.28151 | 0.30329 | 0.22516 | 0.10201 | 0.21783 | 0.18607 | 0.25828 | 0.06196 | 0.22259 | 0.1181
wbc | 0.41909 | 0.36207 | 0.18477 | 0.55319 | 0.24709 | 0.52414 | 0.27334 | 0.5301 | 0.15626 | 0.51664
TABLE IV: Precision @ N Performances - DL methods

In general, the SDL method achieve better results than the DL method, but this
is not always true. Depending on how the random selection of signals is made,
there are chances that abnormal signals to be used during the training
procedure. This is possible for the datasets with a very high percentage of
outliers or small datasets. The same statement is valid for the KDL vs SKDL
comparison. On the other hand, comparing the standard methods with the kernel
methods, we notice that the second ones obtain better results. Moreover, the
selective strategy improves the invariance of dictionaries to representing
abnormal signals. The RKDL-D and SRKDL-D methods often improve the results. In
general, the trained dictionary, $\bm{D}$, is better adapted for the
representation of the normal signals. However, it is likely that the trained
dictionaries contain atoms that are beneficial in the nonlinear representation
of all signals, including the outliers.

The execution time of DL methods is usually larger than that of the PyOD
methods. For example, for the musk dataset, which is among the largest, DL and
SDL take about 6 seconds, i.e., not much more than MCD, which needs about 4
seconds; RKDL algorithms take between 9 and 11 seconds, while SRKDL variant
are slightly faster, with 7-10 seconds. The other PyOD algorithms are at least
$10$ times faster than the methods presented in the article.

##  VI Conclusions

In this paper we have presented a novel unsupervised method for outlier
detection, based on Dictionary Learning and Kernel Dictionary Learning. We
have introduced a reduced kernel DL version that is suitable for problems with
large datasets. The kernel reduction technique is based on choosing a small
sample of signals from the original dataset, which will further be used for
the nonlinear extension. Another way to represent the kernel is to use a
dictionary initially trained in with the standard DL algorithm. Both methods
are accompanied by improved versions based on a random selection of the data
used in the training procedure. This ensures invariance in the representation
of normal signals, while the capabilities of the dictionaries for the
representation of abnormal signals decrease.

Based on these results, we demonstrated that sparse learning can easily
isolate the outliers from the normal signals, while obtaining competitive
results with other unsupervised methods. All the developed algorithms were
introduced in an outlier detection toolbox.

## References

  * [1] Fei Tony Liu, Kai Ming Ting, and Zhi-Hua Zhou.  Isolation forest.  In 2008 eighth ieee international conference on data mining, pages 413–422. IEEE, 2008.
  * [2] Johanna Hardin and David M Rocke.  Outlier detection in the multiple cluster setting using the minimum covariance determinant estimator.  Computational Statistics & Data Analysis, 44(4):625–638, 2004\.
  * [3] Peter J Rousseeuw and Katrien Van Driessen.  A fast algorithm for the minimum covariance determinant estimator.  Technometrics, 41(3):212–223, 1999.
  * [4] Bernhard Schölkopf, John C Platt, John Shawe-Taylor, Alex J Smola, and Robert C Williamson.  Estimating the support of a high-dimensional distribution.  Neural computation, 13(7):1443–1471, 2001.
  * [5] Nam H Nguyen, Nasser M Nasrabadi, and Trac D Tran.  Robust multi-sensor classification via joint sparse representation.  In 14th International Conference on Information Fusion, pages 1–8. IEEE, 2011.
  * [6] Sumit Shekhar, Vishal M Patel, Nasser M Nasrabadi, and Rama Chellappa.  Joint sparsity-based robust multimodal biometrics recognition.  In European Conference on Computer Vision, pages 365–374. Springer, 2012.
  * [7] Amir Adler, Michael Elad, Yacov Hel-Or, and Ehud Rivlin.  Sparse coding with anomaly detection.  Journal of Signal Processing Systems, 79(2):179–188, 2015.
  * [8] Bogdan Dumitrescu and Paul Irofti.  Dictionary learning algorithms and applications.  Springer, 2018.
  * [9] Yagyensh Chandra Pati, Ramin Rezaiifar, and Perinkulam Sambamurthy Krishnaprasad.  Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition.  In Proceedings of 27th Asilomar conference on signals, systems and computers, pages 40–44. IEEE, 1993.
  * [10] Michal Aharon, Michael Elad, and Alfred Bruckstein.  K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation.  IEEE Transactions on Signal Processing, 54(11):4311–4322, 2006\.
  * [11] Ron Rubinstein, Michael Zibulevsky, and Michael Elad.  Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit.  Technical report, Computer Science Department, Technion, 2008.
  * [12] Jayaraman J Thiagarajan, Karthikeyan Natesan Ramamurthy, and Andreas Spanias.  Multiple kernel sparse representations for supervised and unsupervised learning.  IEEE Transactions on Image Processing, 23(7):2905–2915, 2014.
  * [13] Hien Van Nguyen, Vishal M Patel, Nasser M Nasrabadi, and Rama Chellappa.  Design of non-linear kernel dictionaries for object recognition.  IEEE Transactions on Image Processing, 22(12):5123–5135, 2013.
  * [14] Yue Zhao, Zain Nasrullah, and Zheng Li.  PyOD: A Python Toolbox for Scalable Outlier Detection.  Journal of Machine Learning Research, 20(96):1–7, 2019.