11institutetext: Aya Bamba 22institutetext: The University of Tokyo, Hongo
7-3-1, Bunkyo-ku, Tokyo, Japan 113-0033, 22email<EMAIL_ADDRESS>33institutetext: Brian J. Williams 44institutetext: NASA Goddard Space Flight
Center, X-ray Astrophysics Laboratory, Greenbelt, MD 20771 44email:
<EMAIL_ADDRESS>
# Supernova remnants: Types and evolution

Aya Bamba and Brian J. Williams corresponding author

###### Abstract

Although only a small fraction of stars end their lives as supernovae, all
supernovae leave behind a supernova remnant (SNR), an expanding shock wave
that interacts with the surrounding medium, heating the gas and seeding the
cosmos with elements forged in the progenitor In this chapter, we introduce
the basic properties of galactic and extragalactic SNRs (§2). We summarize how
SNRs evolve throughout their life cycles over the course of $\sim 10^{6}$ yrs
(§3). We discuss the various morphological types of SNRs and discuss the
emission processes at various wavelengths.(§4).

##  1 Keywords

Supernova remnants (SNRs) — free expansion phase — adiabatic expansion phase —
shell-type SNR — plerion-type SNR — mixed morphology SNR

##  2 Introduction

Supernova (SN) explosions are among the most energetic events in the universe
since the Big Bang, releasing more energy ($\sim 10^{51}-10^{53}$ ergs) than
the Sun will release over its entire lifetime. They are the cataclysmic ends
of certain types of stars, and are responsible for seeding the universe with
the material necessary to form other stars, planets, and life itself. We owe
our very existence to generations of stars that lived and died billions of
years ago, before the formation of our Sun and solar system. The processes of
stellar evolution continue to occur today, with typical galaxies like the
Milky Way hosting a few supernovae per century, on average. The study of
supernovae and the role they play in shaping the evolution of star systems and
galaxies is truly an exploration of our own origins. Supernova remnants
(SNRs), the expanding clouds of material that remain after the explosion,
spread elements over volumes of thousands of cubic light-years, and heat the
interstellar medium through fast shock waves generated by the ejecta from the
star.

There are currently $\sim 300$ SNRs cataloged in our Galaxy (Green, 2019),
although radio surface brightness studies predict that only half of SNRs have
been identified (Case & Bhattacharya, 1998). Recent X-ray observations with
good spatial resolution and large effective area allow us to detect many SNR
samples in nearby galaxies as well, such as Magellanic CloudsMaggi et al.
(2016), M31Sasaki et al. (2012), and M33Garofali et al. (2017).

##  3 Evolution of supernova remnants

The general picture of an SNR is that of a cloud of expanding material ejected
from the star, rich in heavy elements like O, Si, and Fe, that drive a shock
wave into the interstellar medium (ISM). This shock waves heat the
interstellar gas it encounters to millions of degrees, creating a highly-
ionized plasma. Figure 1 shows Cassiopeia A in X-rays, an example of a young
($\sim 350$ yr) SNR in our own galaxy. Although Cassiopeia A is known to have
resulted from a core-collapse (CC) SN (see Chapter on SN), it is generally
difficult to tell the type of SN only by looking at the remnant. SNRs remain
visible for thousands, often tens or hundreds of thousands of years before
dissipating their energy into the ISM. The life of a SNR can be thought of as
consisting of four phases.

Figure 1: Cassiopeia A, a young remnant in our Galaxy, as seen at X-ray
energies with Chandra. Credit: NASA/CXC/SAO

###  3.1 Free expansion phase

The ejected gas mainly determines the evolution of SNRs in the initial phase.
Since the pressure of surrounding matter is negligible in comparison, the
exploded ejecta expands effectively without deceleration. As a result, this
phase is called the “free expansion” phase.

The characteristic (or typical) speed of the ejecta expansion, $v_{s}$, is
given by

$v_{s}=\sqrt{\frac{2E}{M_{ej}}}\ \ ,$ (1)

where $E$ and $M_{ej}$ are the kinetic energy of the explosion and the ejecta
mass, respectively. With fiducial values of $E_{kin}=1\times 10^{51}$ erg and
$M_{ej}=1M_{\odot}$, a typical velocity of the ejecta is $\sim 10^{4}$ km s-1.
While this velocity for the ejecta is quite high, it is also a bulk velocity
with very little randomization in the velocity of the various fluid elements
without deceleration (This is the difinition of “free expansion”). Thus, there
is no heating and the ejecta temperature is actually quite cold in the
beginning, and the remnant is generally faint at all wavelengths. The
exception is the radiation from nonthermal synchrotron emission from high
energy electrons accelerated in the shock (see Chapter “acceleration”).

Virtually all SNRs emit in radio waves, where the emission is entirely
nonthermal in origin, resulting from synchrotron radiation from relativistic
electrons spiraling around magnetic fields in the shock. Synchrotron emission
is characterized by a featureless power-law spectrum, where the radio flux,
$S_{\nu}$ is given by $S_{\nu}\propto\nu^{-\alpha}$, where $\nu$ is the
frequency and $\alpha$ is the spectral index, which depends on the energy
distribution of the electron population. For more detail, see Chapter
“acceleration”. In Figure 2, we show an example of a remnant, G1.9$+$0.3,
which is still in the free expansion phase Reynolds et al. (2008a). The
acceleration of electrons up to GeV energies is sufficient to produce radio
synchrotron emission.

Figure 2: Left: SNR G1.9+0.3, the youngest remnant in the Galaxy, seen in
radio waves from the VLA in 1985. Right: The remnant, in 2007, seen in X-rays
with Chandra. In the 22 years between these images, the remnant expanded
substantially, allowing the time since the original supernova explosion (about
140 years) to be estimated. Credit: X-ray (NASA/CXC/NCSU/S.Reynolds et al.);
Radio (NSF/NRAO/VLA/Cambridge/D.Green et al.).

A small population of SNRs are apparently able to accelerate electrons to much
higher energies ($>10$ TeV), sufficient to produce synchrotron radiation at
X-ray energies Koyama et al. (1995); Aharonian et al. (2004). G1.9+0.3 also
falls into this category, and Figure 2 also shows an image of the remnant that
is essentially entirely from nonthermal synchrotron emission. For more details
regarding the particle acceleration mechanisms, see the review by Reynolds et
al. (2008b).

###  3.2 Adiabatic expansion phase

At the end of the free expansion phase, the expansion begins to decelerate due
to the forward shock sweeping up a non-negligible amount of material in the
ISM. An estimate for when this transition occurs can be obtained by
calculating the radius, $R_{S}^{0}$, at which the swept-up mass equals the
ejecta mass:

$\displaystyle\frac{4}{3}\pi{R_{s}^{0}}^{3}\rho_{0}$ $\displaystyle=$
$\displaystyle M_{ej}\ \ ,$ (2) $\displaystyle R_{s}^{0}$ $\displaystyle=$
$\displaystyle
5.8{\left(\frac{M_{ej}}{10M_{\odot}}\right)}^{1/3}{\left(\frac{n_{0}}{0.5\
{\rm cm}^{-3}}\ \right)}^{-1/3}\ ({\rm pc})\ \ ,$ (3)

where $\rho_{0}$ and $n_{0}$ are mass and number densities in the cgs unit.
Assuming typical values of supernova explosion and ISM, it generally takes a
few hundred years for the radius of the SNR to reach $R_{s}^{0}$ (though this
can vary by orders of magnitude if the SN explodes into particularly low or
high density regions of the ISM).

This phase of deceleration of the forward shock wave also leads to the
formation of a “reverse shock,” which propagates back into the ejecta towards
the center of the SNR. Similarity solutions for the reverse shock are
discussed in Hamilton & Sarazin (1984a, b).

In this phase, the evolution of the SNR radius as a function of time can be
written as

$R(t)\propto t^{m}\ \ ,$ (4)

where $m$ is the expansion parameter defined by Hughes (1999). At the end of
the free expansion phase, $m\approx 1$, by definition, but during this second
phase, $m$ gradually decreases to $m=0.4$ Sedov (1959).

When the SNR shock starts decelerating, radiative cooling is still negligible
and the expansion is adiabatic. For the most simple case of a uniform ISM, the
evolution can be written as the results of the similarity solution by Sedov
(1959) (thus, this phase is often referred to as the “Sedov phase”).

