\section{Introduction}

The explanation of the accelerated expansion of the Universe is one of the most challenging problems in modern cosmology \cite{Perlmutter:1998np,Riess:1998cb}. From the mathematical point of view, the simplest way to treat this problem is to consider the cosmological constant term  \cite{Carroll:2000fy}. Nonetheless, this model faces some difficulties such as the coincidence problem and the cosmological constant problem.
The latter dictates a huge discrepancy between the observed values of the vacuum energy density and the theoretical large value of the zero-point energy suggested by quantum field theory \cite{Carroll:2000fy}.
There are some alternative models proposed to overcome these problems such as modified gravity \cite{Nojiri:2010wj,Sotiriou:2008rp,Capozziello:2011et,Avelino:2016lpj,Lobo:2008sg,Bamba:2015uma,Nojiri:2017ncd}, mysterious energy-momentum sources \cite{Copeland:2006wr,Horndeski:1974wa,Deffayet:2011gz}, such as quintessence
\cite{Wetterich:1987fm,Ratra:1987rm,Caldwell:1997ii,Zlatev:1998tr} and $k$-essence \cite{ArmendarizPicon:1999rj,ArmendarizPicon:2000dh,ArmendarizPicon:2000ah} fields, and complex equations of state storing the missing energy of the dark side of the Universe \cite{Bento:2002ps,Arun:2017uaw}.
In general, models with varying dark energy candidates may be capable of overcoming all the mathematical and theoretical difficulties, however, the main underlying question is the origin of these terms.
One proposal is to relate this behavior to the energy of quantum fields in vacuum through the holographic principle which allows us to reconcile infrared (IR) and ultraviolet (UV) cutoffs \cite{Cohen:1998zx} (see Ref.~\cite{Tawfik:2019dda} for a review on various attempts to model the dark side of the Universe).

However, an alternative to these dark energy models, as mentioned above, is modified gravity \cite{Nojiri:2010wj,Sotiriou:2008rp,Capozziello:2011et,Avelino:2016lpj,Lobo:2008sg,Bamba:2015uma,Nojiri:2017ncd}. Here, one considers generalizations of the Hilbert-Einstein Lagrangian specific curvature
invariants such as $R^{2}$, $R_{\mu \nu }R^{\mu \nu }$, $R_{\alpha
\beta \mu \nu }R^{\alpha \beta \mu \nu }$, $\varepsilon ^{\alpha
\beta \mu \nu }R_{\alpha \beta \gamma \delta }R_{\mu \nu }^{\gamma
\delta }$, etc. A popular theory that has attracted much attention is $f(R)$ gravity, where one may tackle the problem through several approaches, namely, the metric formalism \cite{Nojiri:2010wj,Sotiriou:2008rp,Capozziello:2011et,Avelino:2016lpj,Lobo:2008sg,Bamba:2015uma,Nojiri:2017ncd}, which considers that the metric is the fundamental field, or the Palatini formalism \cite{Olmo:2011uz}, where here one varies the action with respect to the metric and an independent connection. However, one may also consider a hybrid combination of these approaches that has recently been proposed, namely, the hybrid metric-Palatini gravitational theory \cite{Harko:2011nh}, where the metric Einstein-Hilbert action is supplemented with a metric-affine (Palatini) correction term. The hybrid metric-Palatini theory has the ability to avoid several of the problematic issues that arise in the pure metric and Palatini formalisms.
For instance, metric $f(R)$ gravity introduces an additional scalar degree of freedom, which must possess a low mass in order to be relevant for the large scale cosmic dynamics. However, the presence of such a low mass scalar field would influence the dynamics on smaller scales, such as at the level of the Solar System, and since these small scale effects remain undetected, one must resort to screening mechanisms \cite{Capozziello:2007eu,Khoury:2003rn}. Relative to the Palatini formalism, no additional degrees of freedom are introduced, as the scalar field is an algebraic function of the trace of the energy-momentum tensor. It has been shown that this fact entails serious consequences for the theory, leading to the presence of infinite tidal forces on the surface of massive astrophysical type objects \cite{Olmo:2011uz}.

In this context, the hybrid metric-Palatini gravity theories were proposed initially in \cite{Harko:2011nh} in order to circumvent the above shorcomings in the metric and Palatini formalisms of $f(R)$ gravity. One of the main advantages of the hybrid metric-Palatini theory is that in its scalar-tensor representation a long-range force is introduced that automatically passes the Solar System tests, and thus no contradiction between the theory and the local measurements arise. In a cosmological context, it has also been shown that hybrid metric-Palatini gravity may also explain the cosmological epochs \cite{Lima:2014aza,Lima:2015nma}
(we refer the reader to \cite{Capozziello:2015lza,Harko:2018ayt,Harko:2020ibn} for more details).
In fact, hybrid metric-Palatini gravity has attracted much attention recently, and has a plethora of applications, namely, in the galactic context \cite{Capozziello:2012qt,Capozziello:2013yha,Capozziello:2013uya}, in cosmology \cite{Capozziello:2012ny,Carloni:2015bua,Capozziello:2013wq,Boehmer:2013oxa,Lima:2014aza,Lima:2015nma}, in braneworlds \cite{Fu:2016szo,Rosa:2020uli}, black holes \cite{Bronnikov:2019ugl,Chen:2020evr,Rosa:2020uoi,Bronnikov:2020vgg}, stellar solutions \cite{Danila:2016lqx}, cosmic strings \cite{Bronnikov:2020zob,Harko:2020oxq}, tests of binary pulsars \cite{Avdeev:2020jqo}, and wormholes \cite{Rosa:2018jwp,Capozziello:2012hr}, amongst others.

Here, we explore the possibility that evolving wormhole geometries may be supported by hybrid metric-Palatini gravity. These compact objects, which are theoretical shortcuts in spacetime, have been shown to be threaded by an exotic fluid that violates the null energy condition (NEC), at least for the static case \cite{Morris:1988cz,Morris:1988tu}. However, it has been shown that evolving wormholes are able to satisfy the
energy conditions in arbitrary but finite intervals of time \cite{Kar:1994tz,Kar:1995ss}, contrary to their static counterparts \cite{Visser:1995cc,Lobo:2017oab}. One way to study this subject is to embed a wormhole in a Friedmann–Lema\^{i}tre–Robertson–Walker (FLRW) metric, which
permits the geometry to evolve in a cosmological background \cite{Roman:1992xj,Anchordoqui:1997du,Arellano:2006ex,Ebrahimi:2009ih,Ebrahimi:2009rv,Bordbar:2010zz,Sajadi:2012zz,Cataldo:2013sma,Setare:2015iqa, Bhattacharya:2015oma,Ovgun:2018uin,Zangeneh:2014noa,Cataldo:2008pm,Cataldo:2008ku,Cataldo:2012pw,Maeda:2009tk,KordZangeneh:2020jio}. Due to the somewhat problematic nature of the energy condition violations, an important issue is whether traversable wormholes can be constructed from normal matter throughout the spacetime, or at least partially. Different wormhole structures have been explored from this point of view \cite{Zangeneh:2014noa,Mehdizadeh:2015jra,Mehdizadeh:2015dta,Zangeneh:2015jda,KordZangeneh:2020jio,Parsaei:2019utg}. Indeed, it has been shown that in modified gravity, it is possible to impose that the matter threading the wormhole throat satisfies the energy conditions \cite{Capozziello:2013vna,Capozziello:2014bqa}, and it is the higher order curvature terms that support these nonstandard wormhole geometries \cite{Lobo:2007qi,Lobo:2009ip,Garcia:2010xb,MontelongoGarcia:2010xd,Harko:2013yb,Korolev:2020ohi}. Furthermore, different wormhole structures have been studied extensively in the context of alternative theories of gravity \cite{Antoniou:2019awm, Tangphati:2019pxh, Papantonopoulos:2019ugr, Godani:2019kgy, Banerjee:2020uyi, Fayyaz:2020jzh, Restuccia:2020wls,Tangphati:2020mir,Singh:2020rai,Ibadov:2020btp,Richarte:2009zz} and also from different aspects \cite{Lazov:2017tjs, Savelova:2019lye, Bak:2019nnu, Xu:2020wfm, Jusufi:2020rpw, Berry:2020tky, Fallows:2020ugr, Maldacena:2020sxe,Lobo:2020kxn,Parsaei:2020hke,DeFalco:2020afv,Moti:2020whf,Wielgus:2020uqz}.

In this paper, we study the evolution of traversable wormholes in a FLRW
universe background in the scalar-tensor representation of hybrid metric-Palatini gravity.
Furthermore, we explore the energy conditions for matter which threads these
wormhole geometries. The paper is organized in the following manner: In Sec. \ref{secII}, we briefly present the action and field equations of hybrid metric-Palatini gravity, we consider the spacetime metric and explore a barotropic equation of state for the background fluid. In Sec. \ref{secIII}, we analyse evolving traversable wormhole geometries for specific values of the barotropic equation of state parameter, as well as evolving wormholes with a traceless energy-momentum tensor (EMT), and study the energy conditions for the solutions obtained. Finally, in Sec. \ref{conclusion}, we summary our results and conclude.

\section{Evolving wormholes in hybrid metric-Palatini gravity} \label{secII}

\subsection{Action and field equations}

Here, we briefly present the hybrid metric-Palatini gravitational theory. The action is given by
\begin{equation}  \label{eq:S_hybrid}
S= \frac{1}{2\kappa^2}\int d^4 x \sqrt{-g} \left[ R + f(\mathcal{R})\right]
+S_m \ ,
\end{equation}
where $\kappa^2\equiv 8\pi G$, $R$ is the metric Ricci scalar, and the Palatini curvature is $\mathcal{R} \equiv g^{\mu\nu}\mathcal{R}_{\mu\nu} $, with the Palatini Ricci tensor, $\mathcal{R}_{\mu\nu}$, defined in terms of an independent connection, $\hat{\Gamma}^\alpha_{\mu\nu}$, given by
\begin{equation}
\mathcal{R}_{\mu\nu} \equiv \hat{\Gamma}^\alpha_{\mu\nu ,\alpha} - \hat{%
\Gamma}^\alpha_{\mu\alpha , \nu} + \hat{\Gamma}^\alpha_{\alpha\lambda}\hat{%
\Gamma}^\lambda_{\mu\nu} -\hat{\Gamma}^\alpha_{\mu\lambda}\hat{\Gamma}%
^\lambda_{\alpha\nu}\,,
\end{equation}
and $S_m$ is the matter action.

However, the scalar-tensor representation of hybrid metric-Palatini gravity provides a theoretical framework which is easier to handle from a computational point of view \cite{Harko:2011nh}, where the equivalent action   to (\ref{eq:S_hybrid}) is provided by
\begin{equation}  \label{eq:S_scalar2}
S=\int \frac{d^4 x \sqrt{-g} }{2\kappa^2}\left[ (1+\phi)R +\frac{3}{2\phi}%
\partial_\mu \phi \partial^\mu \phi -V(\phi)\right]+S_m .
\end{equation}
Note a similarity with the the action of the $w=-3/2$ Brans-Dicke theory version of the Palatini approach to $f(R)$ gravity. However, the hybrid theory exhibits an important and subtle difference appearing in the scalar field-curvature coupling, which in the $w=-3/2$ Brans-Dicke theory is of the form $\phi R$ \cite{Harko:2011nh}.