Assuming that a large amount of energy is instantaneously released in a small
volume, expansion can be characterized by only two parameters, $\rho_{0}$ (the
mass density of the ISM) and $E_{kin}$ (the kinetic energy of the explosion).
We can define a non-dimensional parameter $\zeta$ as

$\zeta=R_{s}{\left(\frac{Et^{2}}{\rho_{0}}\right)}^{-1/5}\ \ ,$ (5)

where $R_{s}$ and $t$ represent the radius and age of the SNR. $R_{s}$ can be
written as

$R_{s}=\zeta_{0}{\left(\frac{Et^{2}}{\rho_{0}}\right)}^{1/5}\ \ ,$ (6)

where $\zeta_{0}$ is determined from the energy conservation equation. Using
$\zeta_{0}=1.17$ for ideal gas ($\gamma=5/3$) Landau & Lifshitz (1959),
$R_{s}$, the temperature $T$, and the shock velocity $v_{s}$ are written as

$\displaystyle R_{s}$ $\displaystyle=$ $\displaystyle
5.0{\left(\frac{E_{51}}{n_{0}}\right)}^{1/5}t_{3}{}^{2/5}\ \ {\rm(pc)}\ \ ,$
(7) $\displaystyle T_{s}$ $\displaystyle=$ $\displaystyle 1.5\times
10^{10}\left(\frac{E_{51}}{n_{0}}\right)R_{s}^{-3}\ \ {\rm(K)}$ (8)
$\displaystyle=$ $\displaystyle 1.2\times
10^{8}\times{\left(\frac{E_{51}}{n_{0}}\right)}^{2/5}t_{3}{}^{-6/5}\ \
{\rm(K)}\ \ ,$ (9) $\displaystyle v_{s}$ $\displaystyle=$
$\displaystyle\frac{dR_{s}}{dt}=2.1\times
10^{8}\times{\left(\frac{E_{51}}{n_{0}}\right)}^{1/5}t_{3}{}^{-3/5}\ \
{\rm(cm\ s^{-1})}\ \ ,$ (10)

where $E_{51}$, $n_{0}$, and $t_{3}$ are the explosion energy in the unit of
$10^{51}$ erg, the number density of the ambient gas in the unit of cm-3, and
the SNR age in the unit of $10^{3}$ years. During the adiabatic expansion
phase, $\sim 70\%$ of the initial explosion energy is converted into thermal
energy of swept-up matter Chevalier (1974).

Remnants in the adiabatic or Sedov expansion phase are typically bright at
X-ray wavelengths, owing to both the hot ISM gas shocked by the forward shock
and the hot ejecta shocked by the reverse shock. In Figure 3, we show a
canonical example of a remnant in this phase: DEM L71, a remnant in the LMC
believed to be a few thousand years post-explosion Hughes et al. (2002);
Rakowski et al. (2003).

Figure 3: SNR DEM L71, located in the LMC. This X-ray image shows a clear
separation between the forward shocked material in the ring on the outside and
the reverse-shocked material in blue on the inside. Credit:
NASA/CXC/Rutgers/J.Hughes et al.

High-mass stars typically eject mass as stellar winds before their explosions,
forming a circumstellar medium (CSM). The winds make a density gradient around
the progenitor star, making the deceleration of the shock slower. A modified
similarity solution for this case is shown in Truelove & McKee (1999).

These models introduced here ignored radiative cooling under the assumption of
an ideal gas. In the case that the system is not truly adiabatic, energy is
robbed from the shock and $v_{s}$ becomes smaller. Such a case can be caused
by efficient acceleration of cosmic rays at the shock front (Berezhko et al.,
2002, for example).

###  3.3 Snowplow phase

When the shock velocity has decelerated to $\sim$200 km s-1, the temperature
of the gas drops below $\sim 10^{6}$ K and radiative cooling begins to affect
the dynamics of the shock evolution. Cooling of the gas is a runaway process,
as the more it cools, the more the cooling rate increases. The density of the
shock front becomes highest where the timescale of radiative cooling becomes
shortest; as a result, a cool, dense shell is formed just behind the shock
front, whereas the inner gas still has high temperature and high pressure due
to the low density. It is during this phase that the remnant becomes bright in
optical emission from recombining atoms, most notably H$\alpha$, [S II], [N
II], and [O III]. As an example of a remnant in this phase, we show an HST
image of a portion of the Cygnus Loop in Figure 4. The Cygnus Loop is the
remains of a SN believed to be 10,000-20,000 yr post explosion.

Figure 4: A portion of the Cygnus Loop SNR, seen with HST. Narrowband images
of [O III], [S II], and H$\alpha$, plus V and I continuum bands combine to
make this image. Credit: NASA/ESA/Hubble Heritage Team.

The cool shell continues the expansion further collecting ISM gas in a manner
reminiscent of a snowplow (hence the name ”snowplow phase”). At first, the gas
expands adiabatically, then $pV^{5/3}$ becomes a constant. The shock expands
with the time dependency of $R_{s}\propto t^{2/7}$. As the temperature cools
down further, the pressure can be ignored and the cool shell expands at a
constant radial momentum with a time dependency of roughly $R_{s}\propto
t^{1/4}$.

###  3.4 Dissipation phase

As the SNR gets even older, the expansion velocity becomes smaller. When the
shock slows down to a velocity comparable to the sound velocity, the SNR
blends into the interstellar medium. This is called the “Dissipation phase,”
and marks the end of the SNR phase. The lifetime of typical SNRs is $\sim
10^{6}$ yrs, during which time they will expand to radii of dozens to hundreds
of parsecs, filling a volume large enough to encompass hundreds to thousands
of star systems.

##  4 Types of supernova remnants

###  4.1 Shell-type SNRs

The majority of remnants are classified as “shell-type” SNRs, where the
forward shock is delineated by emission at wavelengths across the electro-
magnetic spectrum. The remnants shown in Figures 1$-$4 are all shell-type
SNRs. The physical processes responsible for emitting at various wavelengths
differ. In radio regime, emission is dominated by synchrotron emission from
energetic electrons (accelerated by the shock wave to GeV energies) spiraling
around the turbulent magnetic field at the shock front. Historically, most
SNRs have been discovered and categorized based on their radio emission
(Green, 2019).

At infrared wavelengths, young SNRs (those in the free expansion or adiabatic
expansion phases) are generally dominated by thermal emission from warm dust
grains; see Dwek & Arendt (1992) for a review. Dust grains themselves do not
generally feel the passage of the shock wave, but are suddenly immersed in the
hot plasma behind the shock, where ion and electron temperatures can be tens
of millions of degrees. Grains are collisionally heated by these particles,
and reach temperatures of 50-200 K. Dust at these temperatures radiates in the
mid-IR, with thermal spectra peaking in the 20-100 $\mu$m range.

Generally speaking, IR emission from young shell-type SNRs only shows up in
emission from the shell of the remnant. This emission results from warmed
interstellar grains; searches for dust associated with the ejecta in SNRs has
turned up very little (Borkowski et al., 2006; Williams et al., 2006). As an
example, in Figure 5, we show the IR emission from DEM L71, the remnant whose
X-ray emission is shown in Figure 3.

Figure 5: Left: The Spitzer 24 $\mu$m emission from DEM L71. Middle: The
remnant in H$\alpha$. Right: The remnant in X-rays. Figure modified from
Figure 1 in Borkowski et al. (2006).

At optical wavelengths, emission from the shell of a shell-type remnant can be
seen in one of two ways. If there is at least partially neutral material
present ahead of the shock (often taken as an indication that the object is
the remnant of a Type Ia SN; see Williams et al. (2011)), then charge exchange
can take place in the post-shock environment giving rise to a strong H$\alpha$
component (Chevalier et al., 1980), as shown in the middle panel of Figure 5.
The temperatures in these remnants is too hot to allow recombination lines in
other elements. In older remnants where the temperatures have cooled
significantly, recombination lines from various elements begin to dominate
(see Section 3).

At X-ray wavelengths, shell-type SNRs may show emission from either or both of
the forward-shocked ISM/CSM and the reverse-shocked ejecta, depending on the
evolutionary state of the remnant. The physical processes causing the X-ray
emission can be either thermal, non-thermal, or both. The thermal emission is
generally a combination of both Bremmsstrahlung continuum produced when hot
ions and electrons interact, as well as line emission from highly-ionized
atoms of various metals, such as O, Fe, Ne, Si, and S. Nonthermal emission in
SNRs arises from synchrotron emission identical to that seen in radio waves,
but from much more energetic ($\sim 10-100$ TeV) electrons. Because electrons
can only be accelerated to energies this high very close to the blast wave,
synchrotron emission is most often seen from the outer edges of the shell in
shell-type SNRs Bamba et al. (2003, 2005).