By varying the action (\ref{eq:S_scalar2}) with respect to the metric provides the following gravitational field equation
\begin{equation}
G_{\mu \nu }=\kappa ^{2}\left( \frac{1}{1+\phi }T_{\mu \nu }+T_{\mu \nu
}^{(\phi )}\right) ,  \label{einstein_phi}
\end{equation}%
where $T_{\mu \nu }$ is the standard matter energy-momentum tensor, and
\begin{eqnarray}
T_{\mu \nu }^{(\phi )} &=&\frac{1}{\kappa ^{2}}\frac{1}{1+\phi }\Bigg[\nabla
_{\mu }\nabla _{\nu }\phi -\frac{3}{2\phi }\nabla _{\mu }\phi \nabla _{\nu
}\phi  \notag \\
&&\text{\hskip0.75cm}+\Bigg(\frac{3}{4\phi }\nabla _{\lambda }\phi \nabla
^{\lambda }\phi -\square \phi -\frac{1}{2}V\Bigg)g_{\mu \nu }\Bigg] ,
\label{tens_perfect}
\end{eqnarray}%
is the energy-momentum tensor of the scalar field of the theory.

Varying the action with the scalar field yields the following second-order differential equation,
\begin{equation}  \label{eq:evol-phi}
-\square \phi+\frac{1}{2\phi}\partial_\mu \phi \partial^\mu \phi +\frac{\phi[%
2V-(1+\phi)V_\phi]} {3}=\frac{\phi\kappa^2}{3}T \ .
\end{equation}
The above equation shows that in hybrid metric-Palatini gravity the scalar field is dynamical. This
represents an important and interesting difference with respect to the
standard Palatini case \cite{Olmo:2011uz}.

\subsection{Spacetime metric}

The metric of a time-dependent wormhole geometry is given by
\begin{equation}
ds^{2}=-e^{2\Phi \left( r\right) }dt^{2}+a^2\left( t\right)\left[ \frac{%
dr^{2}}{1-b\left( r\right) /r}+r^{2}d\Omega ^{2}\right] ,  \label{met}
\end{equation}%
where $\Phi \left( r\right)$ and $b\left(r\right) $, respectively, are the redshift and shape functions, which are $r$-dependent, and $a\left( t\right) $ is the time-dependent scale factor, and $d\Omega ^{2}=d\theta ^{2}+\sin^{2} \theta d\varphi ^{2}$ is the linear element of the unit sphere.

To correspond to a wormhole solution, the shape function has the following restrictions \cite{Morris:1988cz}: (i) $b(r_{0})=r_{0}$, where $r_{0}$ is the wormhole throat, corresponding to a minimum radial coordinate, (ii) $b(r)\leq r$, and (iii) the flaring-out condition given in the form $rb^{\prime }(r)-b\left( r\right) <0$. The latter condition is the fundamental ingredient in wormhole physics, as taking into account the Einstein field equation, one verifies that this flaring-out condition imposes the violation of the NEC. In fact, it violates all of the
pointwise energy conditions \cite{Morris:1988cz,Morris:1988tu,Visser:1995cc,Lobo:2017oab,Lobo:2004wq}. In order to avoid the presence of event horizons, so that the wormhole is traversable, one also imposes that the redshift function, $\Phi(r)$, be finite everywhere. However, throughout this work, we consider a zero redshift function, $\Phi =0$, which simplifies the calculations significantly.

\subsection{Evolving embedding analysis}

It is also interesting to explore how the time evolution of the scale factor affects the wormhole itself. To this effect, in the following analysis, we will follow closely the analysis outlined in Refs. \cite{Morris:1988cz,Roman:1992xj,KordZangeneh:2020jio}, and analyse how the (effective) radius of the throat or length of the wormhole changes it time. As mentioned above, throughout this paper, we consider $\Phi =0$ for simplicity. Thus, in order to analyse the time-dependent dynamic wormhole geometry, we need to choose a specific $b(r)$ to provide a reasonable wormhole at $t=0$, which is assumed to be the onset of the evolution. Note that the radial proper length between any two points $A$ and $B$, through the wormhole, at any $t={\rm const}$
is given by $l(t)=\pm \, a(t) \int_{r_A}^{r_B} {(1-b/r)^{-1/2}}\,dr$, which is merely the initial radial proper separation multiplied by the scale factor.

Now, to analyse how the ``wormhole'' form of the metric is maintained throughout the evolution, we consider a $t={\rm const}$ and $\theta=\pi/2$ slice of the spacetime (\ref{met}). Thus, the metric reduces to
\begin{equation}\label{slice}
ds^2=\frac{a^2(t)\,dr^2}{1-b(r)/r} + a^2(t)\,r^2\,d\varphi^2\,,
\end{equation}
and is embedded in a flat 3-dimensional Euclidean space given by
\begin{equation}
ds^2=d{\bar{z}}^2+d{\bar{r}}^2+{\bar{r}}^2\,{d\varphi}^2\,.
\label{barredslice}
\end{equation}
Comparing the angular coefficients, yields
\begin{eqnarray}
\bar{r}&=&{a(t)\,r}\big|_{t={\rm const}} \,,
       \label{coef1:phi}      \\
{d\bar{r}}^2&=&a^2(t)\,{dr}^2\big|_{t={\rm const}}  \,.
\label{coef2:phi}
\end{eqnarray}
Note, that when considering derivatives, these relations do not represent a ``coordinate transformation'', but a ``rescaling'' of the radial coordinate $r$ along each $t={\rm constant}$ slice.

The ``wormhole'' form of the metric will be preserved, with respect to the ${\bar{z}},{\bar{r}},\varphi$ coordinates, if the embedded slice has the following metric
\begin{equation}\label{WHslice}
ds^2=\frac{d{\bar{r}}^2}{1-\bar{b}(\bar{r})/\bar{r}} +
                      {\bar{r}}^2{d\varphi}^2\,,
\end{equation}
and the shape functions $\bar{b}(\bar{r})$ has a minimum at a specific $\bar{b}(\bar{r}_0)=\bar{r}_0$. Now, the embedded slice  (\ref{slice}) can be readily rewritten in the form of Eq. (\ref{WHslice}) by using Eqs. (\ref{coef1:phi}) and (\ref{coef2:phi}), and the following relation
\begin{equation}
\bar{b}(\bar{r})=a(t)\,b(r).
    \label{bar:b}
\end{equation}
Thus, the evolving wormhole will have the same overall shape and size and relative to the ${\bar{z}},{\bar{r}},\varphi$ coordinate system, as the initial wormhole had relative to the initial $z,r,\varphi$ embedding space coordinate system.
In addition to this, using Eqs. (\ref{barredslice}) and (\ref{WHslice}), one deduces that
\begin{equation}
{{d{\bar{z}}}\over{d{\bar{r}}}}
=\pm\left({{\bar{r}}\over{\bar{b}(\bar{r})}}-1\right)^{-1/2}
={{dz}\over{dr}} \,,
            \label{barredembedding}
\end{equation}
which implies
\begin{eqnarray}
\bar{z}(\bar{r})&=&\pm\int{{d\bar{r}}\over{(\bar{r}/{\bar{b}(\bar{r})}-1)^{1/2}}}
	\nonumber \\
&=& \pm \, a(t)\,\int{\left(\frac{r-b}{b}\right)^{-1/2}} \,dr
	\nonumber  \\
&=&\pm \, a(t)\,z(r)\,.
            \label{embed:relation}
\end{eqnarray}
Thus, the relation between the initial embedding space at $t=0$ and the embedding space at any time $t$ is, from Eqs. (\ref{coef2:phi}) and (\ref{embed:relation})
\begin{equation}
ds^2=d{\bar{z}}^2+d{\bar{r}}^2+{\bar{r}}^2\,{d\varphi}^2
        =a^2(t)\,[dz^2+dr^2+r^2{d\varphi}^2]\,.
\end{equation}
The wormhole will always remain the same size, relative to the ${\bar{z}},{\bar{r}},\phi$ coordinate system, as the scaling of the embedding space compensates for the evolution of the wormhole.
Nevertheless, the wormhole will change size relative to the initial $t=0$ embedding space.