Shell-type remnants with sufficient spatial extent can even be identified in
gamma-rays. For example, Katagiri et al. (2011) used Fermi observations to
show that the Cygnus Loop is spatially-resolved into a shell in the 0.5-10 GeV
energy range. The various physical processes involved in gamma-ray emission
from SNRs is beyond the scope of this chapter; for a thorough review, see
Reynolds (2008).

Figure 6: X-ray image of the Crab nebula taken with Chandra. (credit:
NASA/CXC/SAO)

###  4.2 Plerion-type SNRs

Young neutron stars born in core-collapse SNe often form pulsar wind nebulae.
Pulsar wind nebulae emit bright nonthermal emission from radio to very high
energy gamma-rays, via synchrotron or inverse Compton processes. These are
referred to as “plerion-type,” or “plerionic” SNRs Weiler & Panagia (1978).

The typical plerion-type is the Crab nebula. It is the remnant of SN 1054, and
is now one of the brightest celestial sources at most wavelengths. Figure 6
shows the X-ray image of the Crab nebula taken by Chandra. Jets extend out to
northwest and southeast from the central source (the pulsar), while two torii
are seen emanating outward from the center. The inner radius of the torus is
0.01–0.1 pc, which corresponds the shock of the relativistic winds from the
pulsar terminated by ambient material Rees & Gunn (1974); Kennel & Coroniti
(1984); Bamba et al. (2010). For more details, please see the chapter of this
work concerning PWNe.

Some plerion-type SNRs are surrounded by shells of emission resulting from the
forward shock interacting with the surrounding medium. Such remnants are often
called “composite type.” One of the typical composite type SNRs is G21.5$-$0.9
Slane et al. (2000). Figure 7 is the X-ray image of G21.5$-$0.9. One can see
both the central PWN and the outer shell.

The Crab nebula is a case of a remnant without its shell; it is still an open
issue why there is no detected shell surrounding the Crab nebula. Recent
observations by the Hitomi satellite showed that the thermal X-ray emission
around the Crab nebula is significantly fainter than other young SNRs Hitomi
Collaboration (2018). The total ejecta mass is estimated to be less than $\sim
1M_{\odot}$, implying that the progenitor explosion was relatively low energy.
One possibility is that the Crab nebula is a remnant of an electron capture
supernova Miyaji et al. (1980). More detailed observations for the Crab and
more samples are needed to understand the origin of plerion-type SNRs.

Figure 7: X-ray image of the plerion-type SNR, G21.5$-$0.9 taken with Chandra.
(credit: Credit: Heather Matheson Samar Safi-Harb (Univ. Manitoba), CXC, NASA)

###  4.3 Mixed-Morphology SNRs

Some relatively evolved SNRs show different morphology in radio and X-ray
bands; shell-like structure in radio, whereas the centrally-filed in X-rays.
They are called “Mixed-Morphology SNRs” or MM SNRs (Rho & Petre, 1998). The
X-rays are thermal, thus it is not due to its central pulsar and/or pulsar
wine nebula like plerion cases.

The mechanism by which mixed-morphology remnants are formed is still unclear.
Several hypotheses have been discussed to explain how a remnant with a radio
shell can appear centrally filled in X-rays. The first possibility is “fossil
emission” of the remnant Seward (1985); the outer shell cools down due to the
expansion and radiative cooling, and as a result, the shock temperature
becomes too low to emit thermal X-rays. Synchrotron emission from the
accelerated electrons in the shock disappear from the high energy side due to
the synchrotron cooling. Thus the shells emitting synchrotron disappear within
a few thousand years in the X-ray band, whereas they remain for more than ten
thousand years in the radio band.

Another possibility is due to “evaporated cloud”; the enhanced interior X-ray
emission could arise from gas evaporated from clouds White & Long (1991) in
the cases that shocks of SNRs propagates through a cloudy interstellar medium.
Dense and small size clouds can pass the shock and later evaporate and are
heated, and as a result makes thermal X-ray emission inside the shells.

##  5 Conclusions

Supernova remnants are an incredibly diverse class of objects. From simply a
purely morphological standpoint, their are influenced by a number of factors,
including mass loss from the progenitor system, the energy and isotropy of the
explosion itself, and the structure of the surrounding ISM/CSM. The radiation
they emit also varies from source to source: some remnants are bright in radio
waves, some in optical, some in X-rays. In general, the remnant’s state of
dynamical evolution (which can be roughly thought of as the product of age and
the density being encountered by the forward shock) is the most important
factor for determining the spectral energy distribution for a particular SNR.
Other chapters within this work will go into more detail on many aspects of
SNR study.

## References

  * Green (2019) Green, D. A. 2019, Journal of Astrophysics and Astronomy, 40, 36. doi:10.1007/s12036-019-9601-6
  * Case & Bhattacharya (1998) Case, G. L. & Bhattacharya, D. 1998, THe Astrophysical J., 504, 761. doi:10.1086/306089
  * Maggi et al. (2016) Maggi, P., Haberl, F., Kavanagh, P. J., et al. 2016, A&A, 585, A162. doi:10.1051/0004-6361/201526932
  * Sasaki et al. (2012) Sasaki, M., Pietsch, W., Haberl, F., et al. 2012, A&A, 544, A144. doi:10.1051/0004-6361/201219025
  * Garofali et al. (2017) Garofali, K., Williams, B. F., Plucinsky, P. P., et al. 2017, MNRAS, 472, 308. doi:10.1093/mnras/stx1905
  * Reynolds et al. (2008a) Reynolds, S. P., et al. 2008, ApJ, 680, 41
  * Koyama et al. (1995) Koyama, K., Petre, R., Gotthelf, E. V., et al. 1995, Nature, 378, 255
  * Aharonian et al. (2004) Aharonian, F. A., Akhperjanian, A. G., Aye, K.-M., et al. 2004, Nature, 432, 75
  * Reynolds et al. (2008b) Reynolds, S. P. 2008, ARAA, 46, 89
  * Hamilton & Sarazin (1984a) Hamilton, A. J. S. & Sarazin, C. L. 1984, The Astrophyical J., 287, 282. doi:10.1086/162687
  * Hamilton & Sarazin (1984b) Hamilton, A. J. S. & Sarazin, C. L. 1984, The Astrophysical J., 281, 682. doi:10.1086/162145
  * Hughes (1999) Hughes, J. P. 1999, The Astrophysical J., 527, 298. doi:10.1086/308082
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics, New York: Academic Press, 1959
  * Landau & Lifshitz (1959) Landau, L. D. & Lifshitz, E. M. 1959, Course of theoretical physics, Oxford: Pergamon Press, 1959
  * Chevalier (1974) Chevalier, R. A. 1974, The Astrophysical J., 188, 501. doi:10.1086/152740
  * Hughes et al. (2002) Hughes, J. P., Ghavamian, P., Rakowski, C. E., & Slane, P. O. 2003, ApJ, 582, 95
  * Rakowski et al. (2003) Rakowski, C. E., Ghavamian, P., & Hughes, J. P. 2003, ApJ, 590, 846
  * Truelove & McKee (1999) Truelove, J. K. & McKee, C. F. 1999, The Astrophysical J. Suppliment, 120, 299. doi:10.1086/313176
  * Berezhko et al. (2002) Berezhko, E. G., Petukhov, S. I., & Taneev, S. N. 2002, Astronomy Letters, 28, 632. doi:10.1134/1.1505508
  * Rees & Gunn (1974) Rees, M. J. & Gunn, J. E. 1974, MNRAS, 167, 1
  * Kennel & Coroniti (1984) Kennel, C. F. & Coroniti, F. V. 1984, ApJ, 283, 694
  * Bamba et al. (2010) Bamba, A., Mori, K., & Shibata, S. 2010, ApJ, 709, 507
  * Slane et al. (2000) Slane, P., Chen, Y., Schulz, N. S., et al. 2000, ApJL, 533, L29
  * Hitomi Collaboration (2018) Hitomi Collaboration, Aharonian, F., Akamatsu, H., et al. 2018, PASJ, 70
  * Rho & Petre (1998) Rho, J. & Petre, R. 1998, ApJL, 503, L167. doi:10.1086/311538
  * Dwek & Arendt (1992) Dwek, E., & Arendt, R. 1992, ARA&A, 30, 11
  * Borkowski et al. (2006) Borkowski, K., et al. 2006, ApJ, 642, 141
  * Williams et al. (2006) Williams, B.J., et al. 2006, ApJ, 652, 33
  * Williams et al. (2011) Williams, B.J., et al. 2011, ApJ, 741, 96
  * Chevalier et al. (1980) Chevalier, R.A., Kirshner, R.P., & Raymond, J.C. 1980, ApJ, 235, 186
  * Bamba et al. (2003) Bamba, A., Yamazaki, R., Ueno, M., et al. 2003, ApJ, 589, 827
  * Bamba et al. (2005) Bamba, A., Yamazaki, R., Yoshida, T., et al. 2005, ApJ, 621, 793
  * Katagiri et al. (2011) Katagiri, H., et al. 2011, ApJ, 741, 44
  * Reynolds (2008) Reynolds, S.P. 2008, ARA&A, 46, 89
  * Weiler & Panagia (1978) Weiler, K.W. & Panagia, N. 1978, A&A, 70, 419
  * Miyaji et al. (1980) Miyaji, S., Nomoto, K., Yokoi, K., & Sugimoto, D. 1980, PASJ, 32, 303
  * Seward (1985) Seward, F. D. 1985, Comments on Astrophysics, 11, 15
  * White & Long (1991) White, R. L. & Long, K. S. 1991, ApJ, 373, 543


# Lightweight, Uncertainty-Aware Conformalized Visual Odometry

Alex C. Stutts, Danilo Erricolo, Theja Tulabandhula, and Amit Ranjan Trivedi
Authors are with the University of Illinois Chicago (UIC), Chicago, IL, Email:
<EMAIL_ADDRESS>
###### Abstract

Data-driven visual odometry (VO) is a critical subroutine for autonomous edge
robotics, and recent progress in the field has produced highly accurate point
predictions in complex environments. However, emerging autonomous edge
robotics devices like insect-scale drones and surgical robots lack a
computationally efficient framework to estimate VO’s predictive uncertainties.
Meanwhile, as edge robotics continue to proliferate into mission-critical
application spaces, awareness of model’s the predictive uncertainties has
become crucial for risk-aware decision-making. This paper addresses this
challenge by presenting a novel, lightweight, and statistically robust
framework that leverages conformal inference (CI) to extract VO’s uncertainty
bands. Our approach represents the uncertainties using flexible, adaptable,
and adjustable prediction intervals that, on average, guarantee the inclusion
of the ground truth across all degrees of freedom (DOF) of pose estimation. We
discuss the architectures of generative deep neural networks for estimating
multivariate uncertainty bands along with point (mean) prediction. We also
present techniques to improve the uncertainty estimation accuracy, such as
leveraging Monte Carlo dropout (MC-dropout) for data augmentation. Finally, we
propose a novel training loss function that combines interval scoring and
calibration loss with traditional training metrics–mean-squared error and KL-
divergence–to improve uncertainty-aware learning. Our simulation results
demonstrate that the presented framework consistently captures true
uncertainty in pose estimations across different datasets, estimation models,
and applied noise types, indicating its wide applicability.

##  I Introduction

As machine learning (ML) becomes more prevalent in privacy, safety, and
mission-critical robotics, the ability to quantitatively and visually assess
predictive uncertainties of ML models is becoming essential for risk-aware
control and decision-making. Two types of uncertainty can result from data-
driven learning: epistemic and aleatoric. Epistemic uncertainty arises from
the variance in the training data and can often be reduced with more data. On
the other hand, aleatoric uncertainty arises from random distortions in the
data, such as blurriness, occlusions, overexposure, etc., which additional
training data cannot mitigate.

Although mitigating epistemic uncertainty is challenging, detecting,
explaining, and handling aleatoric uncertainty is even more difficult.
Therefore, for mission-critical risk-aware robotics, it is desired that
predictive models under aleatoric uncertainty must provide confidence measures
and express input/model-dependent predictive uncertainties, along with the
point (mean) prediction. Additionally, the uncertainty estimates must be
extracted with minimal additional computing cost since many edge robotic
systems have limited computing and storage capacity due to cost, footprint,
legacy design, battery power, and other factors.

Meanwhile, given their computational intensity, traditional methods for
uncertainty-aware predictions, such as Bayesian ML models, are unsuited for
edge robotics. For example, in Bayesian ML models, a posterior distribution of
weights is learned from the training data using Bayesian principles. Weights
are then sampled from the posterior, and the network output is statistically
computed and weighed against each sample’s probability. Generating
uncertainty-aware predictions, thus, involves sampling many model parameters
from the posterior and evaluating predictions at each sample, which can be
impractical under typical time and resource constraints for various edge
robotics applications.

Figure 1: Lightweight uncertainty-aware visual odometry (VO): In this paper,
we present a lightweight framework for uncertainty-aware visual odometry for
edge robotics where computing resources are limited while the predictions need
to be made in real time. Our framework exploits conformal inference and
presents four novel methodologies to extract point prediction and upper/lower
bounds under the designated uncertainty coverage.

A more computationally-efficient alternative to Bayesian inference is
Variational inference (VI) [1]. VI considers a family of approximate densities
$F$, from which a member $q^{(w)}$ is learned by minimizing the Kullback-
Leibler (KL) divergence to the posterior. However, selecting a flexible family
that closely matches the true posterior in practice can be challenging.
Additionally, VI struggles to approximate complex posteriors with multiple
modes or sharp peaks.

Recently, conformal inference (CI) of ML models, also known as conformal
prediction, was developed to overcome the above limitations of traditional
uncertainty-aware prediction frameworks[2, 3, 4, 5, 6, 7]. Unlike classical
statistical inference, which relies on data distribution to capture prediction
uncertainty and can be sensitive to model misspecification, CI is a
distribution-free approach that guarantees the statistical validity of
uncertainty intervals given a finite training sample. CI can also assess the
degree of conformity of each new observation to the available data and uses
this information to construct an uncertainty interval calibrated to a user-
specified coverage rate. Moreover, CI can be combined with any underlying
model with an inherent notion of uncertainty to provide uncertainty
quantification that is both statistically valid and model-agnostic. One of the
challenges in making CI practical is that the uncertainty estimates are often
conservative, and our work overcomes this and other limitations, providing a
compelling case for their use in robotics applications with stringent resource
usage requirements.

In particular, leveraging these advantages of CI for edge robotics, in this
work, we investigate the conformalization of visual odometry (VO), an
essential task for autonomous navigation that estimates the position and
orientation (i.e., pose) of a camera mounted on a moving vehicle relative to
the environment as the camera moves. The resulting motion estimates from VO
are used for typical autonomy objectives such as three-dimensional (3D)
reconstruction, mapping, localization, etc. Especially since VO relies solely
on cameras, which are passive, low power, and have a small footprint, VO-based
motion estimates are suitable for ego-motion tracking under stringent
footprint and battery power constraints. Fig. 1 overviews the proposed
techniques.

Exploring CI for lightweight uncertainty estimates in VO, our work makes the
following key contributions:

  * •

We present four frameworks for extracting VO-uncertainty bands by univariate
and multivariate conformalization. The presented frameworks vary in computing
workload, uncertainty accuracy, interval adaptiveness, data dependency, etc.,
aiming to provide an entire spectrum of solutions for extracting VO’s
uncertainty estimates under varying computing resources, training data, and
timing constraints for edge robotics.

  * •

We present a novel loss function for joint training of point (mean) prediction
and heteroskedastic uncertainty bounds. The loss function combines an interval
score function and combined calibration loss function with mean squared error
(MSE) loss and Kullback-Leibler (KL) divergence. Our novel loss function
improves the accuracy of uncertainty coverage and enables computationally
efficient predictive architecture for the simultaneous extraction of point
predictions and uncertainties.

  * •

We also discuss data augmentation for uncertainty-aware VO by providing
additional training set by Monte Carlo Dropout (MC-Dropout), thus enabling the
distillation of uncertainty estimates from sophisticated, computationally
expensive frameworks to conformal quantile bounds.

##  II Related Works

Under VO, the camera’s ego motion is estimated based on the changes in visual
information captured by the camera. The traditional methods utilized
techniques that involved identifying and tracking distinctive visual features,
such as corners or edges, for VO [8]. Direct methods [9] estimated the
camera’s motion by analyzing the changes in the intensity values of the pixels
in the images captured by the camera. Hybrid methods for VO combined feature-
based and direct methods; for example, feature-based techniques were used for
initial camera motion and then refined using direct methods [10]. Structure
from motion (SfM) is popularly used for simultaneously estimating the 3D
structure of the environment and camera’s motion [11].