We also note that the ``flaring out condition'' for the evolving wormhole is given by
\begin{equation}
{{d\,^2{\bar{r}(\bar{z})}}\over{d{\bar{z}}^2}}>0\,,
\end{equation}
at or near the throat, so that taking into account Eqs. (\ref{coef1:phi}),
(\ref{coef2:phi}), (\ref{bar:b}), and (\ref{barredembedding}), it
follows that
\begin{equation}
{{d\,^2{\bar{r}(\bar{z})}}\over{d{\bar{z}}^2}}
       =\frac{1}{a(t)}\,{{b-b'r}\over{2b^2}}
       =\frac{1}{a(t)}\,{{d\,^2r(z)}\over{dz^2}}>0\,,
       \label{barred:flareout}
\end{equation}
at or near the throat. Using Eqs. (\ref{coef1:phi}),
(\ref{bar:b}), and
\begin{equation}
{\bar{b}}'(\bar{r})={{d\bar{b}}\over{d\bar{r}}}
          =b'(r)={{db}\over{dr}}\,,
\end{equation}
the right-hand side of Eq. (\ref{barred:flareout}), relative to the ${\bar{z}},{\bar{r}},\phi$ coordinate system, may be written as
\begin{equation}
{{d\,^2{\bar{r}(\bar{z})}}\over{d{\bar{z}}^2}}
=\left({{\bar{b}-{\bar{b}}'\bar{r}}\over{2{\bar{b}}^2}}\right)>0\,,
         \label{barred:flareout2}
\end{equation}
at or near the throat. Thus, we have verified that the flaring out condition (\ref{barred:flareout2}), using the barred coordinates, has the same form as for the static wormhole.

\subsection{Gravitational field equations}

We also consider an anisotropic matter energy-momentum tensor given by $T_{\nu }^{\mu }={\rm diag}(-\rho,-\tau, p, p)$, where $\rho $, $\tau $
and $p$ are the energy density, the radial tension (which is equivalent to a negative radial pressure) and the tangential pressure, respectively. Using the metric (\ref{met}), the gravitational
field equations (\ref{einstein_phi}) provide the following energy-momentum profile:
\begin{equation}
\rho (t,r) =\left( 3H^{2} +\frac{b^{\prime }}{r^{2}a^{2}}%
\right) \left( 1+\phi \right) +3\dot{\phi}H+\frac{3\dot{\phi}^{2}}{4\phi }-%
\frac{1}{2}V(\phi ),
~  \label{rhoeq}
\end{equation}
\begin{eqnarray}
\tau \left( t,r\right) &=&\left( H^{2}+2\frac{\ddot{a}}{a} +\frac{b}{%
r^{3}a^{2}}\right) \left( 1+\phi \right)  \notag \\
&&+2\dot{\phi}H-\frac{3\dot{\phi}^{2}}{4\phi }+\ddot{\phi}%
-\frac{1}{2}V(\phi ),
\end{eqnarray}
\begin{eqnarray}
p\left( t,r\right) &=&\left( -H^{2}-2\frac{\ddot{a}}{a} -\frac{rb^{\prime}-b}{%
2r^{3}a^{2}}\right) \left( 1+\phi \right)  \notag \\
&&-2\dot{\phi}H+\frac{3\dot{\phi}^{2}}{4\phi }-\ddot{\phi}%
+\frac{1}{2}V(\phi ),  \label{pteq}
\end{eqnarray}
respectively, where $H=\dot{a}/a$ and $\phi =\phi \left( t\right) $, the overdot and prime denote derivatives with respect to $t$ and $r$, respectively, and for notational simplicity, we have considered $\kappa =1$. One recovers the standard field equations of the Morris-Thorne wormhole \cite{Morris:1988cz}, by fixing the scale factor $a$ to unity and excluding the background time-dependent evolution and the scalar field $\phi $ contribution.

Furthermore, we will also analyse the null and weak energy conditions for the solutions obtained below. The weak energy condition (WEC) is defined as $T_{\mu \nu
} u^{\mu }u^{\nu }\geq 0$, where $u^{\mu }$ is a timelike vector, and is
expressed in terms of the energy density $\rho $, radial tension $\tau $ and
tangential pressure $p$ as $\rho \geq 0$, $\rho -\tau \geq 0$ and $\rho
+p\geq 0$, respectively. The last two inequalities, i.e., $\rho -\tau \geq 0$
and $\rho +p\geq 0$ correspond to the NEC, which is defined as $T_{\mu \nu
}k^{\mu }k^{\nu }\geq 0$, where $k^{\mu }$ is \textit{any} null vector.

The scalar field equation of motion (\ref{eq:evol-phi}) yields
\begin{equation}
\ddot{\phi}+3\dot{\phi}H-\frac{\dot{\phi}^{2}}{2\phi }-\frac{1}{3}\phi \left[
T +\left( 1+\phi \right) \frac{dV}{d\phi }-2V(\phi )\right]
=0,  \label{phieom}
\end{equation}%
where the trace of the energy-momentum tensor $T=T\left( t,r\right) =T_{\mu }^{\mu }=-\rho -\tau +2p$ is given by
\begin{eqnarray}
T\left( t,r\right) = -2\left( \frac{b^{\prime }}{a^{2}r^{2}}+3H^{2}+3\frac{%
\ddot{a}}{a}\right) (1+\phi )  \notag \\
-3\ddot{\phi}+\frac{3\dot{\phi}^{2}}{2\phi }-9\dot{\phi}%
H+2V(\phi ).  \label{trace}
\end{eqnarray}%
In order to solve the differential equation (\ref{phieom}) for $\phi \left(
t\right) $, the trace $T\left( t,r\right) $, given by Eq. (\ref{trace}),
should be independent of $r$.  This condition leads to $b^{\prime }/r^{2}=CH_{0}^{2}$
where $H_{0}$ is the present value of the Hubble parameter and{\ $C$ is an
arbitrary dimensionless constant. Thus, the shape function is given by $b(r)=r_{0}+C{H_{0}^{2}}%
(r^{3}-r_{0}^{3})/3$, which satisfies $b(r_{0})=r_{0}$, and the flaring-out condition at the throat imposes $C H_{0}^{2} r_{0}^{2} < 1$.

To solve Eq. (\ref{phieom}) numerically for $\phi$, it is useful to rewrite it in terms of dimensionless functions of the scale factor $a$. Thus, we consider the definitions:
\begin{eqnarray}
\ddot{\phi}=\ddot{a}\phi ^{\prime }(a)+\dot{a}^{2}\phi ^{\prime \prime }(a), \qquad  \dot{\phi}=\dot{a}\phi ^{\prime }(a),
	\nonumber \\
\ddot{a}=H\dot{a}+\dot{H}a, \qquad  \dot{H}=\dot{a}H^{\prime }(a), \qquad \dot{a}=Ha	 \,.
\nonumber
\end{eqnarray}
Note that in the above definitions, the prime denotes a derivative with respect to the scale factor. In addition to this, consider $U=V/3H_{0}^{2}$
and $E=H/H_{0}$, so that Eq. (\ref{phieom}) finally takes the following form
\begin{eqnarray}
\phi ^{\prime \prime }(a)-\frac{\phi ^{\prime 2}}{2\phi (a)}+\frac{%
4(a\phi ^{\prime }(a) + \phi(a))}{a^2}+\frac{E^{\prime }(a)\phi ^{\prime }(a)}{E(a)}  \notag
 \\
+\frac{2C\phi (a)}{3a^{4}E(a)^{2}}+\frac{2\phi
(a)E^{\prime }(a)}{aE(a)}-\frac{\phi (a)}{a^{2}E(a)^{2}}\frac{dU}{d\phi } =0.
\label{phidleom}
\end{eqnarray}

\subsection{Barotropic equation of state}

To solve the differential equation (\ref{phidleom}) for $\phi \left( a\right) $, we need to deduce $E$. To this effect, it is also useful to define the following background quantities:
\begin{equation}
\rho _{b}\left( t\right) =3H^{2}, \qquad  \tau _{b}\left( t\right) =H^{2}+2\frac{\ddot{a}}{a}\,,
\end{equation}
and consider a background barotropic equation of state given by $\tau _{b}=-\omega _{b}\rho_{b} $. From this condition, we find $E=a^{-3\left( \omega _{b}+1\right) /2}$. One could solve this equation of state for scale factor as a function of the time coordinate $t$, which yields $a(t) \propto t^{2/3(1+\omega_{b})}$. For $\omega_{b} = -1$, the scale factor takes the exponential form in terms of $t$.
Furthermore, we set $C=0$ which satisfies the flaring-out condition at the throat (see discussion below Eq. (\ref{trace})), and in addition simplifies the solution, so that $b(r)=r_{0}$.
Thus, taking into these conditions and using the dimensionless definitions given above, Eqs. (\ref{rhoeq})-(\ref{pteq}) lead to the following relations
\begin{equation}
\frac{\rho }{3H_{0}^{2}} =
\frac{ 1+\phi (a)}{a^{3\left( \omega _{b}+1\right) }}
+\frac{\phi ^{\prime }(a)%
}{a^{3\omega _{b}+2}}
+\frac{\phi ^{\prime 2}(a)}{4a^{3\omega _{b}+1}\phi (a)}-\frac{1}{2}U(\phi
), \label{rho}
\end{equation}
\begin{eqnarray}
&&\frac{\rho -\tau }{3H_{0}^{2}} =
\left( \frac{\omega _{b}+1}{a^{3\left(\omega _{b}+1\right) }}-\frac{r_{0}}{3a^{2}H_{0}^{2}r^{3}}\right) \left( 1+\phi(a)\right)
  \nonumber \\
&&\qquad +\frac{\left(\omega _{b}+1 \right) \phi ^{\prime }(a)}{2a^{3\omega _{b}+2}}%
+\frac{\phi ^{\prime 2}(a)}{2a^{3\omega _{b}+1}\phi (a)}-\frac{\phi ^{\prime
\prime }(a)}{3a^{3\omega _{b}+1}},  \label{rhotau}
\end{eqnarray}
\begin{eqnarray}
&&\frac{\rho +p}{3H_{0}^{2}} =
\left( \frac{\omega _{b}+1}{a^{3\left( \omega_{b}+1\right) }}+\frac{r_{0}}{6a^{2}H_{0}^{2}r^{3}}\right) \left( 1+\phi (a)\right)
\nonumber \\ &&
\qquad +\frac{\left(\omega _{b}+1 \right) \phi ^{\prime }(a)}{2a^{3\omega _{b}+2}}%
+\frac{\phi ^{\prime 2}(a)}{2a^{3\omega _{b}+1}\phi (a)}-\frac{\phi ^{\prime
\prime }(a)}{3a^{3\omega _{b}+1}},  \label{rhop}
\end{eqnarray}%
respectively. It is interesting to note that $\rho$ is independent of the radial
coordinate $r$. To keep the terms dimensionless, we will consider the wormhole throat as $r_{0}= A H_{0}^{-1}$, where $A$ is a dimensionless constant. In what follows, we set $A$ to unity.

In the following, we shall consider specific solutions for the parameters $\omega _{b}=(-1, 1/3,0)$, respectively. These parameters range tentatively correspond to the inflationary, radiation and matter epochs.

\section{Specific evolving wormhole solutions and the energy conditions}\label{secIII}

\begin{figure*}[htb!]
\centering%
\subfigure[~Specific case of $\omega_b=-1$  with $U(\phi)=\phi^{0.5}$] {\label{fig1a}%
\includegraphics[width=.32\textwidth]{fig1a.eps} }
\subfigure[~Specific case of $\omega_b=1/3$ with
$U(\phi)=\phi^{2}$] {\label{fig1b}\includegraphics[width=.32\textwidth]{fig1b.eps} }
\subfigure[~Specific case of $\omega_b=0$ with $U(\phi)=\phi^{2}$] {\label{fig1c}\includegraphics[width=.32\textwidth]{fig1c.eps}}
\caption{The behaviors of $\phi$ and $\phi'(a)$ vs $z$ for the specific case of  $\omega_b=-1$, $\omega_b=1/3$, and $\omega_b=0$, respectively. Note that both horizontal and vertical axes are logarithmic.}
\label{fig1}
\end{figure*}

\begin{figure*}[htb!]
\centering\subfigure[~$\rho/3H_0^2$ vs $z$] {\label{fig2a}%
\includegraphics[width=.32\textwidth]{fig2a.eps} }
\subfigure[~$ \abs{\rho-\tau}/3H_0^2$ vs $z$] {\label{fig2b}\includegraphics[width=.32\textwidth]{fig2b.eps}}
\subfigure[~$(\rho+p)/3H_0^2$ vs $z$] {\label{fig2c}%
\includegraphics[width=.32\textwidth]{fig2c.eps} }
\caption{The behaviors of $\protect\rho $, $\protect\rho -\protect\tau $ (which is negative) and
$\protect\rho +p$, respectively, versus $z$ for different values of $r$ for $\protect\omega_b=-1$ with $U(\protect\phi)=\protect\phi^{0.5}$. Note
that both horizontal and vertical axes are logarithmic, except for the vertical axis of Fig. \ref{fig2a}.}
\label{fig2}
\end{figure*}

\begin{figure*}[htb!]
\centering
\subfigure[~$\rho/3H_0^2$ vs $z$] {\label{fig3a} \includegraphics[width=.32\textwidth]{fig3a.eps}}
\subfigure[~$(\rho-\tau)/3H_0^2$ vs $z$] {\label{fig3b}\includegraphics[width=.32\textwidth]{fig3b.eps}
} \subfigure[~$(\rho+p)/3H_0^2$ vs $z$] {\label{fig3c}\includegraphics[width=.32\textwidth]{fig3c.eps} }
\caption{The behaviors of $\protect\rho $, $\protect\rho -\protect\tau $ and
$\protect\rho +p$, respectively, versus $z$ for different values of $r$ in
the specific case of $\protect\omega_b=1/3$ with $U(\protect\phi)=%
\protect\phi^{2}$. Note that both the horizontal and vertical
axes are logarithmic.}
\label{fig3}
\end{figure*}

\begin{figure*}[htb!]
\centering\subfigure[~$\rho/3H_0^2$ vs $z$] {\label{fig4a}%
\includegraphics[width=.32\textwidth]{fig4a.eps} }
\subfigure[~$(\rho-\tau)/3H_0^2$ vs $z$] {\label{fig4b}\includegraphics[width=.32\textwidth]{fig4b.eps}
} \subfigure[~$(\rho+p)/3H_0^2$ vs $z$] {\label{fig4c}\includegraphics[width=.32\textwidth]{fig4c.eps} }
\caption{The behaviors of $\protect\rho $, $\protect\rho -\protect\tau $ and
$\protect\rho +p$, respectively, versus $z$ for different values of $r$ in
the specific case of $\protect\omega_b=0$ with $U(\protect\phi)=%
\protect\phi^{2}$. Both horizontal and vertical axes are
logarithmic.}
\label{fig4}
\end{figure*}

\begin{figure*}[htb!]
\centering
\subfigure[~$\phi$ vs $z$] {\label{fig5a}\includegraphics[width=.32\textwidth]{fig5a.eps}}
\subfigure[~$\phi'(a)$ vs $z$] {\label{fig5b}\includegraphics[width=.32\textwidth]{fig5b.eps}}
\subfigure[~$E$ vs $z$] {\label{fig5c}\includegraphics[width=.32\textwidth]{fig5c.eps}}
\caption{The behaviors of $\phi$, $\phi'(a)$ and $E$ vs $z$ for the traceless EMT case with $U(\protect\phi)=\protect\phi^{2}$. Note that both horizontal and vertical axes are logarithmic, except for the vertical axis of Fig. \ref{fig5a}.
See the text for more details.}
\label{fig5}
\end{figure*}

\begin{figure*}[htb!]
\centering\subfigure[~$\rho/3H_0^2$ vs $z$] {\label{fig6a}\includegraphics[width=.32\textwidth]{fig6a.eps} }
\subfigure[~$(\rho-\tau)/3H_0^2$ vs $z$] {\label{fig6b}\includegraphics[width=.32\textwidth]{fig6b.eps}
} \subfigure[~$(\rho+p)/3H_0^2$ vs $z$] {\label{fig6c}\includegraphics[width=.32\textwidth]{fig6c.eps} }
\caption{The behaviors of $\protect\rho $, $\protect\rho -\protect\tau $ and
$\protect\rho +p$, respectively, versus $z$ for different values of $r$ for a wormhole with traceless EMT with $U(\protect\phi)=\protect\phi^{2}$. Both horizontal and vertical axes are logarithmic.}
\label{fig6}
\end{figure*}

In this section, we analyse evolving traversable wormhole geometries for the specific parameters $\omega _{b}=(-1, 1/3,0)$ considered above, as well as evolving wormholes with a traceless EMT, and study the NEC and WEC for the solutions obtained.

Before going on further, we would like to comment on the potential for the scalar field $U(\phi)$ and its effects on the results. As Eqs. (\ref{rhotau}) and (\ref{rhop}) (and (\ref{rhoeq})-(\ref{pteq})) show, the quantities $\rho-\tau$ and $\rho+p$ (and consequently the NEC) are independent of the scalar field potential. However, it affects the behavior of $\rho$. We will consider a power law potential for scalar field $U(\phi)=\phi^\alpha$. There are some reasonable choice for $\alpha$, such as $\alpha=2$ and $\alpha=4$, inspired by a mass term and the Higgs potential, respectively. In fact, we emphasize that all the results below hold qualitatively for potentials of the form $\phi^2$ and $\phi^4$. For the sake of economy, we will present the results for the quadratic case, $U(\phi)=\phi^2$. In addition to this, to verify that other choices for $\alpha$ could also present reliable results, we consider the specific choice $\alpha=0.5$ in one of the cases.

\subsection{Specific case: $\omega _{b}=-1$}

Here, we consider the specific case $\omega _{b}=-\tau _{b}/\rho _{b}=-1$ with $U\left( \phi \right) =V\left(\phi \right) /3H_{0}^{2}=\phi^{0.5}$. Figure \ref{fig1a} depicts the behaviors of $\phi$ and $\phi^{\prime}$ with respect to the redshift $z$ $(=1/a-1)$.
In Fig. \ref{fig2}, the behaviors of $\rho $, $\rho -\tau $ and $\rho +p$ versus $z$ for
different values of $r$ are shown.
As one can see, $\rho - \tau$ is negative throughout the entire evolution. We emphasize that we have imposed different possible initial values and obtained the same result. Note that one may justify this outcome in the following manner. For this case with $\omega _{b}=-1$, Eq. (\ref{rhotau}) reduces to
\begin{equation}
\frac{\rho -\tau }{3H_{0}^{2}} =
- \frac{\left( 1+\phi (a)\right)r_{0}}{3a^{2}H_{0}^{2}r^{3}}
+a^{2}\left( \frac{\phi ^{\prime 2}(a)}{2\phi (a)}-\frac{\phi ^{\prime
\prime }(a)}{3}\right).
\end{equation}
Since $a\leq 1$, the term with a coefficient $1/a^2$ is dominant. This term is therefore negative, rendering $\rho -\tau $ negative. Thus, for this specific case, both the NEC and the WEC are violated.

\subsection{Specific case: $\omega _{b}=1/3$}

Figure \ref{fig1b} illustrates the behavior of  $\phi$ and $\phi^{\prime}$ versus $z$ for the specific case $\omega _{b}=1/3$. For this case, $E=a^{-2}$, $a(t) \propto t^{1/2}$ and we consider $U\left( \phi \right)=\phi^{2}$. The behaviors of $\rho $, $\rho -\tau $ and $\rho +p$ with respect to the redshift $z$ for different values of $r$ are displayed in Fig. \ref{fig3}, which shows explicitly
that these quantities decrease as time evolves. However, they remain positive at all times and consequently
the NEC and WEC are always satisfied. This occurs for the wormhole throat as
well as for other wormhole radii. It is interesting to note that at a specified
time/redshift, the quantity $\rho -\tau $ increases for increasing values of
the radius, and the minimum value corresponds to the throat, as depicted by
Fig. \ref{fig3b}. On the other hand, $\rho+p$ decreases for increasing
values of the radius, and has a maximum at the throat (Fig. \ref{fig3c}).}

\subsection{Specific case: $\omega _{b}=0$}

Here we intend to study the energy conditions for the case $\omega _{b}=0$, for which $E=a^{-3/2}$ and $a(t)\propto t^{2/3}$. One can verify the behaviors of the scalar field $\phi$ and its derivative with respect to $a$ as functions of the redshift $z$ in Fig. \ref{fig1c}, where we consider the power law potential corresponding to scalar field as $U\left( \phi \right) =\phi ^{2}$.
The energy conditions are depicted in Fig. \ref{fig4}, where one verifies that $\rho $, $\rho -\tau $ and $\rho +p$ are positive at all times and for different values of spatial coordinate $r$. These quantities are also decreasing (increasing) functions with respect to time (redshift). Thus, the NEC and WEC are always satisfied for all wormhole radii including the wormhole throat. The behaviors of $\rho -\tau $ and $\rho + p $ at a specified time/redshift are analogous to the previous case with $\omega _{b}=1/3$, namely, the quantity $\rho -\tau $ ($\rho + p $) increases (decreases) for increasing values of the radius (see Figs. \ref{fig4b} and \ref{fig4c}).

\subsection{Wormholes with traceless EMT}

In this subsection, we consider the traceless energy-momentum tensor, with $T=
-\rho -\tau +2p=0$, which taking into account the dimensionless quantities defined above, provides the following differential equation
\begin{eqnarray}  \label{tracelessEMTb}
 -3 a^2 E(a)^2 \left[\phi ^{\prime \prime }(a)-\frac{\phi ^{\prime}(a)^{2}}{2 \phi (a)}\right]
-3 a E(a) \big[ a E^{\prime }(a)
\notag \\
+4 E(a) \big] \phi^{\prime}(a)
	-6 E(a) \left[ a  E^{\prime }(a)+2 E(a)\right]
\notag \\
\times \left[1+\phi (a) \right]	+6 U(\phi)=0.
\end{eqnarray}
Here, we solve the system of coupled differential equations (\ref{phidleom}) and (\ref%
{tracelessEMTb}), for $\phi $ and $E$, numerically. The behaviors of $\phi$, $\phi^{\prime}$ and $E$ versus $z$, taking into account the choice $U\left( \phi \right) =\phi^{2}$, are shown in Fig. \ref{fig5}. Note that $E$ tends to unity at the present time ($z \rightarrow 0$), as expected (Fig. \ref{fig5c}).
In Fig. \ref{fig6}, the behaviors of $\rho $, $\rho -\tau $ and $\rho +p$, with respect to $z$ for
different values of $r$ are displayed.
Note that $\rho $, $\rho -\tau $ and $\rho +p$ are decreasing (increasing)
functions of time (redshift) and are positive for all radii, as time evolves.
Thus, the NEC and WEC are satisfied at all times.
Moreover, as shown in Figs. \ref{fig6b} and \ref{fig6c}, the qualitative behaviors of $\rho -\tau $ and $\rho + p$ with respect to $r$ is analogous to the previous two cases, namely, $\omega _{b}=1/3$ and $\omega _{b}=0$, where the quantity $\rho -\tau $ ($\rho+p$) has a minimum (maximum) value at throat.

\section{Discussion and Conclusion}\label{conclusion}

In the present paper, we have studied the evolution of dynamic traversable wormhole geometries in a FLRW background in the context of hybrid metric-Palatini gravity. This theory, which recently attracted much attention, consists of a hybrid combination of metric and Palatini terms and is capable of avoiding several of the problematic issues associated to each of the metric or Palatini formalisms (for more details, we refer the reader to Refs. \cite{Capozziello:2015lza,Harko:2018ayt,Harko:2020ibn}). For the evolving wormholes, we presented the components of the energy-momentum tensor that supports these geometries in terms of the model's functions, namely, the scalar field, the scale factor and the shape function (we considered a zero redshift function, for simplicity). Furthermore, we found specific wormhole solutions by considering a barotropic equation of state for the background matter, i.e., $\tau _{b} = - \omega _{b} \rho _{b}$, and considered particular equation of state parameters.

More specifically, we showed that for the specific cases of $\omega _{b}=1/3$ and $\omega _{b}=0$, the entire wormhole matter satisfies the NEC and WEC for all times. The latter cases are similar to the wormhole geometries analysed in the presence of pole dark energy \cite{KordZangeneh:2020jio}, however, there the WEC is violated at late times, i.e., the energy density becomes negative. Thus, the present results outlined in this work strengthen the varying dark energy models and may suggest that hybrid metric-Palatini gravity is a rather more promising model to explore.
In addition to the barotropic equation of state, we also studied evolving wormhole geometries supported by the matter with a traceless EMT. For this specific geometry, we discovered that both the NEC and the WEC are satisfied at all times as well. These results are extremely promising as they build on previous work that consider that the energy conditions may be satisfied in specific flashes of time \cite{Kar:1994tz,Kar:1995ss}.

An interesting astrophysical observational aspect on how to detect these wormholes would be to analyse the physical properties and characteristics of matter forming thin accretion disks around the wormhole geometries analysed in this work, much in the spirit of the analysis carried out in Refs. \cite{Harko:2008vy,Harko:2009xf,Harko:2009gc,Harko:2017fra}. In fact, specific signatures could appear in the electromagnetic spectrum, thus leading to the possibility of distinguishing these wormhole geometries by using astrophysical observations of the emission spectra from accretion disks. This would be an interesting avenue of research to explore.

\begin{acknowledgements}
We thank the referees for the constructive comments that helped us to significantly improve the paper.
MKZ would like to thank Shahid Chamran University of Ahvaz, Iran for
supporting this work. FSNL acknowledges support from the Funda\c{c}\~{a}o
para a Ci\^{e}ncia e a Tecnologia (FCT) Scientific Employment Stimulus
contract with reference CEECIND/04057/2017, and thanks funding from
the research grants No. PTDC/FIS-OUT/29048/2017, No. CERN/FIS-PAR/0037/2019 and
No. UID/FIS/04434/2020.
\end{acknowledgements}


\section{Introduction}
\label{sec:introduction}

There has been rapid progress in deep learning models for a variety of computer vision tasks; such models have increasingly been introduced into real-world applications. It is crucial to quantify the error of model predictions in safety-critical domains such as assisted driving systems and health care. For instance, an autonomous car steering system first processes the sensory inputs, then feeds them to the deep learning models (e.g., depth estimation and semantic segmentation), and finally conducts higher-level decision making based on the models' output ~\cite{bojarski2016end,huval2015empirical}. If the models make nearly random guesses when faced with uncertainty, the result could be a fatal accident. The systems can make reliable high-level decisions if it is given information about the uncertainty.

Previous studies on Bayesian modeling divided uncertainty into epistemic and aleatoric types~\cite{gal2016uncertainty}. Epistemic uncertainty reflects a lack of data and knowledge: for instance, a model trained on a depth dataset of rural areas might perform poorly on highways, owing to its lack of knowledge about such environments. Epistemic uncertainty can be reduced given sufficient data. By contrast, aleatoric uncertainty is intrinsic to an observation, and adding more data may not reduce it.

\begin{figure}[!t]
\begin{center}
    \includegraphics[width=0.97\columnwidth]{./images/introduction_yu.png}
    \caption{{\bf Outline of the network of  our method}, using the example of age estimation. Here $y$ and $w$ are learned to predict the targets and their certainty, respectively. $L_{t}$ is the loss for the target estimation, and is calculated using the target $y^*$ and its prediction $y$. $L_{u}$ is the loss for the certainty estimation, where we propose using the virtual residual $\tilde{r}_i$ to optimize the certainty estimation $w_i$.}
    \label{fig:introduction}
  \end{center}
\end{figure}

We focus on the aleatoric uncertainty, which cannot be mitigated by increasing the amount of data. Aleatoric uncertainty is further subdivided into homoscedastic uncertainty and heteroscedastic uncertainty. The former is constant for different inputs, while the latter is not. Thus, predicting the heteroscedastic uncertainty from the inputs is particularly important~\cite{NIPS2017_7141}. Unless specified otherwise, the term ``uncertainty'' in this paper refers to heteroscedastic aleatoric uncertainty, i.e., the data-dependent uncertainty from observation noise.

Previous methods for aleatoric uncertainty estimation have used a Bayesian deep learning framework, in which neural networks predict not only regressed values $y$, but also the uncertainty in the model's prediction~\cite{gal2016uncertainty,NIPS2017_7141}. These methods optimize the estimation of the target and uncertainty at the same time -- the uncertainty estimation is defined based on the target error for the training data. However, usually the error for the test data is far larger than that for the training data, which the model tends to overfit. Thus, in Bayesian formulations, the overfitting of the signal estimation is problematic: when overfitting occurs, the model underestimates errors. There are many ways to mitigate the overfitting, but practically it cannot be removed entirely\cite{overfitting, dropout, stopped}.

In this paper, we present an optimization framework for aleatoric uncertainty in regression problems by separating the estimation of the target from its uncertainty. In contrast to previous studies, our method allows the network to learn to predict the signal $y$ and certainty $w$ by optimizing the separated loss $L(y, w)=L_t(y)+\lambda L_u(w)$, as shown in Figure~\ref{fig:introduction}. Here, $L_t(y)$ and $L_u(w)$ are the losses of target estimation and certainty estimation respectively, and $\lambda$ is a parameter to balance the weighting of the two terms. Note that the certainty $w$ is estimated instead of the uncertainty itself. We introduce the strict definition of $w$ in Section\ref{sec:related_work}. Our method addresses the shortcomings of existing approaches in two respects. First, we can adjust the two performances of regression and uncertainty estimation by controlling $\lambda$. Second, unlike previous studies of uncertainty estimation, our method explicitly avoids the effect of overfitting. By making use of other models trained with specific samples, we successfully avoid underestimating the error. To demonstrate the effectiveness of our proposed framework, we conduct experiments on three different tasks: simple simulation, age estimation, and depth estimation. Our contributions are as follows:
\begin{itemize}
\item We propose a novel optimization framework for aleatoric uncertainty estimation by a separate formulation, which allows adjusting the weight between the predictive task and the uncertainty estimation, and avoids the effect of overfitting.
    \item We introduce new metrics that measure the goodness of error quantification.
    \item We evaluate our approach on different tasks including age estimation and depth estimation. We find that our method outperforms the state-of-the-art existing methods in terms of error estimation.
\end{itemize}
In our earlier preliminary work, we have proposed a separable formulation without virtual residuals~\cite{asai}. In the last section of this paper, we argue for the benefit of introducing the virtual residuals by comparing our method to the separable formulation without virtual residuals.

\section{Related Work}
\label{sec:related_work}
Although numerous studies have estimated epistemic uncertainty~\cite{gal2016dropout,teye2018bayesian,lakshminarayanan2017simple,Postels2019ICCV}, aleatoric uncertainty has been given less attention. Previous studies have conducted aleatoric uncertainty estimation using Bayesian deep learning frameworks~\cite{gal2016uncertainty,NIPS2017_7141} that train networks to predict data-dependent noise variance from inputs. In Bayesian deep learning, given the input $x$ and the model parameter $\theta$, the network outputs the signal estimation and uncertainty estimation.
The loss function is given by
\begin{equation}
  \mathcal{L}(x, \theta) = \frac{1}{N}\sum_{i=1}^N \left\{ \frac{|y_i^* - y_i(x; \theta)|^2}{2\sigma_i(x, \theta)^2} + \frac{1}{2}\log{\sigma_i(x; \theta)^2} \right\} ,
  \label{eq:kendall_gaussian_original}
\end{equation}
where $y^*$ is the ground truth of $y$. Minimizing this loss function means maximizing the log likelihood of a Gaussian distribution. We call this type of formulation a ``joint formulation'' (JF). Such formulations have previously been used to study applications of heteroscedastic uncertainty ~\cite{bloesch2018codeslam, feng2018leveraging}.

Maximizing the log likelihood of a Laplace distribution rather than the Gaussian distribution of Eq.~\ref{eq:kendall_gaussian_original}, we can use the loss function as below:
\begin{equation}
  \mathcal{L}(x, \theta) = \frac{1}{N}\sum_{i=1}^N \left\{ \frac{|y_i^* - y_i(x; \theta)| }{b_i(x, \theta)} + \log{b_i(x; \theta)} \right\} ,
  \label{eq:kendall_laplace_original}
\end{equation}
Here, $b$ is the scale of the distribution. Kendall et al.~\cite{NIPS2017_7141} used this optimization in their experiments, because it applies $L1$ loss on the target estimation residuals, and seems to outperform the optimization by Eq.~\ref{eq:kendall_gaussian_original}, which applies $L2$ loss.

Eq.~\ref{eq:kendall_laplace_original} can be re-written as Eq.~\ref{eq:kendall_laplace} below, and the network outputs the signal estimation $y$ and the certainty estimation $w$:
\begin{equation}
  \mathcal{L}(x, \theta) = \frac{1}{N}\sum_{i=1}^N \left\{ \exp(w_i(x; \theta)) {r_i(x; \theta)} - w_i(x; \theta)\right\}.
  \label{eq:kendall_laplace}
\end{equation}
In this equation, we introduce two variables: $r_i$ is the absolute value of residual, and $w_i$ is the certainty estimation:
\begin{align}
r_i &= |y_i(x; \theta) - y^*_i|, \label{eq:r_definition} \\
w_i &= -\log b_i. \label{eq:w_definition}
\end{align}

\section{The Influence of Overfitting}
\label{sec:overfitting}
The certainty estimation can be used to estimate the expected value of the error of the target estimation. From the partial derivative of Eq.~\ref{eq:kendall_laplace} at $w_i$, the optimal $w^{opt}_i$ satisfies the following equation:
\begin{equation}
  r_i(x) = \exp(-w^{opt}_i(x)).
  \label{eq:uncertainty_optimization}
\end{equation}
Thus, when we get $y_i$ and $w_i$ as outputs, we can estimate the error of $y_i$ using Eq.~\ref{eq:uncertainty_optimization}. Error estimation plays an important role when we want to know how much we can trust the estimation. However, overfitting hinders using $w_i$ for this purpose.

\begin{figure}[!t]
\begin{center}
    \includegraphics[width=0.75\columnwidth]{./images/error_compare.png}
    \caption{{\bf The distribution of errors in an age estimation task}. The blue and red bars visualize the distribution of the estimation errors of the training and test data, respectively. The estimation was conducted by ResNet50~\cite{he2016deep} and optimized by $L1$ loss. The specifics of the implemaentation are the same as in~\ref{subsec:age_estimation}.}
    \label{fig:error_compare}
  \end{center}
\end{figure}

Fig.~\ref{fig:error_compare} shows the difference of the distribution of the errors between the training and test data for an age estimation task. The estimation errors of the training data are very small -- the majority of them have errors less than one. By contrast, estimation errors are large in the test data. This clearly shows the effect of overfitting, which makes the error of the training data smaller than that of the test data. In its training, we used the validation data, paying careful attention to avoid increasing the error of the validation data. We finish training before the validation error starts ascending. Although we carefully trained the model, the effect of overfitting remained.

Once overfitting occurs, $r_i$, the error of the target estimation during training, becomes smaller than it should be. Overfitting will make the uncertainty estimation too small if $w_i$ is optimized based on the value of $r_i$ as in the existing method~\cite{NIPS2017_7141}. The aim of our method is to estimate uncertainty suppressing the effect of overfitting.

\section{Proposed Method}
\label{sec:proposed_method}

\subsection{Separate Formulation of Aleatoric Uncertainty}
\label{subsec:separate_formulation}
In this section, we present an optimization framework for uncertainty estimation in regression problems. We focus on training a network that predicts the signal value $y_i$ and its certainty $w_i$.

We cannot balance the weighting of the target estimation and the uncertainty estimation in Eq.~\ref{eq:kendall_laplace} because the two variables for optimization are inseparable. Thus, we propose decoupling the loss function into two terms:
\begin{equation}
  \mathcal{L}(x, \theta) = \frac{1}{N}\sum_{i=1}^N \left\{ L_t(y_i) + \lambda L_u(w_i) \right\}.
  \label{eq:ours_abstract}
\end{equation}
Here, $L_t(y_i)$ is the loss for the target estimation, and $L_u(w_i)$ is that for the uncertainty estimation. By using this formulation, we can balance the target and uncertainty estimation terms by adjusting $\lambda$ appropriately. We define $L_t(y_i) = r_i(x; \theta)$, and $L_u(w_i) = \exp(w_i(x; \theta)) \tilde{r}_i(x; \theta) - w_i(x; \theta)$, resulting in the loss function
\begin{equation}
\begin{split}
& \mathcal{L}(x, \theta) =
\\
& \frac{1}{N}\sum_{i=1}^N \left\{ r_i(x; \theta) + \lambda \left\{ \exp(w_i(x; \theta)) \tilde{r}_i(x; \theta) - w_i(x; \theta) \right\} \right\}.
  \label{eq:ours_laplace}
  \end{split}
\end{equation}
Here, as in Eq.~\ref{eq:kendall_laplace}, we use $L1$ loss and the log likelihood of a Laplace distribution. $\tilde{r}_i$ is a virtual residual which we explain later. The certainty estimation $w_i$ depends on $\tilde{r}_i$ in Eq.9 while it depends on the target residual $r_i$ in Eq.3. By treating the virtual residual $\tilde{r}_i$ as a constant value, $L_t$ works only for the target optimization, and $L_u$ works only for the uncertainty optimization. The virtual residual $\tilde{r}_i$ enables us not only to adjust the balance of the two terms, but also to mitigate the effect of overfitting.

\begin{figure}[!t]
	\begin{center}
  \begin{subfloat}[{\bf}]{
          \includegraphics[clip, width=0.97\columnwidth]{./images/proposal_a_yu.png}
          \label{fig:proposal_a}}
        \end{subfloat}
        \\
\begin{subfloat}[{\bf}]{
          \includegraphics[clip, width=0.97\columnwidth]{./images/proposal_b_yu.png}
          \label{fig:proposal_b}}
        \end{subfloat}
	\end{center}

	\caption{{\bf An overview of the use of the virtual residual in our method}, using an example of age estimation. The training data is separated into $m$ subsets $(X_1, Y_1), (X_2, Y_2), ...... , (X_m, Y_m)$. Similarly, we define $m$ models $M_1, M_2, ...... , M_m$. We train each model $M_j$ using the data $\cup\left\{(X_k, Y_k) | k \neq j \right\}$ as shown in the top of Figure~\ref{fig:proposal_a}. $M_j$ is used for calculating the virtual residuals $\tilde{R}_j$ corresponding to the dataset $(X_j, Y_j)$. The virtual residuals $\tilde{R}_1, \tilde{R}_2, ...... , \tilde{R}_m$ are then used when training the model $M$ (which predicts the age and certainty) as shown in Figure~\ref{fig:proposal_b}.
}
	\label{fig:proposal}
\end{figure}

\subsection{Virtual Residuals}
\label{subsec:virtual_residual}

Next, we discuss the definition of the virtual residuals $\tilde{r}_i$. Fig.~\ref{fig:proposal} illustrates our method. We produce and save the virtual residuals beforehand in the manner described below.

First, in addition to the model $M$ for the target and uncertainty estimation, we define $m$ models $M_1, M_2, ...... , M_m$. Here, $M_j (1 \leq j \leq m)$ has the same construction of layers as $M$, except that $M_j$ outputs only the signal estimation, not the uncertainty estimation. We then randomly separate the training data $D$ into $m$ subsets, namely, $D_1, D_2, ...... , D_m$. Next, we train each model $M_j$ using the data $\cup\left\{D_k | k \neq j \right\}$. Now, we have $m$ trained models, and we let $M_j$ estimate the target values of data subset $D_j$. We then save the error of the signal estimation using the signal estimation from $M_j$. After predicting all $D_j$ by $M_j$, we compute the error $\tilde{r}_i$ for all the training data. Finally, we use $\tilde{r}_i$ in the training by applying Eq.~\ref{eq:ours_laplace}, and train the model $M$ by using all the training data.

We use this virtual residual $\tilde{r}_i$ because it helps us to predict the error, removing the effect of overfitting. When the model estimation of the signal $y_i$ overfits the training data, the residual $r_i$ of training data becomes smaller than that of test data. Therefore if we use Eq.~\ref{eq:kendall_laplace} for training, $w_i$ shows a high certainty, resulting in inaccurate error estimation. In our method, when we optimize the certainty estimation with a sample in $D_j$ for training, we use $\tilde{r}_i$, which is a residual produced from the target prediction of $M_j$. The samples in $D_j$ is not used when we train $M_j$, and thus, the predicted value of $M_j$ is not affected by overfitting to the samples. We can then avoid the certainty estimation being large due to overfitting. Therefore, the certainty estimation $w_i$ of the model $M$ is trained without the effect of overfitting.

\section{Experiments}
\label{sec:experiments}

We applied our method and the existing one to three different tasks. First, we trained the models on toy datasets to demonstrate specific features of the models. Second, we applied the uncertainty estimation to age estimation on the UTKFace dataset. Finally, we estimated uncertainty on an RGB-based depth prediction task, which is inherently difficult. In each of these experiments, we derived the loss function using a Laplacian prior, i.e., Eq.~\ref{eq:kendall_laplace} and Eq.~\ref{eq:ours_laplace}. When $\exp(w_i)$ in Eq.~\ref{eq:kendall_laplace} or Eq.~\ref{eq:ours_laplace} gets too large, we suppress it by replacing the term to $\exp(\frac{w_i}{c})$ with a constant value $c$. We introduce the experimental settings and their results of the three tasks respectively. We find that our method is successful in estimating the error and avoiding underestimation. We also discuss the influence of the virtual residuals by comparing our method with a separable formulation without virtual residuals.

\begin{figure}[!t]
\begin{center}
    \includegraphics[width=0.97\columnwidth]{./images/pir_definition.pdf}
    \caption{{\bf Visualization of the PiR($\eta$)}. When plotting $(r_i, \exp(-w_i))$ for each sample, PiR($\eta$) shows the rate of the sample that is included in the blue area.}
    \label{fig:pir_definition}
  \end{center}
\end{figure}

\begin{table}[!t]
\small
  \begin{center}
     \caption{{\bf Simulation results}. Showing performance according to metrics (RMSE, eRMSE, and PiR) for simple $L1$ loss estimation, the joint formulation (JF), and our model for various values of $\Delta$ in Eq.~\ref{eq:simulation_function_right}.}
      \begin{tabular}{|c|c||c||c|c|} \hline
    $\Delta$ & Loss & Target & eRMSE & PiR ($0.5$) \\
    & & RMSE ($\downarrow$) & ($\downarrow$) & ($\uparrow$) \\ \hline \hline
  \multirow{4}{*}{$0$}
  & $L1$ & 1.39 & - & -  \\ \cline{2-5}
  & JF~\cite{NIPS2017_7141} & 1.40 & 0.905 & 0.498 \\ \cline{2-5}
  & our method ($\lambda = 0.1$) & 1.40 & {\bf 0.857} & {\bf 0.554} \\ \hline
 \multirow{4}{*}{$1$}
  & $L1$ & 1.42 & - & - \\ \cline{2-5}
  & JF & 1.39 & 0.890 & 0.504 \\ \cline{2-5}
  & our method ($\lambda = 0.02$) & 1.42  & {\bf 0.858} & {\bf 0.550} \\ \hline
 \multirow{4}{*}{$5$}
  & $L1$ & 1.40 & - & - \\ \cline{2-5}
  & JF & 1.43 & 0.930 & 0.474 \\ \cline{2-5}
  & our method ($\lambda = 0.1$) & 1.41 & {\bf 0.853} & {\bf 0.566} \\ \hline
  \end{tabular}
\label{tbl:simulation_results}
  \end{center}
\end{table}

\begin{figure}[!t]
	\begin{center}
  \begin{subfloat}[{\bf $\Delta = 0$.}]{
          \includegraphics[clip, width=0.47\columnwidth]{images/pro_pro_pir.png}}
        \end{subfloat}
\begin{subfloat}[{\bf $\Delta = 1$.}]{
          \includegraphics[clip, width=0.47\columnwidth]{images/sta_pro_pir.png}}
        \end{subfloat}
\begin{subfloat}[{\bf $\Delta = 5$.}]{
          \includegraphics[clip, width=0.47\columnwidth]{images/hyp_pro_pir.png}}
        \end{subfloat}
	\end{center}

	\caption{{\bf PiR of the simulation.}
Each image shows the value of PiR($\eta$) as $\eta$ decreases by $0.1$ from $1$ to $0$. The green and red curves are the results of the joint formulation ~\cite{NIPS2017_7141} and our method, respectively.
}
	\label{fig:simulation_pir}
\end{figure}

\begin{figure*}[!ht]
	\begin{center}
  \begin{subfloat}[{The signal estimatino with JF while $\Delta = 0$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/pro_pro_kendall_t.png}}
        \end{subfloat}
\begin{subfloat}[{The signal estimation with JF while $\Delta = 1$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/sta_pro_kendall_t.png}}
        \end{subfloat}
\begin{subfloat}[{The signal estimation with JF while $\Delta = 5$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/hyp_pro_kendall_t.png}}
        \end{subfloat}
        \begin{subfloat}[{The error estimation with JF while $\Delta = 0$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/pro_pro_kendall_u.png}}
        \end{subfloat}
\begin{subfloat}[{The error estimation with JF while $\Delta = 1$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/sta_pro_kendall_u.png}}
        \end{subfloat}
\begin{subfloat}[{The error estimation with JF while $\Delta = 5$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/hyp_pro_kendall_u.png}}
        \end{subfloat}
         \begin{subfloat}[The signal estimation with our method while $\Delta = 0$.]{
        \includegraphics[clip, width=0.31\columnwidth]{images/pro_pro_ours_t.png}}
        \end{subfloat}
\begin{subfloat}[{The signal estimation with our method while $\Delta = 1$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/sta_pro_ours_t.png}}
        \end{subfloat}
\begin{subfloat}[{The signal estimation with our method while $\Delta = 5$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/hyp_pro_ours_t.png}}
        \end{subfloat}
        \begin{subfloat}[{The error estimation with our method while $\Delta = 0$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/pro_pro_ours_u.png}}
        \end{subfloat}
\begin{subfloat}[{The error estimation with our method while $\Delta = 1$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/sta_pro_ours_u.png}}
        \end{subfloat}
\begin{subfloat}[{The error estimation with our method while $\Delta = 5$.}]{
          \includegraphics[clip, width=0.31\columnwidth]{images/hyp_pro_ours_u.png}}
        \end{subfloat}
	\end{center}
  \caption{{\bf Signal and uncertainty predictions during the simulation.}  The upper figures show the joint formulation (JF) results~\cite{NIPS2017_7141}, while the lower figures show those of our method. In the six figures on the left side ((a), (b), (c), (g), (h), and (i)), the $x$-axis corresponds to the input, and the $y$-axis to the signal. Blue points indicate the training data; red points are the models' predictions for the test data. In the other six figures on the right side, the $x$-axis corresponds to the input, and the $y$-axis corresponds to the error of the signal estimation. Blue points indicate the ideal mean error given by $\sigma(x)$, and red points are estimations of the error by the models calculated by Eq.~\ref{eq:uncertainty_optimization}. Each point corresponds to one test datum. }
  \label{fig:simulation_visualization}
\end{figure*}

\subsection{Metrics for Uncertainty Estimation}
\label{subsec:metrics}

First of all, we introduce the metrics that we use in our experiments. Sparsification plots are widely used for accurate evaluation of the performance of uncertainty estimation~\cite{wannenwetsch2017probflow,bruhn2006confidence,mac2012learning,eddy2018opticalflow}. In these plots, samples are removed in order of the uncertainty. The plots themselves do not show the uncertainty estimation performance numerically. In this work, we focus on how accurately the model can estimate the value of the error. We evaluate the performance of the error estimation using the following two metrics:
\begin{itemize}
    \item {\bf eRMSE}: To evaluate the goodness of the certainty value $w_i$ based on Eq.~\ref{eq:uncertainty_optimization},
we introduce a metric eRMSE defined as:
\begin{equation}
  eRMSE = \sqrt{E\left[| \exp(-w_i(x)) - r_i(x) |^2\right]}.
  \label{eq:ermse_definition}
\end{equation}
This metric indicates how well the uncertainty is estimated.
\item {\bf PiR}: We check whether each sample satisfies
\begin{equation}
  \eta \exp(-w_i(x)) \leq r_i(x) < \frac{1}{\eta} \exp(-w_i(x)),
  \label{eq:pir_definition}
\end{equation}
where $\eta$ is a fixed value such that $0 < \eta \leq 1$. We define the rate of the samples satisfying Eq.~\ref{eq:pir_definition} as PiR($\eta$). Considering Eq.~\ref{eq:uncertainty_optimization}, when we plot a point $(r_i, \exp(-w_i))$ for each sample on a two-dimensional plane, ideally the points are located near the line $\exp(-w) = r$. PiR($\eta$) indicates the rate of the points that is located between two lines $\frac{1}{\eta} \exp(-w) = r$ and $\eta \exp(-w) = r$ as Figure~\ref{fig:pir_definition} shows, and thus, it works as another metric indicating how well the uncertainty is estimated.
\end{itemize}
In addition, we compare the performance of target estimation using the root mean square error (RMSE).

\begin{table}[!t]
\small
  \begin{center}
      \caption{{\bf Age estimation results on the UTKFace dataset}.}
\label{tbl:age_results}
      \begin{tabular}{|c||c||c|c|} \hline
    Loss & Target & eRMSE & PiR($\eta$) \\
    & RMSE ($\downarrow$) & ($\downarrow$) & ($\uparrow$) \\ \hline \hline
    $L1$  & 7.37 & - & - \\ \hline
    JF~\cite{NIPS2017_7141} & 7.51 & 6.10 & 0.352 \\ \hline
    our method ($\lambda = 0.02$) & 7.27 & {\bf 4.99} & {\bf 0.523} \\ \hline
  \end{tabular}
  \end{center}
\end{table}

\begin{figure}[!t]
\begin{center}
    \includegraphics[width=0.5\columnwidth]{./images/age_pir.png}
    \caption{{\bf PiR of the age estimation.}
This plot shows the value of PiR($\eta$), where $\eta$ decreases by $0.1$ from $1$ to $0$. The green and red line illustrate the results of the joint formulation~\cite{NIPS2017_7141} and our method, respectively.
}
    \label{fig:age_pir}
  \end{center}
\end{figure}

\subsection{Simulation}
\label{subsec:simulation}

In this experiment, we perform a one-dimensional regression task to clearly visualize the behavior of the signal and uncertainty estimations.

\paragraph{Datasets}

We randomly generate input $\{x | 0 \leq x < 1\}$ from a uniform distribution.
The target $y(x)$ is generated from the signal $f(x)$ and noise $n(x)$ as below:
\begin{equation}
  y(x) = f(x) + n(x).
  \label{eq:simulation_target_definition}
\end{equation}
$f(x)$ is defined as
\begin{align}
  f(x) &= 2x - 1 (0 \leq x < 0.5), \label{eq:simulation_function_left} \\
  f(x) &= 2x - 1 + \Delta (0.5 \leq x < 1).  \label{eq:simulation_function_right}
\end{align}
We choose three different $\Delta$ values: $0$, $1$, and $5$.
Note that $f(x)$ is a linear function without a step when we choose $0$ for $\Delta$.

\begin{table}[!t]
\begin{center}
  \caption{{\bf Examples of age and error estimation by two methods (joint formulation and ours)}.}
  \label{tbl:age_examples}
  \begin{tabular}{|c|c||c||c|c|c|}
    \hline input & age & Loss & age & prediction & error \\
    & & & estimation & error & estimation \\ \hline
    \hline
  \multirow{2}{*}{
\begin{minipage}{8mm}
      \centering
      \scalebox{0.08}{\includegraphics{./images/utkface_2858.jpg}}
    \end{minipage}
  }
  &
  \multirow{2}{*}{20}
	& JF~\cite{NIPS2017_7141} & 23.131 & 3.131 & 1.499 \\ \cline{3-6}
  & & our method & 22.692 & 2.692 & 2.655 \\ \hline \hline
  \multirow{2}{*}{
\begin{minipage}{8mm}
      \centering
      \scalebox{0.08}{\includegraphics{./images/utkface_21579.jpg}}
    \end{minipage}
  }
  &
  \multirow{2}{*}{68}
	& JF & 59.401 & 8.599 & 2.487 \\ \cline{3-6}
  & & our method & 58.005 & 9.995 & 6.553 \\ \hline \hline
  \end{tabular}
  \end{center}
\end{table}

Noise $n(x)$ is dependent on the input and generated by the Gaussian function $\mathcal{N}(0, \sigma(x)^2)$.
For $\sigma(x)$, we choose a linear function that has a large value when $x$ is large:
\begin{equation}
  \sigma(x) = 1.99x + 0.01.
  \label{eq:simulation_variance}
\end{equation}
In our experiments, given a noisy target $y_i(x)$, we estimate the signal $f_i(x)$ and certainty $w_i(x)$; the latter is taken to be $-\log(|n_i(x)|)$.

We generate 128 samples for training and 500 samples for testing. We train models with relatively small datasets to simulate a finite dataset for large multidimensional tasks and use more test samples for better visualization of the results.

\begin{figure}[!t]
	\begin{center}
  \begin{subfloat}[The age estimation with JF.]{
          \includegraphics[clip, width=0.47\columnwidth]{images/scatter_kendall_t.png}}
        \end{subfloat}
\begin{subfloat}[The error estimation with JF.]{
          \includegraphics[clip, width=0.47\columnwidth]{images/scatter_kendall_u.png}}
        \end{subfloat}
\begin{subfloat}[The age estimation with our method.]{
          \includegraphics[clip, width=0.47\columnwidth]{images/scatter_ours_t.png}}
        \end{subfloat}
        \begin{subfloat}[The error estimation with our method.]{
          \includegraphics[clip, width=0.47\columnwidth]{images/scatter_ours_u.png}}
        \end{subfloat}
	\end{center}

  \caption{{\bf Model outputs for each image in the age estimation experiment.} The upper figures show the results of the joint formulation~\cite{NIPS2017_7141}, and the lower figures show the results of our method. The figures on the left side show the ground truth ($x$-axis) and predicted age ($y$-axis) of each sample in the test data.
The figures on the right side show the error of the age prediction ($x$-axis) and the prediction of the error ($y$-axis) of each sample.
}
\label{fig:age_scatter}
\end{figure}

\paragraph{Implementation Details}

We use a simple network composed of eight fully connected layers. We first feed the input into the layers, which have outputs of $32, 64, 128, 128, 64$, and $32$, respectively. After the last layer follows a pair of layers having outputs of $16$ and $1$; these predict the signal and uncertainty. Each layer except for the last two is followed by a ReLU. We divide the training data into 10 subsets for the training of each $M_j$.

\paragraph{Results}

Table~\ref{tbl:simulation_results} compares the values of the estimation metrics, and Figure~\ref{fig:simulation_pir} compares the PiR of the methods.
From Table~\ref{tbl:simulation_results}, it can be seen that our method has the best performance on the two metrics of uncertainty estimation. When $\Delta$ increases, which means that signal estimation becomes more difficult, the joint formulation becomes worse for all metrics. By contrast, our method seems to predict both the signal and uncertainty almost independently of the value of $\Delta$. Figure~\ref{fig:simulation_pir} compares the PiR($\eta$) of the frameworks, and it also shows that our method outperforms the others across a wide range of $\eta$.

Figure~\ref{fig:simulation_visualization} compares ourmethod to the joint formulation. When we look at the signal estimation, we see that the larger $\Delta$ is, the more unstable the signal estimation becomes in all the methods, and also that the task becomes difficult when $\Delta$ is large. By contrast, when we focus on the uncertainty estimation in Figure~\ref{fig:simulation_visualization}, the joint formulation is unstable with a large $\Delta$. We think that this instability is a result of being unable to balance the weight between the target and uncertainty estimation in the joint estimation. In the case of our method, the error estimation is stable and closer to the ideal mean error given by Eq.~\ref{eq:simulation_variance}. The error estimation is larger when $x$ is larger, but this is because the target estimation is inappropriate for this noisy area.

\subsection{Age Estimation}
\label{subsec:age_estimation}

In this experiment, we predict ages and their uncertainty from facial images.

\paragraph{Datasets}

We use the UTKFace dataset~\cite{zhifei2017cvpr} for the age estimation. The UTKFace dataset is a large-scale facial dataset consisting of more than 20,000 data. It is composed of cropped images of the faces of people ranging in age from 0 to 116 years. They used DEX~\cite{Rothe2016IJCV} algorithm to guess the ground truth of the ages, and then checked it manually. We randomly choose 80\% of the data for training, 10\% for validation, and the remaining 10\% for testing. For the inputs of the networks, we resize the images into a resolution of $64 \times 64$.

\begin{table}[!t]
\small
  \begin{center}
  \caption{{\bf Depth estimation results on NYU Depth Dataset V2}. All metrics were calculated for each image; the table shows averages over all images.}
  \label{tbl:depth_results}
  \begin{tabular}{|c||c||c|c|} \hline
    loss & Target & eRMSE & PiR($0.5$) \\
    & RMSE ($\downarrow$) & ($\downarrow$) & ($\uparrow$) \\ \hline \hline
	$L1$ & 0.529 & - & - \\ \hline
    JF~\cite{NIPS2017_7141}  & 0.531 & 0.482  & 0.364 \\ \hline
    our method ($\lambda = 0.05$)  & 0.522 & {\bf 0.399} & {\bf 0.401} \\ \hline
  \end{tabular}
  \end{center}
\end{table}

\begin{figure}[!t]
\begin{center}
    \begin{subfloat}{
          \includegraphics[clip, width=0.47\columnwidth]{images/depth_pir.png}}
        \end{subfloat}
\begin{subfloat}{
          \includegraphics[clip, width=0.47\columnwidth]{images/depth_pir_magnify.png}}
        \end{subfloat}
    \caption{{\bf PiR of the depth estimation.}
The left image shows the value of PiR($\eta$), where $\eta$ decreases by $0.1$ from $1$ to $0$. The enlarged detail on the right shows the region where $\eta$ decreases from $0.7$ to $0.3$. The green and red curves are the results of the joint formulation~\cite{NIPS2017_7141} and our method, respectively.
}
    \label{fig:depth_pir}
  \end{center}
\end{figure}

\paragraph{Implementation Details}

We use ResNet50~\cite{he2016deep}, which is pre-trained on ImageNet~\cite{deng2009imagenet}. This pre-trained model is provided on PyTorch. We change the last layer of the provided model to regress two values: the age and the certainty of the estimation. During this experiment, we divide the training data into 20 subsets for the training of each $M_j$.

\begin{figure*}[!t]
	\begin{center}
  \begin{subfloat}[Input]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_input.jpg}}
    \end{subfloat}
\begin{subfloat}[GT]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_gt.jpg}}
    \end{subfloat}
    \begin{subfloat}[Pred. (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_kendall_d.jpg}}
        \end{subfloat}
        \begin{subfloat}[Error (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_kendall_e.jpg}}
        \end{subfloat}
        \begin{subfloat}[Unc. (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_kendall_u.jpg}}
        \end{subfloat}
\begin{subfloat}[Pred. (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_ours_d.jpg}}
    \end{subfloat}
\begin{subfloat}[Error (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_ours_e.jpg}}
        \end{subfloat}
\begin{subfloat}[Unc. (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_5_ours_u.jpg}}
        \end{subfloat}
        \vspace{2mm}
\begin{subfloat}[Input]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_input.jpg}}
    \end{subfloat}
\begin{subfloat}[GT]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_gt.jpg}}
    \end{subfloat}
    \begin{subfloat}[Pred. (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_kendall_d.jpg}}
        \end{subfloat}
        \begin{subfloat}[Error (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_kendall_e.jpg}}
        \end{subfloat}
        \begin{subfloat}[Unc. (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_kendall_u.jpg}}
        \end{subfloat}
\begin{subfloat}[Pred. (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_ours_d.jpg}}
    \end{subfloat}
\begin{subfloat}[Error (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_ours_e.jpg}}
        \end{subfloat}
\begin{subfloat}[Unc. (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_16_ours_u.jpg}}
        \end{subfloat}
        \vspace{2mm}
  \begin{subfloat}[Input]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_input.jpg}}
    \end{subfloat}
\begin{subfloat}[GT]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_gt.jpg}}
    \end{subfloat}
    \begin{subfloat}[Pred. (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_kendall_d.jpg}}
        \end{subfloat}
        \begin{subfloat}[Error (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_kendall_e.jpg}}
        \end{subfloat}
        \begin{subfloat}[Unc. (JF)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_kendall_u.jpg}}
        \end{subfloat}
\begin{subfloat}[Pred. (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_ours_d.jpg}}
    \end{subfloat}
\begin{subfloat}[Error (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_ours_e.jpg}}
        \end{subfloat}
\begin{subfloat}[Unc. (ours)]{
          \includegraphics[clip, width=0.23\columnwidth]{images/nyu_66_ours_u.jpg}}
        \end{subfloat}
	\end{center}

	\caption{{\bf Examples of depth estimation.}
From left to right: input RGB images (Input), ground-truth (GT) depth maps, joint formulation (JF)~\cite{NIPS2017_7141} predicted depth map (Pred.), JF absolute errors (Error), JF expected uncertainty map (Unc.), and the three corresponding images from our method. In the ground truth and predicted depth maps, the black or red regions have large depth values, whereas the blue or green regions contain small depth values. In the error and uncertainty maps, the black or red regions have larger values, whereas the white or yellow regions contain smaller values.
}
	\label{fig:depth_examples}