A recent trend for VO is to use deep neural networks (DNN) to directly learn
the relation between the camera’s visual field and its ego-motion from data.
PoseNet [12] demonstrated the feasibility of training convolutional neural
networks for end-to-end tracking of the ego-motion of monocular cameras.
PoseLSTM [13] utilized long-short-term memories (LSTM) to improve ego-motion
accuracy compared to frame-based methods such as PoseNet. DeepVO [14]
demonstrated the combined strength of recurrent convolutional neural networks
(RCNNs) in learning features and modeling sequential relationships between
consecutive camera frames. UnDeepVO [15] proposed an unsupervised deep
learning approach capable of absolute scale recovery.

While many prior studies have focused on improving the accuracy of VO, only a
few have addressed the challenge of extracting its predictive uncertainties,
especially under time and computing resource constraints; for instance,
integrated deep learning-based depth predictions with particle filtering to
achieve uncertainty-aware visual localization [16]. However, particle
filtering-based uncertainty quantification can be computationally expensive
for large or high-dimensional state spaces. Kalman filters are a more
computationally efficient alternative to particle filters, but their usage is
limited to linear systems with Gaussian noise assumptions. By integrating
Monte-Carlo Dropout (MC-dropout) [17] with deep learning-based predictors,
such as PoseNet, [18] showed uncertainty-aware VO estimations. However, the
method requires performing a sufficient number of dropout iterations and
evaluating predictions for each, making it challenging to implement under time
and computing resource constraints. Besides, MC-dropout, as a basis for
uncertainty, lacks statistical rigor and ultimately leads to conservative
prediction intervals that are not particularly adaptable to the data. D3VO
[19] also utilized Bayesian techniques for uncertainty quantification,
demonstrating state-of-the-art accuracy but with a large computational
expense. With similar computational overhead, MDN-VO [20] combined an RCNN and
mixture density network for uncertainty estimations based on maximizing the
likelihood of pose estimations across a sequence of images. In contrast, the
authors in [21] presented a lightweight deterministic uncertainty estimator as
a small neural network that could be applied to VO or any data-driven deep
neural network model. The method uncovers spatial and semantic model
uncertainty with significantly less computation but lacks statistical
guarantees.

Unlike the above, our methods build on conformal inference, which can provide
statistically valid uncertainty intervals while not requiring heavy
computational budgets. Combining CI with an underlying model with an inherent
notion of uncertainty (e.g., conditional quantile regression) guarantees true
value coverage within the uncertainty intervals even when the model’s
predictions are poor. However, CI alone is not highly adaptable to
heteroskedastic data, as prediction sets can appear fixed and weakly dependent
on the model’s predictors. For example, conformal classification can produce
rather conservative prediction sets given their sole reliance on softmax
scores, as mentioned in [5]. For regression problems, conformalized quantile
regression (CQR) was introduced in [22]. CQR produces adaptive and
distribution-free prediction intervals with a guaranteed marginal coverage
rate, such as 90%. While typical full conformal prediction assumes that
samples are drawn exchangeably (i.e., conditionally, i.i.d, on the joint
probability distribution function between input and labels), CQR departs from
this assumption to achieve a finite sample coverage guarantee through a
technique called split conformal prediction. Various forms of CI, including
conformal classification and CQR, were initially developed for one-dimensional
(1D) data. They have since been extended to multivariate cases, making them
particularly valuable for our application to VO.

Table I: Qualitative comparison of proposed techniques on VO-uncertainty band
extraction
---
Method |  CQR (Sec. III-A) |  CSP (Sec. III-B) |  MCQR w/ MC-dropout (Sec. IV-A) |  CJP (Sec. IV-B)
Uncertainty Bands |  Univariate |  Univariate |  Multivariate |  Multivariate
Disjoint Bands |  No |  Yes |  No |  No
Computational Complexity |  Low |  Low |  High |  Medium
Interval Calibration |  Fixed |  Fixed |  Fixed |  Tunable
Interval Adaptiveness |  Conservative |  Semi-tunable |  Flexible |  Tunable
Uncertainty Accuracy |  Low |  Medium |  High |  High
Training Time |  Short |  Medium |  Long |  Medium
Data Dependency |  Low |  High |  High |  Low

##  III Extracting Univariate Bands of Uncertainties for Visual Odometry (VO)

This section introduces two methods for extracting predictive uncertainty
bands in VO: univariate conformalized quantile regression and conformalized
set prediction. The methods presented here are lightweight but produce
rectangular uncertainty bands, which can be too conservative. The first method
uses quantile regression on each coordinate of the position and orientation
output vectors of VO, which are combined to produce an overall quantile
region. The second method maps the regression problem in VO to a
classification problem and performs a set prediction, which can produce
disjoint uncertainty bands. The presented techniques are compared in Table I
and discussed below:

###  III-A VO-Uncertainty by Univariate CQR

We employ conformalized quantile regression (CQR) [22] to extract univariate
uncertainty bands for each output coordinate in VO. To test our approach, we
use PoseNet [12] with a ResNet34-based feature extractor for conformalization.
However, the presented techniques are generalizable to other feature-based and
direct predictors for VO. The PoseNet model outputs translational elements
$x$, $y$, and $z$, as well as rotational elements $p$, $q$, and $r$ of angle
$w$ in radians. To reduce the number of response variables from seven to six,
we convert the orientation quaternion ($w,p,q,r$) to Euler angles (roll
$\phi$, pitch $\theta$, yaw $\psi$). After pre-training, we apply CQR to
obtain prediction intervals along each one-dimensional variable, which reveals
the predictive uncertainty. Fig. 2 depicts the process for this method.

Figure 2: VO uncertainty estimation using univariate CQR: Features of the
input image and extracted and regressed to six-dimensional pose. Subsequently,
each pose dimension’s predictions are conformalized using the held-out
calibration data. Finally, multivariate uncertainty regions are predicted by
multiplying 1D uncertainty bands.

The main goal of CQR is to produce prediction intervals where the probability
of true realizations falling within these intervals is guaranteed to be
equivalent to or greater than some chosen coverage rate $p$, i.e.,
$\mathbb{P}\\{Y\in C(X)\mid X=x\\}\geq p$. This property is known as marginal
coverage, and is calculated over all training and test set samples. To compute
$C(X)$, we use Random Forest-based quantile regression, which is less prone to
over-fitting, requires fewer data and computing, and results in lower and
upper bound percentiles that capture heteroskedasticity by expanding and
contracting based on the underlying uncertainty in the prediction. The
quantiles are then conformalized (corrected) by a calibration set that is
split a priori from the training set. The following equations represent the
three steps of the procedure:

$\begin{split}\\{{Q_{l},Q_{h}}\\}&\leftarrow QR(\\{(X_{i},Y_{i})\\}),\ i\in
I_{train}\end{split}$ (1)
$\begin{split}E_{i}&\leftarrow\max\\{{Q_{l}(X_{i})-Y_{i},Y_{i}-Q_{h}(X_{i})}\\},\
i\in I_{cal}\end{split}$ (2)
$\begin{split}C_{cal}(X)&=[Q_{l}(X)-Q_{cal}(E,I_{cal}),\\\
&Q_{h}(X)+Q_{cal}(E,I_{cal})],\textrm{ where}\end{split}$ (3a)
$\begin{split}Q_{cal}&(E,I_{cal}):=(1-\alpha)(1+1/|I_{cal}|)-\textrm{th}\\\
&\textrm{quantile of }\\{E_{i}:i\in I_{cal}\\}\end{split}$ (3b)

Here, given a certain miscoverage rate $\alpha$, the quantile regression
algorithm ($QR$) fits two conditional quantile functions $Q_{l}$ and $Q_{h}$
on the training set $I_{train}$, which contains samples $X$ and labels $Y$.
The effectiveness of the initial prediction interval $[Q_{l}(X),Q_{h}(X)]$ in
covering $Y_{i}$ is then measured using conformal scores $E_{i}$, which are
evaluated on the held out calibration set $I_{cal}$. If $Y_{i}$ is outside the
boundaries, $E_{i}$ measures the distance from the nearest boundary. If
$Y_{i}$ falls within the desired boundaries, $E_{i}$ measures the larger of
the two distances, accounting for both undercoverage and overcoverage.
Finally, the calibrated prediction interval $C_{cal}$ is formed using
$Q_{cal}(E,I_{cal})$, which is the empirical $(1-\alpha)(1+1/|I_{cal}|)$-th
quantile of ${{E_{i},\ i\in I_{cal}}}$. Using a calibration set ensures that
the resulting intervals have the desired miscoverage rate.