\end{figure*}

\begin{figure}[!t]
\begin{center}
    \includegraphics[width=0.47\columnwidth]{./images/scatter_asais_u.png}
    \caption{{\bf Plot of model outputs for each image in the age estimation}. The error of the age prediction of each sample is shown on the $x$-axis, and the prediction of the error on the $y$-axis. The model is trained with the separate formulation without virtual residuals.}
    \label{fig:age_estimation_ablation}
  \end{center}
\end{figure}

\paragraph{Results}

Table~\ref{tbl:age_results} compares the age estimation in the two frameworks of uncertainty estimation, and in a simple $L1$ loss estimation without uncertainty estimation. Figure~\ref{fig:age_pir} compares the PiR($\eta$) of the frameworks. They show that our method outperforms the existing methods in terms of eRMSE and PiR.

To analyze the differences, we plot two sets of pairs in Figure~\ref{fig:age_scatter}: ground truth and predicted ages, and actual and estimated prediction errors. A comparison of the two figures for the age estimation shows that the models estimate the age similarly; both approximately follow the line $y = x$. By contrast, the figures for the error estimation are different. The existing method, JF, uniformly predicts the error to be small, even when the real error is large, whereas our method predicts a larger error in response to a larger real estimation error. The models (including ours) always estimate the error as less than 30 years. This is probably because the number of images that have such large errors is few.

Table~\ref{tbl:age_examples} shows two examples of age estimation. The existing method, JF, tends to underestimate the error. When the error is larger, this tendency is clearer. The underestimation of large errors is especially mitigated by our method.

\subsection{Depth Estimation}
\label{subsec:depth_estimation}

Finally, we compare the models based on a depth estimation task using RGB images. This task is inherently difficult~\cite{Ma2018SparseToDense}.

\paragraph{Datasets}

We use a popular RGB-based depth prediction task, NYU Depth V2~\cite{Silberman:ECCV12}.
This dataset consists of RGB and depth images collected from 464 different indoor scenes using Microsoft Kinect. We use the official split of the data: 249 scenes are used for training, and the remaining 215 are used for testing, as in~\cite{Ma2018SparseToDense,laina2016deeper,eigen2014depth}.

\paragraph{Implementation Details}

We implement our proposed method and the existing ones at the top of the official sparse-to-dense~\cite{Ma2018SparseToDense} PyTorch implementation\footnote{\url{github.com/fangchangma/sparse-to-dense.pytorch}}. Note that we do not use any sparse depth sample as an input. We increase the number of epochs to 30, because it is difficult in this task for the models to achieve convergence. We implement the network layers in the same way as ~\cite{asai}, using ResNet~\cite{he2016deep} and UpProj~\cite{laina2016deeper}.
(For more details, please refer to ~\cite{asai}.)
During this difficult task, we divide the training data into 30 subsets for the training of each $M_j$.