Figure 3: Uncertainty-aware VO by conformal set prediction (CSP): (a) The
architecture of CSP for VO. The output of the feature extraction unit is fed
to parallel 1D softmax classifiers. Subsequently, classifier outputs are
conformalized for the predictive set generation along each dimension. Next,
the sets are multiplied and projected to 3D space to generate uncertainty-
aware rectangular predictive regions. (b) Our approach also leverages density-
aware classification boundaries along each dimension. Based on the training
set, frequently visited locations are packed into classes with a narrow gap
between upper and lower boundaries for high precision determination when
projecting to 3D space.

###  III-B VO-Uncertainty by Conformalized Set Prediction (CSP)

This section describes our method, conformalized set prediction (CSP), for
converting each of the multiple univariate pose regression problems in VO to
pose classification using one-hot encoding and generating conformalized
prediction sets by discretizing the flying space. We use the same feature
extractor as discussed in the previous section. However, instead of relying on
quantile regression to determine the model’s uncertainties, we leverage the
model’s normalized exponential scores (i.e., softmax scores) for conformal set
prediction. A notable advantage of this approach is the ability to predict
disjoint regions of uncertainty, unlike the previous one that only predicts
contiguous regions of uncertainty.

The approach begins by discretizing the drone’s flying space into $K$ sets
along each dimension $x$, $y$, and $z$. To achieve this, a non-uniform space
discretization is followed, wherein the histograms of the training set
trajectories along each dimension are divided into $K$ quantiles to determine
class boundaries. The result is a one-hot encoding matrix based on $K$, and a
lightweight neural network classifier is built with a programmable output
layer of size $K$. For each dimension of the pose, a network is trained to
output reliable softmax scores that can be conformalized. Fig. 3 illustrates
the softmax classifiers for each dimension sharing their feature extraction
head for parameter sharing and computational efficiency. The product of the
predicted set of classes along each dimension produces the net uncertainty
regions. Notably, this method can represent frequently visited locations with
higher precision, represented by a class with a narrower gap between upper and
lower boundaries. Conversely, less frequently visited locations are
represented with lower precision.

Under the above classification treatment of VO, to perform conformal set
prediction, a calibration set is first separated from the training set to
compute conformal scores. It is essential to ensure the statistical guarantee
of conformal prediction by ensuring that the average probability of correct
classifications within the prediction sets $C(X)\subset{1,...,K}$ is almost
exactly $1-\alpha$, where $\alpha$ denotes an arbitrary miscoverage rate
(e.g., 10%). This marginal coverage property can be expressed mathematically
as $1-\alpha\leq\mathbb{P}{Y\in C(X)\mid X=x}\leq 1-\alpha+\frac{1}{n+1}$,
where $n$ is the size of the calibration set. Conformal scores are then
obtained by subtracting the softmax of the correct class for each input from
one. Finally, $\hat{q}$ is computed as the $ceil(\frac{(n+1)(1-\alpha)}{n})$
empirical quantile of these conformal scores, and the conformalized prediction
set $C(X)\mid X=x$ is formed by incorporating classes with softmax scores
greater than $(1-\hat{q})$.

Importantly, since the predicted set along each dimension may include classes
that are not proximal, the above method has the ability to produce disjoint
uncertainty regions. In contrast, the previous method in Sec. III-A (and
multivariate methods presented in the later discussion) can only output
contiguous uncertainty regions. The class labels in the above conformal set
prediction can also be determined by discretizing the entire 3D space;
however, the necessary classes for matching precision to dimension-wise
discretization would grow exponentially. We also found that the above
conformalization procedure may sometimes be too strict, especially when the
conditional softmax outputs are not representative. However, it can be
improved by using the softmax outputs of all classes in gathering the
conformal scores as in [5].

##  IV Extracting Multivariate Bands of Predictive Uncertainties in Visual
Odometry (VO)

This section presents two more approaches for extracting predictive
uncertainty bands in VO: multivariate conformalized quantile regression with
MC-dropout, and conformalized joint prediction. The first approach
demonstrates the novel usage of MC-dropout as a data augmentation technique in
improving the performance of a conditional variational autoencoder (CVAE). The
second approach jointly trains predictive uncertainty bands and pose
estimation using a new loss function, thus reducing computational resources
for applications where uncertainty estimates are always required.

###  IV-A VO-Uncertainty by Multivariate CQR (MCQR)

To generate more informative VO-uncertainty bands that consider the
correlation among pose dimensions, we adopt the multivariate conformalized
quantile regression algorithm in [23] and enhance it with MC-dropout (MCQR w/
MC-dropout). This method entails utilizing a CVAE to acquire a proper latent
representation of the sample distribution $Y\mid X$ and applying an extension
of directional quantile regression (DQR) to create quantile regions. The
architecture of the CVAE, shown in Fig. 4(a), is a variant of a variational
autoencoder (VAE) that conditions the generative model on labels during both
the encoding and decoding steps. This additional information improves the
model’s ability to produce targeted outputs. Instead of directly fitting
quantile functions to the data with simple regression, a separate neural
network is developed to learn the best conditions for sampling the latent
representation $Z\mid X$ of the complete multivariate response variable based
on $Y\mid X$. At the same time, a quantile region is generated in this
dimension-reduced latent space, which has a more informative, flexible, and
arbitrary shape when propagated through the decoder.

Compared to the procedure for extracting univariate bands in Sec. III-A,
extracting multivariate bands of uncertainty here using a VAE predictor
requires more training data due to the model’s complexity and higher number of
parameters. Therefore, to improve the accuracy of VAE-based uncertainty bands,
we utilize MC-dropout to generate additional training data. Training VAEs with
the additional data enables uncertainty distillation from the MC-dropout
procedure into lightweight conformalized quantile bands. However, the
underlying pose estimation model’s performance heavily influences this
approach’s effectiveness. Poor predictions of mean estimates can hinder
learning with additional data.

###  IV-B VO-Uncertainty by Conformalized Joint Prediction (CJP)

This section discusses the joint training of multivariate VO-uncertainty bands
and the mean pose, i.e., conformalized joint prediction (CJP). The approach
reduces the computational costs of uncertainty-aware VO by sharing parameters
and network layers for mean prediction and uncertainty bands, which is
particularly useful for platforms where uncertainty estimation must always be
ON.

Fig. 4(b) shows the proposed VAE architecture for simultaneously extracting
pose predictions and upper/lower uncertainty bounds based on true conditional
quantiles. In the network, we adopt a parametric rectified linear unit (PReLU)
as an activation function to introduce additional learnable parameters, better
handle negative values and avoid the dying ReLU problem (values become zero
for any input). In Fig. 4, the network first extracts features from the input
image using a pre-trained convolutional neural network (e.g., MobileNet [24]),
followed by encoding to a 3D latent space with a spherical Gaussian
probability distribution. Next, a sample is taken from the latent space and
propagated through the decoder to generate the multivariate pose and
uncertainty-bound estimations at the desired coverage rate.

Figure 4: VAE architecture for joint training of point (mean) prediction and
uncertainty bounds: (a) Predictions from a pose estimation network (such as
PoseNet) are fed to a conditional VAE for pose reconstruction and uncertainty
estimation by first encoding to the latent space and then generatively passing
through a decoder. (b) Comparatively, the network here jointly learns point
prediction and uncertainty estimates directly from the features extracted from
camera images. The novel loss function in Eq. (6) is used for the integrated
training. The simultaneous extraction of mean and uncertainty bounds provides
better computational efficiency for applications where uncertainty awareness
must always be ON.

The pinball loss function, a.k.a. quantile loss, used for conformalized
quantile regression in [22, 23] is inadequate for such integrated training of
the mean prediction and uncertainty bands. Thus, we propose a new training
loss function combining the traditional mean-squared error and Kullback-
Leibler (KL) divergence losses of a VAE, $\text{MSE}_{loss}$ and
$\text{KL}_{loss}$, with two additional loss components
$\text{INTSCORE}_{loss}$ and $\text{COMCAL}_{loss}$. This total loss function,
$\mathcal{L}_{Total}$, capitalizes on two unique findings reported in [25].