\paragraph{Results}

In this experiment, the depth and uncertainty are estimated per pixel, and thus we calculate all metrics per image. Table~\ref{tbl:depth_results} compares the depth estimation of the two uncertainty estimation frameworks and a simple $L1$ loss estimation. Figure~\ref{fig:depth_pir} compares the PiR($\eta$) of the two frameworks. The values of eRMSE and PiR($\eta$) in our framework exceed those in JF. This task is difficult for neural networks to predict ~\cite{Ma2018SparseToDense}, and thus the epistemic uncertainty of the networks is larger. In such a case, quantifying the error from the aleatoric uncertainty is much more difficult, resulting in smallness of the difference of performance between the methods.

Figure~\ref{fig:depth_examples} shows examples of the estimations. In the JF estimation, there are areas where the error estimation is inappropriately small, and there are some white areas in where the errors are large. Additionally, JF is less sensitive to the wall, where the errors tend to be larger. Our method avoids these phenomena.

\subsection{The benefit of the virtual residuals}

Lastly, we argue for the benefit of introducing virtual residuals to optimize $w_i$.

We conduct experiments without using virtual residuals -- we substitute the value $r_i$ directly for $\tilde{r_i}$ in Eq.~\ref{eq:ours_laplace} instead of virtual residuals. We treat $\tilde{r}_i$ as a constant in a minibatch and do not calculate the gradient because we want to balance the weight of the two terms, $L_t$ and $L_u$, and focus on the ablation of the virtual residuals.

Figure~\ref{fig:age_estimation_ablation} shows the prediction error and its estimation in the age estimation case. Comparing this with Figure~\ref{fig:age_scatter}, we can see that the separable formulation without virtual residuals estimates the error to be much smaller than the joint formulation. We consider this to be robustness of the target estimation in Eq.~\ref{eq:kendall_laplace}. The partial derivative of Eq.~\ref{eq:kendall_laplace} with respect to $\theta$ is
\begin{equation}
  \frac{\partial\mathcal{L}}{\partial\theta} = \frac{1}{N}\sum_{i=1}^N \left\{ \exp(w_i)\frac{\partial r_i}{\partial \theta} + \exp(w_i) r_i \frac{\partial w_i}{\partial \theta} - \frac{\partial w_i}{\partial \theta} \right\}.
  \label{eq:kendall_laplace_derivative}
\end{equation}
When given noisy samples with large aleatoric uncertainties, $w_i$ becomes small in the training. Therefore the first term of Eq.~\ref{eq:kendall_laplace_derivative} also becomes small. It implies that the signal estimation by Eq.~\ref{eq:kendall_laplace} is robust. (This effect is also mentioned in ~\cite{NIPS2017_7141} as attenuating loss.) However, this effect is not sufficient as we can see in the experiments. The virtual residuals avoid underestimating.

\section{Conclusion}
\label{sec:conclusion}

In this study, we presented a framework for error quantification based on  heteroscedastic aleatoric uncertainty estimation in regression problems of deep learning. Existing methods cannot estimate the error accurately because of overfitting in the target estimation. We mitigated the effect of overfitting with virtual residuals and avoided underestimating the estimation errors. In the experiments, we evaluated our method on three different tasks, and showed that it estimated both targets and errors, outperforming the existing method on many metrics.

\section*{Acknowledgment}
This work was partially supported by JST JPMJCR19F4 and JSPS 18H03254.

\bibliographystyle{IEEEtran}

\input{arxiv.bbl}

\end{document}