The first finding in the proposed loss function is the unconventional use of
the interval score function (i.e., Winkler score [26, 27]) to train an
uncertainty quantification model to consistently output centered prediction
intervals within a specified quantile range. Notably, the interval score
function is typically used to evaluate prediction intervals, not optimize
them. This approach is advantageous for generating informative prediction
intervals over pose trajectories. We choose two percentiles (e.g.,
$\alpha_{l}=$ 5% and $\alpha_{h}=$ 95%) as uncertainty bounds, compute the
interval score loss shown below for each percentile, and then take their
expected value:

$\begin{split}\text{INTSCORE}_{loss}&=(Q_{h}-Q_{l})+\frac{2}{\alpha}(Q_{l}-y)\mathbb{I}\\{y<Q_{l}\\}\\\
&+\frac{2}{\alpha}(y-Q_{h})\mathbb{I}\\{y>Q_{h}\\}\\\ \end{split}$ (4)

Here, $y$ represents the pose label, $Q_{l}$ and $Q_{h}$ are each dimension’s
lower and upper quantile estimates, and $\mathbb{I}$ is the indicator
function. Including the interval score in $\mathcal{L}_{Total}$ is
instrumental in linking pose reconstruction optimization via
$\text{MSE}_{loss}$ to quantile optimization over each dimension while
honoring the dimensions’ covariance and correlation.

The second finding is combined calibration loss, denoted as
$\text{COMCAL}_{loss}$. $\text{COMCAL}_{loss}$ establishes an explicit and
controllable balance between prediction interval calibration and sharpness
rather than relying on the loosely implicit balance in basic pinball loss.
$\text{COMCAL}_{loss}$ comprises two objectives: one that minimizes the
difference between $p^{cov}_{avg}$ and the chosen marginal coverage rate $p$
($\text{CAL}_{obj}$), and another that minimizes the distance between $Q_{l}$
and $Q_{h}$ ($\text{SHARP}_{obj}$). $p^{cov}_{avg}$ is the estimated average
probability that the pose label values lie within $[\mathit{Q_{l},Q_{h}}]$.
The two objectives may conflict with each other since one aims to minimize
over-coverage and under-coverage, while the other aims to minimize the length
of the prediction interval. This necessitates including a hyper-parameter
$\lambda$ to strike the appropriate balance. Thus, $\text{COMCAL}_{loss}$ is
defined as follows:

$\begin{split}&\text{CAL}_{obj}=\\\
&\mathbb{I}\\{p^{cov}_{avg}<p\\}\times\frac{1}{N}\sum_{i=1}^{N}[(y_{i}-Q_{l,h}(x_{i}))\mathbb{I}\\{y_{i}>Q_{l,h}\\}]+\\\
&\mathbb{I}\\{p^{cov}_{avg}>p\\}\times\frac{1}{N}\sum_{i=1}^{N}[(Q_{l,h}(x_{i})-y_{i})\mathbb{I}\\{y_{i}<Q_{l,h}\\}]\end{split}$
(5a)
$\begin{split}\text{SHARP}_{obj}&=\frac{1}{N}\sum_{i=1}^{N}\begin{cases}Q_{l}(x_{i})-Q_{h}(x_{i}),&\text{$p\leq
0.5$}\\\ Q_{h}(x_{i})-Q_{l}(x_{i}),&\text{$p>0.5$}\end{cases}\end{split}$ (5b)
$\text{COMCAL}_{loss}=(1-\lambda)\times\text{CAL}_{obj}+\lambda\times\text{SHARP}_{obj}$
(5c)

where $x_{i}$ and $y_{i}$ are individual images and labels of a training batch
and $Q_{l,h}$ is meant to denote that $\text{CAL}_{obj}$ is computed
separately for both $Q_{l}$ and $Q_{h}$ and then summed.

The complete loss function combining MSE loss, KL divergence, interval score,
and calibration loss in our approach is thus given as:

$\begin{split}\mathcal{L}_{Total}&=\text{MSE}_{loss}(y,\hat{y})+\text{KL}_{loss}(\mu,log(\sigma^{2}))\\\
&+\text{INTSCORE}_{loss}(y,\mathit{Q_{l},Q_{h}},\\{\mathit{\alpha_{l},\alpha_{h}}\\})\\\
&+\text{COMCAL}_{loss}(y,p^{cov}_{avg},\mathit{Q_{l},Q_{h}})\\\ \end{split}$
(6)

where $\hat{y}$ is the reconstructed pose. $\mu$ and $log(\sigma^{2})$ are the
hyper-parameters of the VAE’s latent space. We estimate $p^{cov}_{avg}$ for
each pose dimension with every training batch to enable flexible calibration
during training.

Traditional CQR involves splitting the training set into two subsets and
conformalizing the quantiles afterward (i.e., split conformal prediction).
This is done to avoid full conformal prediction, which requires many more
calibration steps than just one, sacrificing speed for statistical efficiency.
However, cross-conformal prediction provides a viable alternative for
conformalization by striking a reasonable balance between computational and
statistical efficiency, utilizing multiple calibration steps across the entire
training set [28]. Our approach is similar, representing a unique form of
cross-conformal prediction. We assume that $p^{cov}_{avg}$, which is averaged
over randomized training batches, is sufficiently close to an average over the
entire training set.

As will be demonstrated later in Sec. V, this approach yields superior
outcomes compared to the other methods with a small increase in model
complexity and size. The results underscore the method’s accuracy,
adaptability, and flexibility in producing conformalized prediction intervals
for multivariate pose data. The method maintains its strength even when
swapping the underlying model from a larger, more accurate one, such as
ResNet34, to one better suited for lightweight edge applications like
MobileNetV2. Furthermore, its advantages are consistent across multiple
datasets, implying its robustness to numerous applications.

Figure 5: Comparison of uncertainty-aware prediction intervals over 1D and 2D
pose estimations from the four proposed conformalized visual odometry (CVO)
methods in Secs. III and IV: Differences in average interval length,
contiguity, precision, accuracy, coverage, and adaptiveness are highlighted
across the methods using Scene 2 of RGB-D Scenes Dataset v.2.

##  V Results and Discussions

This section presents a chronological review of the results from the four
frameworks proposed for conformalized visual odometry (CVO) and compared in
Table I. The aim of these methods is to develop reliable prediction intervals
that can capture true uncertainty in data-driven pose estimation with
guaranteed confidence. While each approach is based on PoseNet [12], they can
be applied to other pose estimation models. In highlighting the differences
between the methods, we use the data from scene 2 of RGB-D Scenes Dataset v.2
[29] as visuals of its pose trajectory are easy to comprehend. All results
illustrate uncertainty bands that exhibit at least 90% guaranteed marginal
coverage through CI.

In Sec. III-A, we introduced univariate conformalized quantile regression
(CQR) to generate uncertainty bands. However, as shown in Fig. 5, the
resulting prediction interval is overly cautious and box-shaped despite its
ability to adjust to the data’s heteroscedasticity. This is because the method
employs quantile regression on the distinct variables of the multivariate pose
response and then merges them without considering their covariance and
correlation. As a result, this method lacks sufficient information and is not
particularly useful in most cases, despite its low computational complexity.
This drawback is consistent across all datasets.

In Sec. III-B, we presented conformalized set prediction (CSP) to generate
uncertainty bands. As illustrated in Fig. 5, this approach adapts well to the
data, even with non-continuous uncertainty boundaries. The method’s
adaptiveness is deeply dependent on the model’s softmax scores and the proper
discretization of pose values, which can be adjusted by changing the number of
classes $K$. Like the first method, this approach has difficulty handling
multivariate data, which limits its ability to estimate true uncertainty
accurately. Although the average interval length is significantly shorter, the
conservativeness of the rectangular uncertainty region has only scaled
linearly, limiting its usefulness.

Figure 6: Comparison of uncertainty-aware prediction intervals over 3D pose
estimations from CJP under different noise conditions: After switching from
the ResNet34 feature extractor to MobileNetV2, gaussian, salt and pepper, and
speckle noise were applied to Scene 2 of RGB-D Scenes Dataset v.2 to further
stress test the method and demonstrate its resilience. These conditions are
meant to simulate real-world disturbance.

In Sec. IV-A, we introduced the third method, which utilizes multivariate
conformalized quantile regression with MC-dropout (MCQR w/ MC-dropout) to
generate uncertainty bands. This approach showcases the novel use of MC-
dropout as a data augmentation technique, in addition to regularization and
uncertainty estimation, to enhance learning. This method is the first step
towards addressing the limitations of the previous methods by considering the
relationships between pose dimensions. As shown in Fig. 5, this method has the
tightest uncertainty estimates, but they represent the opposite extreme of the
conservatism observed in the univariate methods. However, this is inconsistent
across all datasets; the method performs much better in some other datasets,
indicating that it is highly data-dependent. This inconsistency in uncertainty
accuracy is particularly noticeable in smaller datasets, as is the case here.
Furthermore, this method has the highest computational complexity, requiring
training and optimizing multiple networks. Although it can potentially produce
flexible and highly accurate uncertainty intervals, its inconsistency across
datasets and computational inefficiency limit its wide applicability.

The final method, described in Sec. IV-B demonstrates a novel combination of
multivariate conformalized quantile regression and the mean prediction to
produce uncertainty bounds using a new loss function that remarkably achieves
multiple objectives simultaneously. We deem this approach conformalized joint
prediction (CJP). Fig. 5 demonstrates its significant improvements over the
other methods. Firstly, the prediction interval optimally expands and
contracts in response to the data, being neither too dispersed nor sharp in
each plot. The increase in interval length (and hence uncertainty) at the
middle and end points in the trajectory corresponds to worse lighting
conditions in the scene’s image sequence. Secondly, with the tunable
calibration parameter $\lambda$ $\in$ $\\{0,1\\}$, we can optimize the network
for sharpness and/or coverage. Here, $\lambda$ was set to 0.5 to showcase the
balanced standard. Additionally, Fig. 6 shows results of this method when
switching the feature extraction model from ResNet34 to MobileNetV2,
effectively reducing its learning capacity by over 80%. It handles the same
image sequence with comparable accuracy even when introducing different forms
of noise, such as gaussian, salt and pepper, and speckle. Notably, the
uncertainty increases with this transition and even more with added noise. The
speckle noise exasperates light variation the most and thus has the largest
average interval length. An increase in the average uncertainty intervals
indicates that the framework can capture the sensor’s uncertainties along with
the model’s predictive uncertainties. The optimality of this method in
extracting true pose estimation uncertainty is consistent across several
datasets and conditions, indicating that it can fit many applications.

##  VI Conclusion

This paper introduced and compared four novel frameworks for capturing
aleatoric uncertainty in data-driven VO through conformal inference while
considering computational resource limitations. The presented framework trade-
off essential metrics of uncertainty-aware inference, such as computational
workload, adaptiveness of uncertainty interval, training time, etc., aiming to
provide a spectrum of solutions against widely varying robustness and resource
constraints on edge robotics. We discussed the architectures of generative
deep neural networks for estimating multivariate uncertainty bands along with
point (mean) prediction, data augmentation techniques to improve the accuracy
of uncertainty estimation, and novel training loss functions that combined
interval scoring and calibration loss with traditional training metrics such
as KL-divergence to improve uncertainty-aware learning.

## References

  * [1] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, “Variational Inference: A Review for Statisticians,” _Journal of the American Statistical Association_ , vol. 112, no. 518, pp. 859–877, 2017.
  * [2] V. Vovk, A. Gammerman, and G. Shafer, _Algorithmic Learning in a Random World_. Berlin, Heidelberg: Springer-Verlag, 2005.
  * [3] J. Lei, M. G’Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman, “Distribution-Free Predictive Inference for Regression,” 2018.
  * [4] G. Shafer and V. Vovk, “A Tutorial on Conformal Prediction,” _J. Mach. Learn. Res._ , vol. 9, p. 371–421, jun 2008.
  * [5] A. N. Angelopoulos and S. Bates, “A Gentle Introduction to Conformal Prediction and Distribution-Free Uncertainty Quantification,” 2021. [Online]. Available: https://arxiv.org/abs/2107.07511
  * [6] B.-S. Einbinder, Y. Romano, M. Sesia, and Y. Zhou, “Training uncertainty-aware classifiers with conformalized deep learning,” in _NeurIPS_ , 2022.
  * [7] S. Messoudi, S. Rousseau, and S. Destercke, “Deep conformal prediction for robust models,” in _Information Processing and Management of Uncertainty in Knowledge-Based Systems_. Cham: Springer International Publishing, 2020, pp. 528–540.
  * [8] D. Scaramuzza and F. Fraundorfer, “Visual odometry [tutorial],” _IEEE robotics & automation magazine_, vol. 18, no. 4, pp. 80–92, 2011.
  * [9] H. Alismail, M. Kaess, B. Browning, and S. Lucey, “Direct visual odometry in low light using binary descriptors,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 444–451, 2016.
  * [10] N. Krombach, D. Droeschel, and S. Behnke, “Combining feature-based and direct methods for semi-dense real-time stereo visual odometry,” in _Intelligent Autonomous Systems 14: Proceedings of the 14th International Conference IAS-14 14_. Springer, 2017.
  * [11] J. Gräter, T. Schwarze, and M. Lauer, “Robust scale estimation for monocular visual odometry using structure from motion and vanishing points,” in _2015 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2015, pp. 475–480.
  * [12] A. Kendall, M. Grimes, and R. Cipolla, “PoseNet: A Convolutional Network for Real-Time 6-DOF Camera Relocalization,” in _2015 IEEE International Conference on Computer Vision (ICCV)_ , 2015.
  * [13] F. Walch, C. Hazirbas, L. Leal-Taixé, T. Sattler, S. Hilsenbeck, and D. Cremers, “Image-Based Localization Using LSTMs for Structured Feature Correlation,” in _2017 IEEE International Conference on Computer Vision (ICCV)_ , 2017, pp. 627–637.
  * [14] S. Wang, R. Clark, H. Wen, and N. Trigoni, “DeepVO: Towards end-to-end visual odometry with deep Recurrent Convolutional Neural Networks,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , 2017, pp. 2043–2050.
  * [15] R. Li, S. Wang, Z. Long, and D. Gu, “UnDeepVO: Monocular Visual Odometry Through Unsupervised Deep Learning,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018.
  * [16] P. Shukla, A. C. Stutts, S. Ravi, T. Tulabandhula, A. R. Trivedi _et al._ , “Robust Monocular Localization of Drones by Adapting Domain Maps to Depth Prediction Inaccuracies,” _arXiv preprint arXiv:2210.15559_.
  * [17] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _International Cnference on Machine Learning_. PMLR, 2016, pp. 1050–1059.
  * [18] G. Costante and M. Mancini, “Uncertainty Estimation for Data-Driven Visual Odometry,” _IEEE Transactions on Robotics_ , 2020.
  * [19] N. Yang, L. von Stumberg, R. Wang, and D. Cremers, “D3VO: Deep Depth, Deep Pose and Deep Uncertainty for Monocular Visual Odometry,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 1278–1289.
  * [20] N. Kaygusuz, O. Mendez, and R. Bowden, “MDN-VO: Estimating Visual Odometry with Confidence,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2021.
  * [21] M. Lee, B. Mudassar, and S. Mukhopadhyay, “Lightweight Model Uncertainty Estimation for Deep Neural Object Detection,” in _2022 International Joint Conference on Neural Networks (IJCNN)_ , 2022.
  * [22] Y. Romano, E. Patterson, and E. J. Candès, _Conformalized Quantile Regression_. Red Hook, NY, USA: Curran Associates Inc., 2019.
  * [23] S. Feldman, S. Bates, and Y. Romano, “Calibrated multiple-output quantile regression with representation learning,” _Journal of Machine Learning Research_ , vol. 24, no. 24, pp. 1–48, 2023.
  * [24] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam, “MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications,” 2017. [Online]. Available: https://arxiv.org/abs/1704.04861
  * [25] Y. Chung, W. Neiswanger, I. Char, and J. Schneider, “Beyond pinball loss: Quantile methods for calibrated uncertainty quantification,” _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * [26] T. Gneiting and A. Raftery, “Strictly Proper Scoring Rules, Prediction, and Estimation,” _Journal of the American Statistical Association_ , vol. 102, pp. 359–378, 03 2007.
  * [27] R. L. Winkler, “A Decision-Theoretic Approach to Interval Estimation,” _Journal of the American Statistical Association_ , 1972.
  * [28] V. Vovk, “Cross-conformal predictors,” _Annals of Mathematics and Artificial Intelligence_ , vol. 74, 08 2012.
  * [29] K. Lai, L. Bo, and D. Fox, “Unsupervised feature learning for 3D scene labeling,” in _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2014, pp. 3050–3057.