# Direct high-order edge-preserving regularization for tomographic image
reconstruction

Daniil Kazantsev1,2, Evgueni Ovtchinnikov3, William R. B. Lionheart5, Philip
J. Withers1,2, Peter D. Lee1,2 1 The Manchester X-Ray Imaging Facility, School
of Materials, The University of Manchester, Manchester, M13 9PL, UK  2 The
Manchester X-Ray Imaging Facility, Research Complex at Harwell, Didcot,
Oxfordshire, OX11 0FA, UK  3 Visualization Group, Research Complex at Harwell,
STFC, Oxfordshire, OX11 0QX, UK  4 School of Mathematics, Alan Turing
Building, The University of Manchester, M13 9PL, UK
<EMAIL_ADDRESS>
###### Abstract

In this paper we present a new two-level iterative algorithm for tomographic
image reconstruction. The algorithm uses a regularization technique, which we
call edge-preserving Laplacian, that preserves sharp edges between objects
while damping spurious oscillations in the areas where the reconstructed image
is smooth. Our numerical simulations demonstrate that the proposed method
outperforms total variation (TV) regularization and it is competitive with the
combined TV-$\ell_{2}$ penalty. Obtained reconstructed images show increased
signal-to-noise ratio and visually appealing structural features. Computer
implementation and parameter control of the proposed technique is
straightforward, which increases the feasibility of it across many tomographic
applications. In this paper, we applied our method to the under-sampled
computed tomography (CT) projection data and also considered a case of
reconstruction in emission tomography The MATLAB code is provided to support
obtained results.

##  1 Introduction

Frequently, in X-ray computed tomography (CT), the amount of collected
projection data is lower that it is required by the Nyquist sampling theorem
[1]. In medical imaging, the restrictions are applied to minimize the ionizing
radiation which can harm living tissue cells [2]. In material science, the aim
is to better resolve temporal resolution via higher frame rate acquisition
[3]. In such cases of limited data, iterative techniques can provide better
reconstructions than analytical methods [4].

Dealing with ill-posed and ill-conditioned inverse problems, iterative methods
require regularization to constrain the space of desirable solutions. Due to
edge-preserving properties, total variation (TV) regularization [5] has been
extensively used in tomographic iterative reconstruction for the last three
decades. However, TV penalty produces piecewise-constant images (so-called
“cartoon” or “staircasing” effect) even if the original object is smooth [6].
This “cartoon” effect can be particularly undesirable in emission tomography
(ET) [2], such as Positron Emission Tomography (PET) or Single-Photon Emission
Computed Tomography (SPECT) where due to limited resolution of the imaging
system and low photon-number statistics, reconstructed images are naturally
blurred. Since in ET the activity distribution is piecewise-smooth and the use
of TV penalty can result in undesirable artifacts, this can potentially bias
the following clinical interpretation of reconstructed images [7].

Various improvements have been made to reduce “cartoon” appearance of the
recovered images [8, 9, 10, 11, 12, 13, 14]. The most straightforward way is
to couple TV term with $\ell_{2}$ norm [8, 9, 12, 13]. The goal of this
approach is to establish a trade-off (by means of regularization constants)
between piecewise-constant and piecewise-smooth regularization terms. This,
however, complicates the optimization problem and for practical reconstruction
purposes it can be a challenging task to establish the required regularization
parameters. Therefore, the use of a single penalty term which can accommodate
the required properties is highly desirable for tomographic reconstruction.
Some single term penalties have been proposed for image denoising and they
based on the edge preserving Laplacian [11, 14] or generalized forms of TV
norm [10].

In this paper, we propose a high-order penalty which has similarities with the
Laplacian based methods [11, 14] but it also possess some of unique qualities.
We adapt our high-order penalty to the tomographic reconstruction problem
using the two-step fixed point iteration algorithm [15]. To demonstrate
applicability of the proposed regularizer to image reconstruction problems we
compare it with the classical TV [5] and combined TV-$\ell_{2}$ methods for CT
and ET image reconstruction. All techniques are compared quantitatively and
visually.

##  2 Method

###  2.1 Image reconstruction problem

Tomographic image reconstruction problem consists of determining the inner
structure of an object based on its X-ray observations from the several
different angular positions. Incoming photons with different energies (due to
absorption) are registered by detectors and information about the path length
a photon has travelled along the line can be decoded [1]. Therefore, by
solving the inverse problem where projection data is given, the degree of
absorption or attenuation coefficient of the object can be recovered. In
mathematical terms, the reconstruction problem can be formulated as the least
squares (LS) problem:

$\displaystyle\hat{\bm{u}}=\operatorname*{argmin}\limits_{\bm{u}}\left\|\mathrm{A}\bm{u}-\bm{b}\right\|_{2}^{2},$
(1)

where $\bm{b}\in\mathbb{R}^{M}$ is a discrete function of the number of the
detector bins and the observation angles describing the projection data
(sinogram), $\bm{u}\in\mathbb{R}^{N}$ is a function of spacial variables
describing the observed object (e.g. density of the object’s material) and
$\mathrm{A}:\mathbb{R}^{{N}}\rightarrow\mathbb{R}^{{M}}$ is a sparse system
projection matrix (discrete approximation of the continuous Radon transform
[1]) mapping the “space of objects” to the “space of observations”. In
continuous space the operator $\mathrm{A}$ is an integral operator, and zero
is a condensation point of its singular values, which makes the problem (1)
ill-posed [16].

Depending on the numbers of spacial grid points, detectors and angles, the
matrix $\mathrm{A}$ can be “fat” (the number of rows is less than the number
of columns) or “tall” (the number of rows is greater than the number of
columns). Irrespective of the shape, the singular values of $\mathrm{A}$ form
a very tight cluster near zero owing to the same property of the integral
operator from which $\mathrm{A}$ derives [16].

###  2.2 Regularization

The quadratic functional
$J(\bm{u})=\left\|\mathrm{A}\bm{u}-\bm{b}\right\|_{2}^{2}$ with its “normal”
form $\mathrm{A}^{*}\mathrm{A}\bm{u}=\mathrm{A}^{*}\bm{b}$ can be minimized by
a suitable iterative minimization algorithm, e.g. Conjugate Gradient Least
Squares (CGLS) algorithm [17], however, the convergence of iterations can be
very slow because of the poor conditioning of the Hessian
$\mathrm{A}^{*}\mathrm{A}$ [15]. The slow convergence of iterations is closely
related to another difficulty in dealing with this kind of problem, which is
the ill-posedness of the problem. Generally, the solution $\hat{\bm{u}}$ of
the problem (1) is not unique (if $\mathrm{A}$ is “fat”), and even if it is,
in practical computation $\hat{\bm{u}}$ is indistinguishable from any
$\hat{\bm{u}}+\bm{h}$ if $\left\|\mathrm{A}\bm{h}\right\|_{2}$ is below the
round-off error level, which, in the case at hand, may happen even if
$\left\|\bm{h}\right\|_{2}$ is substantial.

Both difficulties can be tackled by a regularization technique, whereby (1) is
replaced with

$\displaystyle\hat{\bm{u}}_{\alpha}=\operatorname*{argmin}\limits_{\bm{u}}\psi_{\alpha}(\bm{u}),\quad\psi_{\alpha}(\bm{u})=\frac{1}{2}\left\|\mathrm{A}\bm{u}-\bm{b}\right\|_{2}^{2}+\alpha
R(\bm{u}),$ (2)

where $R(\bm{u})$ is a suitable regularization functional, and $\alpha$ is a
positive scalar parameter. The LS data term in (2) can be replaced by other
fit-to-data functional, for example, Poisson log- likelihood functional, which
is better suited noise model for ET reconstruction [4]. Optimization problem
with regularized Poisson log-likelihood is

$\hat{\bm{u}}_{\alpha,P}=\operatorname*{argmax}\limits_{\bm{u}}\psi_{\alpha,P}(\bm{b},\mathrm{A}\bm{u}),$
$\displaystyle\psi_{\alpha,P}(\bm{b},\mathrm{A}\bm{u})=\sum_{i}^{M}\\{-[\mathrm{A}\bm{u}]_{i}+\bm{b}_{i}\log[\mathrm{A}\bm{u}]_{i}\\}+\alpha
R(\bm{u}).$ (3)

Now we consider various regularization terms $R(\bm{u})$ that can be used in
(2) and (3) cost functions. A prime example of the regularization is known as
Tikhonov’s [15], where

$\displaystyle R_{\ell_{2}}(\bm{u})=\left\|\bm{u}\right\|_{2}^{2}.$ (4)

Regularization makes the minimization problems (2) and (3) well-posed, the
eigenvalues of the Hessian of $\psi_{\alpha}(\bm{u})$ being not less than
$\alpha$. Tikhonov regularization is quadratic, therefore high frequencies
which are related to the object boundaries are penalized more, resulting in a
smooth recovery of $\hat{\bm{u}}$. To preserve boundaries one needs to
consider non-quadratic penalties, e.g. the TV penalty [5] can significantly
improve oscillatory solutions. The differentiable (due to small $\epsilon$
constant) TV penalty is given as:

$\displaystyle
R_{TV}(\bm{u})=\left\||\nabla\bm{u}|_{\epsilon}\right\|_{1}=\left\|\sqrt{\bm{u}_{x}^{2}+\bm{u}_{y}^{2}+\epsilon^{2}}\right\|_{1}.$
(5)

Unlike Tikhonov’s penalty (4), TV can recover function while preserving sharp
discontinuities, however, it also leads to the patchy appearance of the image
[8, 9, 10, 11, 12, 13, 14]. Since TV is unable to recover smooth variations
alone, it can be coupled with $\ell_{2}$-norm based penalty [8, 9, 12, 13]
resulting in the combined functional

$\displaystyle
R_{TV-\ell_{2}}(\bm{u})=\left\||\nabla\bm{u}|_{\epsilon}\right\|_{1}+\frac{\mu}{\alpha}\left\|\bm{u}\right\|_{2}^{2},$
(6)

where $\mu$ is an additional regularization parameter for the quadratic term.
In [8], the following combined functional has been proposed:

$\displaystyle
R_{TV-\ell_{2}}(\bm{u})=\left\||\nabla\bm{u}|_{\epsilon}\right\|_{1}+\frac{\mu}{\alpha}\left\|\frac{(\Delta\bm{u})^{2}}{|\nabla\bm{u}|^{3}}\right\|_{2},$
(7)

where $\Delta\bm{u}=\bm{u}_{xx}+\bm{u}_{yy}$ denotes an isotropic Laplacian.
In this work, we will be using (5) and (7) functionals for comparison with the
proposed penalty.

###  2.3 A note on the regularization error

It is intuitively obvious that the regularization parameter $\alpha$ in (2)
and (3) must not be large. In order to get some further insight into the
issue, let us estimate the difference between $\hat{\bm{u}}_{\alpha}$ and
$\hat{\bm{u}}$ assuming for simplicity that all singular values of
$\mathrm{A}$ are positive.

Let us assume $R(\bm{u})=\|\mathrm{R}\bm{u}\|^{2}_{2}$, where $\mathrm{R}$ is
a square non-degenerate matrix, and denote
$\mathrm{M}=\mathrm{A}^{*}\mathrm{A}$, and
$\mathrm{N}=\mathrm{R}^{*}\mathrm{R}$. In the case at hand,
$\hat{\bm{u}}_{\alpha}$ (where $\alpha$ may be zero) satisfies

$\displaystyle(\mathrm{M}+\alpha\mathrm{N})\hat{\bm{u}}_{\alpha}=\mathrm{A}^{*}\bm{b},$
(8)

which implies the following equation for the regularization error,
$\bm{h}_{\alpha}=\hat{\bm{u}}_{\alpha}-\hat{\bm{u}}$:

$\displaystyle(\mathrm{M}+\alpha\mathrm{N})\bm{h}_{\alpha}=-\alpha\mathrm{N}\hat{\bm{u}}.$
(9)

Multiplication by $\mathrm{M}^{-1}\mathrm{N}\bm{h}_{\alpha}$ yields

$\displaystyle((\mathrm{N}+\alpha\mathrm{N}\mathrm{M}^{-1}\mathrm{N})\bm{h}_{\alpha},\bm{h}_{\alpha})=-\alpha(\mathrm{N}\hat{\bm{u}},\mathrm{M}^{-1}\mathrm{N}\bm{h}_{\alpha}).$
(10)

Now, in the left-hand side of (10) we have

$\displaystyle((\mathrm{N}+\alpha\mathrm{N}\mathrm{M}^{-1}\mathrm{N})\bm{h}_{\alpha},\bm{h}_{\alpha})\geq(\mathrm{N}\bm{h}_{\alpha},\bm{h}_{\alpha})=(\mathrm{R}^{*}\mathrm{R}\bm{h}_{\alpha},\bm{h}_{\alpha})=\left\||\mathrm{R}\bm{h}_{\alpha}\right\|_{2}^{2},$

and in the right-hand side
$(\mathrm{N}\hat{\bm{u}},\mathrm{M}^{-1}\mathrm{N}\bm{h}_{\alpha})=$

$\displaystyle(\mathrm{N}\mathrm{M}^{-1}\mathrm{N}\hat{\bm{u}},\bm{h}_{\alpha})=(\mathrm{R}\mathrm{M}^{-1}\mathrm{N}\hat{\bm{u}},\mathrm{R}\bm{h}_{\alpha})\leq\left\|\mathrm{R}\mathrm{M}^{-1}\mathrm{N}\hat{\bm{u}}\right\|_{2}\left\|\mathrm{R}\bm{h}_{\alpha}\right\|_{2}.$

Thus, (10) implies

$\displaystyle\|\mathrm{R}\bm{h}_{\alpha}\|_{2}\leq\alpha\left\|\mathrm{R}\mathrm{M}^{-1}\mathrm{N}\hat{\bm{u}}\right\|_{2}.$
(11)

###  2.4 Edge-preserving piecewise-smooth regularization

At first glance, the regularization appears to be merely a compromise move
that distorts the problem so that it becomes solvable. While this is certainly
so in the case of Tikhonov’s regularization, an alternative viewpoint can be
offered, which is helpful in designing a proper regularization for the problem
at hand. One observes that the problem (1) can only be solved approximately,
if only because of the inexactness of computer arithmetic. One observes
further that it may have infinitely many approximate solutions that are
indistinguishable in practical computation, as pointed out in the previous
section, if one is only guided by the smallness of the goal data fidelity
functional $\left\|\mathrm{A}\bm{u}-\bm{b}\right\|_{2}^{2}$. The
regularization can be viewed as some kind of additional criterion that helps
to verify whether a particular computed solution is acceptable. This viewpoint
is supported by the fact that in the case where
$R(u)=\left\|\mathrm{R}\bm{u}\right\|_{2}^{2}$, the regularized problem (2) is
equivalent to the original problem (1) for these extended $\mathrm{A}$ and
$\bm{b}$:

$\displaystyle\mathrm{A}_{\alpha}=\left[\begin{array}[]{c}\mathrm{A}\\\ \alpha
R\end{array}\right],\quad\bm{b}_{\alpha}=\left[\begin{array}[]{c}\bm{b}\\\
0\end{array}\right].$ (16)

In the problem (1), $\mathrm{A}$ is a discretization of an integral operator,
owing to which $\left\|\mathrm{A}\bm{h}\right\|_{2}$ is small on oscillating
grid functions $\bm{h}$ with wavelengths that are close to the grid step.
Hence, if one directly applies e.g. CGLS algorithm to the minimization of the
quadratic functional (1), then after sufficiently many iterations one is
likely to end up with a quickly oscillating approximate solution
$\hat{\bm{u}}$. But most images that one encounters in practice do not feature
such oscillations and can be actually represented by piecewise-smooth
functions $\hat{\bm{u}}$. This suggests that the value of $R(\bm{u})$ should
be large on short-wavelength functions $\bm{u}$. At the same time, $R(\bm{u})$
should remain small on the boundaries (walls) between objects constituting the
image, where $\bm{u}$ is discontinuous or has large gradients.

The following regularizer is therefore suggested:

$\displaystyle
R_{EL}(\bm{u})=\left\|\bm{w_{x}}\frac{\partial^{2}\bm{u}}{\partial
x^{2}}\right\|_{2}^{2}+\left\|\bm{w_{y}}\frac{\partial^{2}\bm{u}}{\partial
y^{2}}\right\|_{2}^{2},$ (17)

where the weights $\bm{w_{x}}$ and $\bm{w_{y}}$ are given by

$\displaystyle\bm{w_{x}}=\left(1+\beta\left(\frac{1}{a_{x}}\frac{\partial\bm{u}}{\partial
x}\right)^{2}\right)^{-1},\quad\bm{w_{y}}=\left(1+\beta\left(\frac{1}{a_{y}}\frac{\partial\bm{u}}{\partial
y}\right)^{2}\right)^{-1},$ (18)

$\beta$ is a positive scalar parameter and it determines which points are
considered to be on an edge between objects rather than inside it - the
smaller the value of $\beta$, the smaller the edge area. Constants $a_{x}$ and
$a_{y}$ are the average $x$\- and $y$-derivatives of $\bm{u}$. These averages
are introduced purely for the sake of scale-invariance, and can be computed
e.g. as $a_{x}=2u_{max}/d_{x}$ and $a_{y}=2u_{max}/d_{y}$ where $u_{max}$ is
the maximum of $\bm{u}$ and $d_{x}$ and $d_{y}$ are the sizes of the square
containing the image.

By design, $R_{EL}(\bm{u})$ (where EL stands for the Edge preserving
Laplacian) is large on a short-wavelength functions $\bm{u}$ with small
oscillation amplitude (approaching to isotropic smoothing) inside the
reconstructed objects, and small on the edges between them (owing to the
smallness of weights $\bm{w_{x}}$ and $\bm{w_{y}}$). The edge points are
identified as points where the derivatives of $\bm{u}$ are significantly
higher than their average values.

The regularizer (17) shares some similarities with the high-order penalty
proposed for image denoising in [11]:

$\displaystyle\hat{u}=\operatorname*{argmin}\limits_{u}\int_{\Omega}\omega|\Delta
u|d\bm{\mathrm{x}}+\frac{\lambda}{2}\int_{\Omega}(u-u_{0})^{2}d\bm{\mathrm{x}},$
(19)

where $u_{0}$ is a noisy image and weights defined as
$\omega=\left(1+\beta|\nabla G_{\sigma}\ast u_{0}|\right)^{-1},$ where
$G_{\sigma}$ denotes the Gaussian filter with the kernel width $\sigma$ and
$\ast$ is the convolution. There are few differences of the proposed penalty
(17) from (19). For tomographic reconstruction a good initialization
$\bm{u_{0}}$ is not usually available, therefore our regularizer is built on a
different principle of “catching” prominent edges. The proposed penalty (17)
is independent of the smoothed gradient while remain stable to the presence of
noise (we will justify it with our numerical experiments). Moreover, we
incorporate directional high-order smoothing in our term with consideration of
variant weights (second derivatives are weighted independently), in our
penalty $\bm{w_{x}}\neq\bm{w_{y}}$ (cross partial derivatives can be also
added). Overall the penalty (17) is more rigorous than (19) and it will not
allow any harmonic noise into the solution.

###  2.5 Regularized reconstruction with LS data fidelity term

In order to simplify the computation of the gradient of
$\psi_{\alpha}(\bm{u})$ and dealing with the minimization of a quadratic
functional, we resort to an inner-outer iterative scheme with lagged
diffusivity fixed point iteration [15], with $\bm{w_{x}}$ and $\bm{w_{y}}$
only updated on restarts (this approach helps to deal with non-convexity of
penalty (17) cf. widely used trust region technique [17]).

Let us denote by $\mathrm{L}_{x}$ and $\mathrm{L}_{y}$ the matrices
representing the discretized second partial $x$\- and $y$-derivatives with
Neumann boundary conditions, and by $\mathrm{W}_{x}$ and $\mathrm{W}_{y}$ the
diagonal matrices representing $\bm{w_{x}}$ and $\bm{w_{y}}$. The symmetric
gradient matrix $\mathrm{R}(\bm{u})$, such as $\nabla
R(\bm{u})=\mathrm{R}(\bm{u})\bm{u}$, is given for EL term as (ignoring the
dependence of $\bm{w_{x}}$ and $\bm{w_{y}}$ on $\bm{u}$ – cf. inner-outer
iterative scheme)

$\displaystyle\mathrm{R}_{EL}=\mathrm{L}_{x}^{*}\mathrm{W}_{x}^{*}\mathrm{W}_{x}\mathrm{L}_{x}+\mathrm{L}_{y}^{*}\mathrm{W}_{y}^{*}\mathrm{W}_{y}\mathrm{L}_{y}.$
(20)

We compare the proposed EL penalty with TV regularization term (5) with the
following gradient matrix $\mathrm{R}_{TV}$:

$\displaystyle\mathrm{R}_{TV}=\mathrm{D}_{x}^{*}\Phi(\bm{u})\mathrm{D}_{x}+\mathrm{D}_{y}^{*}\Phi(\bm{u})\mathrm{D}_{y},$
(21)

where $\mathrm{D}_{x}$ and $\mathrm{D}_{y}$ are matrices representing the
discretized first partial $x$\- and $y$-derivatives with Neumann boundary
conditions, $\Phi(\bm{u})=\mathit{diag}(\phi^{\prime}(\bm{u}))$ is a diagonal
matrix whose diagonal elements are $\phi^{\prime}(\bm{u})$,
$\phi(t)=2\sqrt{t+\epsilon^{2}}$. Different choices can be used for $\phi(t)$
function to approximate $\ell_{1}$ norm, for example the Huber function [15,
21].

We also compare the proposed penalty (17) with TV-$\ell_{2}$ functional with
the gradient matrix defined as

$\displaystyle\mathrm{R}_{TV-\ell_{2}}=\mathrm{D}_{x}^{*}\Psi(\bm{u})\mathrm{D}_{x}+\mathrm{D}_{y}^{*}\Psi(\bm{u})\mathrm{D}_{y}+\mathrm{L}_{x}^{*}\Upsilon(\bm{u})\mathrm{L}_{x}+\mathrm{L}_{y}^{*}\Upsilon(\bm{u})\mathrm{L}_{y},$
(22)

where $\Psi(\bm{u})$ is a diagonal matrix with elements
$\alpha/|\nabla\bm{u}|_{\epsilon}$ and $\Upsilon(\bm{u})$ is a diagonal matrix
with elements $2\mu/|\nabla\bm{u}|^{3}_{\gamma}$. Values $\epsilon$ and
$\gamma$ were chosen similar to [8], $\epsilon=10^{-5}u_{max}$ and
$\gamma=u_{max}^{2}$.

Here we used the lagged diffusivity fixed point iteration [15] (see algorithm
1) to solve regularized LS optimization problem with penalties (20,21) and
(22).

Algorithm 1 Lagged Diffusivity Fixed Point Method for Regularized LS

$\nu:=0$

$\bm{u_{0}}:=$ initialization

begin outer (fixed point) iterations

$\bm{g}_{\nu}:=\mathrm{A}^{*}(\mathrm{A}\bm{u}-\bm{b})+\alpha\mathrm{R}\bm{u};$
$\%$ gradient

$\mathrm{H}:=\mathrm{A}^{*}\mathrm{A}+\alpha\mathrm{R};$ $\%$ approximate
Hessian

$\bm{s}_{\nu}:=-\mathrm{H}^{-1}\bm{g}_{\nu};$ $\%$ quasi-Newton step

$\bm{u}_{\nu+1}:=\bm{u}_{\nu}+\bm{s}_{\nu};$ $\%$ update approximate solution

$\nu=\nu+1;$

end outer (fixed point) iterations

The system $\mathrm{H}\bm{s}_{\nu}=\bm{g}_{\nu}$ is solved by usual CG method
[17], let $l$ be an iteration index and $L$ is the maximum number of inner
iterations for CG method. Inner and outer iterations can be terminated earlier
if the following stopping criteria met
$\left\|\bm{s}_{l}-\bm{s}_{l+1}\right\|_{2}^{2}\leq\rho$ and
$\left\|\bm{u}_{\nu}-\bm{u}_{\nu+1}\right\|_{2}^{2}\leq\rho$ respectively. We
chose $\rho=10^{-4}$ to be a small constant.

The largest eigenvalue of the matrix $\mathrm{R}$ is of the order ${\cal
O}(h^{-2})$, $h=\min(h_{x},h_{y})$, and the smallest are of the order ${\cal
O}(1)$. If $\alpha$ is not very small (considerably larger than $h^{4}$), then
the large condition number of $\mathrm{R}^{*}\mathrm{R}$ can slow down the
convergence of CG iterations for the minimization of $\psi_{\alpha}(\bm{u})$.
To alleviate this problem, we introduce preconditioning that consists in the
multiplication of the gradient of $\psi_{\alpha}(\bm{u})$ by the inverse of
$\mathrm{H}=\sigma^{2}\mathrm{I}+\alpha\mathrm{R}$ in the course of CG
iterations, where $\sigma$ is a scalar value of the order of the largest
singular value of $\mathrm{A}$. Since matrix $\mathrm{H}_{\sigma}$ is very
sparse, the application of its inverse to a vector can be efficiently
performed (via the factorization of $\mathrm{H}_{\sigma}$) by modern state-of-
the art sparse direct solvers [18].

###  2.6 Regularized reconstruction with Poisson data fidelity term

To solve regularized Poisson log-likelihood problem (3) we use the splitting
technique introduced in [20] and used for SPECT reconstruction in [21]. The
main idea is to decouple data-fidelity term from the regularizer by solving
two problems independently (see algorithm 2).

Algorithm 2 Splitting algorithm for regularized Poisson term

$\nu:=0$

$\bm{u_{0}}:=$ initialization with ones

begin outer iterations

$\bm{u}_{\nu+\frac{1}{2}}:=\bm{u}_{\nu}./(\mathrm{A}^{*}\bm{1}).*\mathrm{A}^{*}(\bm{b}./\mathrm{A}\bm{u}_{\nu});$
$\%$ MLEM step

$\bm{f_{0}}:=\bm{u}_{\nu+\frac{1}{2}};$

$l:=0;$

$\bm{f_{l+1}}:=\bm{f_{l}}+\tau((\bm{f_{l}}-\bm{f_{0}})+\mathrm{R}\bm{f_{l}});$
$\%$ image denoising

$\bm{u}_{\nu+1}:=\bm{f_{l+1}};$ $\%$ update approximate solution

$\nu=\nu+1;$

end outer iterations

From the structure of algorithm 2, one can see that the optimization problem
with Poisson log-likelihood term is solved independently with Maximum
Likelihood Expectation Maximization (MLEM) algorithm [15]. Regularization is
performed as an additional image denoising step. In MLEM step, $\bm{1}$
denotes the vector of 1’s, and $.*$ and $./$ denote componentwise
multiplication and division, respectively. Additionally the weights
$\bm{w_{x}}$ and $\bm{w_{y}}$ for the proposed term (20) are calculated after
each MLEM step and kept fixed in image denoising updates.

##  3 Numerical Results

In this section we present two different numerical experiments with the
proposed penalty EL (20) and compare it with TV (21) and TV-$\ell_{2}$ (22)
penalties. In the first experiment we model reconstruction of the synthetic
phantom with algorithm 1 and in the second experiment we perform
reconstruction of more realistic phantom for emission tomography with
algorithm 2. The MATLAB code for tomographic reconstruction using TV (21)
regularization and the proposed EL penalty (20) is provided [19].

###  3.1 CT reconstruction

To test the proposed penalty in regularized tomographic reconstruction we
designed an analytical phantom which consists of a smooth (two Gaussians and
two parabolas) and piecewise-constant (one rectangle) functions (see Fig. 1).

Figure 1: (a) 2D analytical phantom which consists of 2 parabolas, 2 Gaussians
and a rectangle; (b) a surface representation of the phantom.

To avoid of reconstructing on the same grid where projection data was
generated (so-called reconstruction with “inverse crime” [16]), we used a
higher resolution of the phantom on a $500\times 500$ isotropic pixel grid to
generate projections with a strip kernel [1]. Then Poisson distributed noise
was applied to the projection data, assuming an incoming beam intensity of
3$\cdot 10^{5}$ (photon count). Reconstructions were calculated on a
$250\times 250$ isotropic pixel grid and with a linear projection model [1].
We used 90 projection angles in 180 degrees (assuming a parallel beam
geometry) to reconstruct the phantom.

For a fair comparison of different regularizing penalties we initially
optimized the regularization parameters (see Fig. 2) with respect to the value
of root-mean-square-error (RMSE), defined as:

$\textrm{RMSE}(\overline{\bm{u}},\hat{\bm{u}})=\frac{\left\|\hat{\bm{u}}-\overline{\bm{u}}\right\|_{2}}{\left\|\overline{\bm{u}}\right\|_{2}},$
(23)

where $\overline{\bm{u}}$ is an exact image and $\hat{\bm{u}}$ is a
reconstructed image.

Figure 2: Optimization procedure to find the optimal regularization parameters
for TV (left), TV-$\ell_{2}$ (middle) and EL (right) algorithms. For
TV-$\ell_{2}$ penalty we optimized the second regularization constant ($\mu$)
while keeping the optimal $\alpha$ fixed (found for TV previously).

We found empirically that $\beta=0.03$ for EL penalty (20) gives good results
for the presented experiments, therefore we will keep it fixed for the rest of
our tests. With fixed optimal regularization parameters (see Fig. 2) we
perform $L=80$ outer (fixed point) iterations and 5 inner iterations of
algorithm 1 with different penalties (see Fig. 3).

Figure 3: Plot of RMSE values to the number of iterations for CGLS, CGLS-TV,
CGLS-TV-$\ell_{2}$ and CGLS-EL methods. CGLS-EL method outperforms other
algorithms (see table 1).

In Fig. 3 one can see that the CGLS method diverges quickly while regularized
methods converge to a steady-state solution. The lowest RMSE value is achieved
with the proposed EL penalty (see table 1). Reconstruction with TV penalty
gives the highest RMSE value (from all compared regularized methods) however
it shows faster convergence on first fixed point iterations. TV-$\ell_{2}$
penalty performs slightly better than TV, but still with higher RMSE than EL
method.

Table 1: RMSE values for CGLS, CGLS-TV, CGLS-TV-$\ell_{2}$ and CGLS-EL methods | CGLS | CGLS-TV | CGLS-TV-$\ell_{2}$ | CGLS-EL
---|---|---|---|---
RMSE | 0.1710 | 0.0919 | 0.0902 | 0.0868

Reconstructed images are presented in Fig. 4. Since CGLS-TV-$\ell_{2}$
reconstruction might look more appealing than CGLS-EL we also show the surface
representations of reconstructed images (see Fig. 5) and horizontal middle
cross-sections (see Fig. 6).

Figure 4: Reconstructed phantoms with (a) CGLS, (b) CGLS-TV, (c) CGLS-
TV-$\ell_{2}$ and (d) CGLS-EL method.

Figure 5: Surface representations of the reconstructed phantoms with (a) CGLS,
(b) CGLS-TV, (c) CGLS-TV-$\ell_{2}$ and (d) CGLS-EL method.

Figure 6: Horizontal middle cross-sections of the phantom and reconstructed
images (see Fig. 4). Left side images show the region selected and the right
side images show the magnified region. It can be clearly seen that the EL
penalty performs much better for smooth objects while slightly more perturbed
in uniform areas (e.g. the top of the rectangle).

One can notice that CGLS reconstruction is very noisy. CGLS-TV method better
suppresses noise, however smooth features are strongly affected by the
“staircasing” effect. CGLS-TV-$\ell_{2}$ method provides reconstruction with
smoother features and CGLS-EL method resolves smooth features even better
(e.g. cone-shaped parabola). Although CGLS-EL method performs very well for
smooth objects one can notice the wave-like variations of intensity in the
background and also at the top of the rectangle (see Fig. 6). This issue can
be explained by the properties of our regularizer, in contrast to TV, our
penalty does not seek the sparsest solution and does not penalize strongly
(pushing to the constant value) a small intensity perturbations. The EL term
tends to preserve all sharp edges while uniform noise is smoothed
isotropically with the Laplacian. In Fig. 6 one can see that the CGLS-EL
method provides better recovery of smooth features while slightly higher
(compare to TV and TV-$\ell_{2}$) perturbations visible in uniform areas (the
top of the rectangle), however, the edges of the rectangle are defined sharper
with the EL penalty.

###  3.2 ET reconstruction

To simulate emission tomography reconstruction we designed a more realistic
phantom from the high-quality X-ray scan of a mice bone. The data was acquired
on a Nikon Metris Custom Bay cone-beam scanner at the Henry Moseley Manchester
X-ray facility, and was reconstructed with the Feldkamp algorithm (see Fig. 7
(left)). We thresholded the obtained reconstruction and added six gaussians
with various kernel widths (see Fig. 7 (middle and right)).

Figure 7: (left) High quality X-ray reconstruction of a mice bone; (middle)
phantom created from the reconstruction (left) with a six Gaussians added to
it; (right) a surface representation of the phantom.

To simulate PET projection data we used NiftyRec [22], a software for
tomographic reconstruction, providing GPU-accelerated reconstruction for
emission and transmission computed tomography. The phantom size is $400\times
400$ pixels and 300 projections was simulated. Poisson noise was added to
projections with an expected number of $10\cdot 10^{6}$ photon counts in
total. Twenty noise realizations were simulated to estimate methods
quantitatively. The point spread function of the PET system was modelled (with
convolution of the sinogram columns with a Gaussian of full width half maximum
of three pixels) in the projection and back-projection operations. No scatter
was simulated in this study. For our experiments (see algorithm 2) we
performed 130 MLEM iterations and 5 inner iterations (denoising step).

Figure 8: Optimization procedure to find the optimal regularization parameters
for MLEM-TV (left), MLEM-TV-$\ell_{2}$ (middle) and MLEM-EL (right)
algorithms. For MLEM-TV-$\ell_{2}$ method we optimized the second
regularization constant ($\mu$) while kept the optimal $\alpha$ fixed (found
for MLEM-TV method previously).

To quantify obtained reconstructions we used averaged over all noise
realizations RMSE (23) values in the bone region (BR) and in Gaussian regions
(GR). All regularization parameters were carefully selected by comparing the
mean of all RMSE values over all noise realizations in GR and BR (see Fig. 8).

After estimation of regularization parameters we performed twenty
reconstructions for each method with various Poisson noise distributions. The
mean values for GR and BR over all noise realizations are shown in Fig. 9.
This result proves that the EL penalty is very successful in resolving smooth
features (six Gaussians in this case) and also quite competitive for the BR
(lower RMSE value than for TV).

Figure 9: Plot of the best RMSE values for MLEM, MLEM-TV, MLEM-TV-$\ell_{2}$
and MLEM-EL methods. MLEM-EL method outperforms other algorithms for GR and
gives better result than TV for BR.

In Fig. 10 and 11 we demonstrate reconstructed images for the best selected
RMSE values.

Figure 10: Reconstructed phantoms with (a) MLEM, (b) MLEM-TV, (c) MLEM-
TV-$\ell_{2}$ and (d) MLEM-EL method.

Figure 11: Surface representations of the reconstructed phantoms with (a)
MLEM, (b) MLEM-TV, (c) MLEM-TV-$\ell_{2}$ and (d) MLEM-EL method.

In Fig. 10 and 11 one can notice that the BR is very smooth for TV and
TV-$\ell_{2}$ penalties and some long-wave oscillations can be seen in the
reconstructed image with EL penalty. This result corresponds to the expected
behaviour of the EL penalty. We note here that the phantoms background (see
Fig. 7) is not as flat as TV and TV-$\ell_{2}$ penalty recovered it.
Furthermore, a small size dot-like feature (approximately in the centre of the
phantom) is almost smoothed out with TV and TV-$\ell_{2}$ recovery. However,
it is visible and well recovered with EL penalty. The sharp features, overall,
are reconstructed very well with MLEM-EL method and seem even sharper compare
to other methods (see the bone outer rim in Fig 10).

##  4 Discussion

In this paper, we emphasize that piecewise-smooth images are more favourable
than piecewise-constant ones, this, however might not be the case for all
reconstructed objects. Therefore, some prior knowledge about the investigated
object is needed to choose an appropriate regularization term. For instance,
the activity distribution in ET is smooth, therefore, penalties like EL are
particularly suitable. The proposed penalty gives more gentle approach to
regularization in uniform areas and does not penalize small perturbations
aggressively. Moreover, the proposed method can recover small features
successfully (e.g. lesions in ET case) with both step or ramp intensity
variations.

In terms of the choice of parameters, complexity of computer implementation
and the speed of computation, our method is very similar to reconstruction
techniques with TV penalty. Our future work will be to explore further the
space of parameters and give some recommendations for automated choice of some
parameters. Additionally we will look into the issue of low-frequency
oscillations of our penalty in the uniform areas.

##  5 Conclusion

In this paper, we presented a novel two-step iterative reconstruction
algorithm with high-order regularization penalty. Our method outperforms in
terms of signal-to-noise ratio the conventional total variation based
reconstruction techniques and it is competitive with other state-of-the-art
high-order based approaches. From the preliminary experiments, the proposed
method is well suited for the limited data problems in X-ray tomography as
well as emission tomography.

This work has been supported by the Engineering and Physical Sciences Research
Council under grants EP/J010553/1, EP/J010456/1 and EP/I02249X/1.

## References

## References

  * [1] Kak A C and Slaney M 2009 Principles of computerized tomographic imaging (New York: IEEE Press)
  * [2] Wernick M N and Aarsvold J N 2004 Emission tomography: the fundamentals of PET and SPECT (Academic Press)
  * [3] Kazantsev D and Van Eyndhoven G and Lionheart W R B and Withers P J and Dobson K J and McDonald S A and Atwood R and Lee P D 2015 Employing temporal self-similarity across the entire time domain in computed tomography reconstruction Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 373(2043) 20140389
  * [4] Qi J and Leahy R M 2006 Iterative reconstruction techniques in emission computed tomography Phys Med. Biol. 51 541–578
  * [5] Rudin L and Osher S and Fatemi E 1992 Nonlinear total variation based noise removal algorithms Physica D 60(1-4) 259–268
  * [6] Wang Y and Yang J and Yin W and Zhang Y 2008 A new alternating minimization algorithm for total variation image reconstruction SIAM Journal on Imaging Sciences 1(3) 248–272
  * [7] Mumcuoglu E Ü and Leahy R M and Cherry S R 1996 Bayesian reconstruction of PET images: methodology and performance analysis Physics in medicine and Biology 41(9) 1777
  * [8] Chan T and Marquina A and Mulet P 2000 High-order total variation-based image restoration SIAM Journal on Scientific Computing 22(2) 503–516
  * [9] Chan R H and Liang H and Wei S and Nikolova M and Tai X C 2015 High-order Total Variation Regularization Approach for Axially Symmetric Object Tomography from a Single Radiograph Inverse Problems and Imaging 9(1) 55–77
  * [10] Bredies K and Kunisch K and Pock T 2010 Total generalized variation SIAM Journal on Imaging Sciences 3(3) 492–526
  * [11] Liu X and Huang L and Guo Z 2011 Adaptive fourth-order partial differential equation filter for image denoising Applied Mathematics Letters 24(8) 1282–1288
  * [12] Papafitsoros K and Schönlieb C-B 2014 A combined first and second order variational approach for image reconstruction Journal of Mathematical Imaging and Vision 48 (2) 308–338
  * [13] Do S and Karl W C and Kalra M K and Brady T J and Pien H 2010 Clinical low dose CT image reconstruction using high-order total variation techniques In SPIE Medical Imaging, International Society for Optics and Photonics 76225D-76225D
  * [14] Lysaker M and Lundervold A and Tai X C 2003 Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time Image Processing, IEEE Transactions on 12(12) 1579–1590.
  * [15] Vogel C R 2002 Computational methods for inverse problems (Vol. 23) (SIAM)
  * [16] Hansen P C 2010 Discrete inverse problems: insight and algorithms 7 (SIAM)
  * [17] Nocedal J and Wright S 2006 Numerical Optimization (Springer)
  * [18] Davis T A 2006 Direct methods for sparse linear systems (Vol. 2) (SIAM)
  * [19] Kazantsev D, 2015, http://ccpforge.cse.rl.ac.uk/gf/project/ccpi_itr/frs/?action=FrsReleaseView&release_id=468 Matlab implementation of tomographic reconstruction with TV and EL penalties
  * [20] Sawatzky A and Brune C and Wubbeling F and Kosters T and Schafers K and Burger M 2008 Accurate EM-TV algorithm in PET with low SNR In Nuclear Science Symposium Conference Record (NSS/MIC) 5133–5137
  * [21] Kazantsev D and Arridge S R and Pedemonte S and Bousse A and Erlandsson K and Hutton B F and Ourselin S 2012, An anatomically driven anisotropic diffusion filtering method for 3D SPECT reconstruction Physics in medicine and biology 57(12) 3793
  * [22] Pedemonte S and Bousse A and Erlandsson K and Modat M and Arridge S and Hutton B F and Ourselin S 2010 GPU accelerated rotation-based emission tomography reconstruction In Nuclear Science Symposium Conference Record (NSS/MIC) 2657–2661


# Harmonic Control Lyapunov Barrier Functions for Constrained Optimal Control
with Reach-Avoid Specifications

Amartya Mukherjee1, Ruikun Zhou1, and Jun Liu1 1Amartya Mukherjee, Ruikun
Zhou, and Jun Liu are with the Department of Applied Mathematics, University
of Waterloo, Waterloo, Ontario, Canada N2L 3G1 (email: a29mukhe, ruikun.zhou,
j.liu@uwaterloo.ca).

###### Abstract

This paper introduces harmonic control Lyapunov barrier functions (harmonic
CLBF) that aid in constrained control problems such as reach-avoid problems.
Harmonic CLBFs exploit the maximum principle that harmonic functions satisfy
to encode the properties of control Lyapunov barrier functions (CLBFs). As a
result, they can be initiated at the start of an experiment rather than
trained based on sample trajectories. The control inputs are selected to
maximize the inner product of the system dynamics with the steepest descent
direction of the harmonic CLBF. Numerical results are presented with four
different systems under different reach-avoid environments. Harmonic CLBFs
show a significantly low risk of entering unsafe regions and a high
probability of entering the goal region.

##  1 Introduction

For safety-critical systems, ensuring stability and safety is of paramount
importance. One approach to address this issue for such systems is to model it
as a constrained optimal control problem. The objective of constrained optimal
control is to design a control system that optimizes a notion of performance
while satisfying a set of constraints. These constraints are usually designed
to prevent catastrophic events, ensure safety, and take into multiple
objectives with bounds. Furthermore, with experiments conducted in real life,
resetting the environment if the agent goes outside a domain can be expensive
[11]. One of the most common scenarios for safety-critical systems is the
reach-avoid problem, which can be regarded as a subset of constrained optimal
control problems where a trajectory aims to reach a goal state while avoiding
a set of unsafe states.

In recent years, there has been a growing interest in reach-avoid problems
among the control theory and reinforcement learning community. Such methods
involve learning control Lyapunov functions to certify reachability and
control barrier functions that satisfy avoidability [9, 3]. In reinforcement
learning, agents learn optimal policies through exploration of the environment
with unknown dynamics [19]. Constraints are imposed using penalty methods [10]
and Lagrangian methods [1]. In these methods, two value functions are computed
to estimate cumulative rewards and cumulative costs.

The problem with using two certificates or two value functions is that it
leads to conflicts in the control policy whether to satisfy reachability or
avoidability at a given point in time. To unify the conflicts, reachability
and safety have been combined into a single certificate called the control
Lyapunov barrier function (CLBF). It is shown in [16] that theoretically there
exists a single Lyapunov barrier function for such reach-avoid type
specifications with stability and safety guarantees, provided that such
specifications are robustly satisfiable. Moreover, [5] introduces robust CLBFs
and derives a QP-based controller, in which the CLBF is computed using neural
networks by encoding the constraints in their loss function. In a similar
manner but with reinforcement learning, the Lyapunov barrier-based actor-
critic method is proposed in [6] to find a controller that satisfies both
certificates with a CLBF that is obtained by solving the Lagrangian relaxation
of the constrained optimization problem. Unfortunately, both methods do not
provide guarantees that their value function $V$ satisfies all the CLBF
constraints, and this is evidenced by their experiments. For example, in 2D
quadrotor control problems, the contour plots of the derived CLBF do not
completely cover the unsafe regions. Our work addresses this issue by
enforcing the CLBF constraints as boundary conditions for harmonic functions.
Numerical results demonstrate that this technique indeed improves safety
rates, compared with the results reported in [5].

Harmonic functions are solutions to the Laplace equation, which is a second-
order linear elliptic partial differential equation (PDE). They satisfy the
maximum principle on a compact set, thus making it easier to impose CLBF
constraints. They do not have any critical points other than saddle points in
the interior of their domain, which makes it easy to derive optimal control
strategies. Harmonic functions have been used in the control literature to
derive potential fields in past works. In the work of [12], harmonic potential
functions with panel methods were used to build potential fields over a task
space. The work of [18] used harmonic functions with complex variables and
conformal mappings to achieve moving obstacle avoidance. Furthermore, [4]
combines harmonic functions using superposition to derive potential fields.

In this paper, we will focus on reach-avoid problems, where an agent must
reach a goal region while avoiding unsafe regions. For example, the unsafe
regions could be regions on a floor that have holes or pillars. We achieve
this by introducing harmonic control Lyapunov barrier functions (harmonic
CLBFs) that use Laplace’s equation to encode the properties of control
Lyapunov barrier functions (CLBFs). Furthermore, Laplace’s equation can be
solved using numerical methods such as finite element methods (FEMs) and do
not be trained using neural networks based on sample trajectories as done in
prior work, thus significantly reducing the computational cost of reach-avoid
experiments.

In all, the main contributions of this paper are threefold. Firstly, we
propose a novel method to unify harmonic functions with CLBFs to provide both
stability and safety guarantees. Secondly, under this framework, we show that
optimal controllers for reach-avoid tasks can be derived using gradient-based
methods directly with no need for searching such controllers by training a
neural network or solving optimal problems. Lastly, we conduct a set of
numerical experiments that show a significantly low risk of trajectories
entering unsafe regions with the proposed approach.

##  2 Preliminaries

Throughout this work, we consider a nonlinear control system of the form

$\dot{x}=f(x,u),\quad x(0)=x_{0},$ (1)

where $x\in\mathcal{S}\subset\mathbb{R}^{n}$ is the state of the system, and
$u\in U\subset\mathbb{R}^{m}$ is the control input.

Let $\mathcal{S}$ be a compact subspace of $\mathbb{R}^{n}$ denoting the space
of all admissible states in a control problem. Let
$\mathcal{S}_{goal}\subset\mathcal{S}$ be an open set denoting the space of
states that the controller aims to reach, and let $\mathcal{S}_{unsafe}$ be a
compact subset of $\mathcal{S}$ denoting the space of states the controller
aims to avoid. Define
$\mathcal{S}_{safe}=\mathcal{S}\backslash\\{\overline{\mathcal{S}_{goal}}\cup\mathcal{S}_{unsafe}\\}$.

###  2.1 Control Lyapunov Barrier Functions

######  Definition 1 (Control Lyapunov Barrier Function [5, 6]).

A function $V:\mathcal{S}\to\mathbb{R}$ is a control Lyapunov barrier function
(CLBF) if, for some constants $c,\lambda>0$:

  1. 1.

$V(s)=0$ for all $s\in\mathcal{S}_{goal}$;

  2. 2.

$V(s)>0$ for all $s\in\mathcal{S}\backslash\mathcal{S}_{goal}$;

  3. 3.

$V(s)\geq c$ for all $s\in\mathcal{S}_{unsafe}$;

  4. 4.

$V(s)<c$ for all $s\in\mathcal{S}\backslash\mathcal{S}_{unsafe}$;

  5. 5.

There exists a controller $\pi:\,S_{safe}\rightarrow U$ such that
$\langle\nabla V(s),f(s,\pi(s))\rangle+\lambda V(s)\leq 0$ for all $s\in
S_{safe}$.

###  2.2 Harmonic Functions

######  Definition 2 (Harmonic Function [8]).

A harmonic function $u$ is a $C^{2}$-function that satisfies the Laplace
equation $\nabla^{2}u=\nabla\cdot\nabla u=0$

In the following theorems, we introduce properties of harmonic functions that
are of importance in the field of elliptic PDEs, differential geometry, and
complex analysis.

######  Theorem 1 (Mean Value Theorem [17]).

If $u$ is harmonic on a domain $\mathcal{S}$, then $u$ satisfies the mean
value property in $\mathcal{S}$. Furthermore, if $x\in\mathcal{S}$ and $r>0$
are such that $\overline{B_{r}(x)}\subset\mathcal{S}$, where $B_{r}(x)$ is a
sphere of radius $r$ centered at $x$, then

$u(x)=\frac{\int_{B_{r}(x)}u(y)dy}{\int_{B_{r}(x)}dy}.$ (2)

This theorem is essential for proving Theorems 2 and 3 that provide bounds on
harmonic functions.

######  Theorem 2 (Strong Maximum Principle) [8]).

If $u$ is harmonic on a domain $\mathcal{S}$ and $u$ has its maximum in
$\mathcal{S}\backslash\partial\mathcal{S}$, then $u$ is constant.

######  Theorem 3 (Weak Maximum Principle [8]).

If $u\in C^{2}(\mathcal{S})\cup C^{1}(\partial\mathcal{S})$ is harmonic in a
bounded domain $\mathcal{S}$, then

$\max_{\overline{\mathcal{S}}}u=\max_{\partial\mathcal{S}}u.$ (3)

This shows that Theorems 2 and 3 can be exploited to derive a CLBF that has
its maximum in unsafe regions and minimum in goal regions.

##  3 Harmonic CLBF

In this work, we explore the intersection between CLBF and harmonic functions.
By exploiting the maximum principle, we can derive a function
$V:S\to\mathbb{R}$ that satisfies properties (1)–(4) for CLBF by imposing
boundary conditions on the CLBF.

######  Definition 3 (Harmonic CLBF).

A function $V\in C^{2}(\mathcal{S})\cup C^{1}(\partial\mathcal{S})$ is a
harmonic CLBF if it satisfies the following conditions:

  1. 1.

$\nabla^{2}V(s)=0$ for all $s\in\mathcal{S}_{safe}$

  2. 2.

$V(s)=0$ for all $s\in\overline{\mathcal{S}_{goal}}$

  3. 3.

$V(s)=c$ for all $s\in\partial\mathcal{S}\cup\mathcal{S}_{unsafe}$.

In other words, $V(s)$ is a solution to the boundary value problem given
above.

######  Theorem 4.

All harmonic CLBFs $V$ satisfy properties (1)–(4) in Definition 1.

###### Proof.

Properties (1) and (3) in Definition 1 are trivially satisfied as they are set
as boundary conditions from properties (2) and (3) in Definition 3.

If $V(s)\leq 0$ for some $s\in\mathcal{S}\backslash\mathcal{S}_{goal}$
($V(s)\geq c$ for some $s\in\mathcal{S}\backslash\mathcal{S}_{unsafe}$), then
$V(s)$ has its minimum (maximum) in the interior of $\mathcal{S}_{safe}$. By
Theorem 2, $V$ is then constant in $\mathcal{S}_{safe}$, which contradicts the
non-constant boundary conditions imposed on it. Thus, harmonic CLBFs satisfy
properties (2) and (4) in Definition 1. ∎

######  Proposition 1.

Every reach-avoid problem admits a unique harmonic CLBF.

###### Proof.

Consider an arbitrary reach-avoid problem characterized by
$S,S_{goal},S_{unsafe}$. Assume for a contradiction that a reach avoid problem
has two distinct harmonic CLBFs, $V_{1}(s)$ and $V_{2}(s)$.

Define the function $U(s)=V_{1}(s)-V_{2}(s)$. The Laplacian of $U(s)$ is
$\nabla^{2}U=\nabla^{2}V_{1}-\nabla^{2}V_{2}=0$ on $S_{safe}$. As a result,
$U$ is a function satisfying $\nabla^{2}U=0$ on $S_{safe}$, $U=0$ on
$\overline{S_{goal}}$, and $U=0$ on
$\partial\mathcal{S}\cup\mathcal{S}_{unsafe}$. By Theorem 3, $U$ has its
maximum and minimum of $0$ on the boundary, so $U=0$ on $S_{safe}$. Thus, the
harmonic CLBF is unique. ∎

We now introduce a sufficient condition for verifying property (5) of the CLBF
conditions in Definition 1.

######  Theorem 5.

Let $V$ be a harmonic CLBF, and consider the deterministic nonlinear system
(1). If there exists a positive constant $\lambda$ such that

$\lambda\leq-\sup_{x\in\mathcal{S}_{safe}}\inf_{u\in U}\frac{\langle
f(x,u),\nabla V(x)\rangle}{V(x)},$ (4)

then property (5) in Definition 1 holds.

###### Proof.

Define an optimal controller $\pi^{*}(x)$ as

$\pi^{*}(x)=\arg\inf_{u\in U}\langle f(x,u),\nabla V(x)\rangle.$

Since $\lambda V(x)\geq 0$ for all $x\in\mathcal{S}_{safe}$, it follows that

$\displaystyle\frac{\langle f(x,\pi^{*}(x)),\nabla V(x)\rangle}{V(x)}$
$\displaystyle\leq\inf_{u\in U}\frac{\langle f(x,u),\nabla V(x)\rangle}{V(x)}$
$\displaystyle\leq\sup_{x\in\mathcal{S}_{safe}}\inf_{u\in U}\frac{\langle
f(x,u),\nabla V(x)\rangle}{V(x)}$ $\displaystyle\leq-\lambda,$

where the last inequality is precisely (4). We have verified property (5) in
Definition 1.

∎

Furthermore, a sufficient condition for a system to avoid
$\mathcal{S}_{unsafe}$ is:

$\sup_{x\in\mathcal{S}_{safe}}\inf_{u\in U}\langle f(x,u),\nabla
V(x)\rangle\leq 0$ (5)

Since $\mathcal{S}_{unsafe}$ is a closed set, this means $V(x)<c$ for all
$x\in\mathcal{S}_{safe}$. Furthermore, since $\frac{\partial}{\partial
t}V(x(t))\leq 0$, this means for all $t\in[0,\infty)$, $x(t)$ will never
approach a point where $V(x)=c$, meaning, it will never approach a point in
$\mathcal{S}_{unsafe}$.

###  3.1 Avoiding the saddle points

As shown in [20], all critical points in $V(x)$ are saddle points. Random
noise is added in [7] to parameters while performing gradient descent to
guarantee that the algorithm converges to a local minimum. The authors of [13]
show that gradient descent with random parameter initialization asymptotically
avoids saddle points. So in this paper, we compare deterministic control
inputs with stochastic control inputs to see which method shows better
convergence to $\mathcal{S}_{goal}$.

At every point $x\in\mathcal{S}_{safe}$, the control input is selected as
$u=\text{argmin}_{u\in U}\langle f(x,u),\nabla V(x)\rangle+z$ where $z$ is a
small noise sampled from a normal distribution $N(0,\sigma^{2})$ to mitigate
local optima.

###  3.2 Conditions on the noise

In order to preserve the properties of the CLBF, we introduce bounds that must
be imposed on $z$ with Taylor series expansion so that the trajectory avoids
unsafe regions when the first-order terms dominate.

Let $\Delta t$ be the time step size used for numerical simulation. Consider
the first-order expansion of $x(t+\Delta t)$ as

$x(t+\Delta t)=x+f(x,u)\Delta t+O(\Delta t^{2}).$ (6)

By adding some noise $z$ to the control input to avoid saddle points, we have

$x(t+\Delta t)=x+f(x,u+z)\Delta t+O(\Delta t^{2}).$ (7)

Then we apply again the first-order expansion of $f(x,u+z)$, yields

$x(t+\Delta t)=x+f(x,u)\Delta t+f_{u}(x,u)^{T}z\Delta t+O(\Delta t^{2}).$ (8)

Also, we consider the first-order expansion of $V(x(t+\Delta t))$:

$\displaystyle V(x(t+\Delta t))=$ $\displaystyle V(x)+\langle\nabla
V,f(x,u)+\nabla_{u}f(x,u)^{T}z\rangle\Delta t$ $\displaystyle+O(\Delta
t^{2}).$ (9)

To have the safety guarantees, the condition $V(x(t+\Delta t))<c$ needs to be
satisfied, that is,

$\langle\nabla V(x),f(x,u)+\nabla_{u}f(x,u)^{T}z\rangle<\frac{c-V(x)}{\Delta
t}.$ (10)

Expanding and simplifying this expression gives the following bounds:

$\langle\nabla_{u}f(x,u)\nabla V(x),z\rangle<\frac{c-V(x)}{\Delta
t}-\langle\nabla V(x),f(x,u)\rangle.$ (11)

Using the Cauchy-Schwarz inequality, a sufficient way to bound the noise $z$
is

$||\nabla_{u}f(x,u)\nabla V(x)||_{2}||z||_{2}<\frac{c-V(x)}{\Delta
t}-\langle\nabla V(x),f(x,u)\rangle.$ (12)

If $f$ is a control-affine function $f(x,u)=f_{1}(x)+f_{2}(x)u$, then we can
bound $z$ as:

$||z||_{2}<\frac{\frac{c-V(x)}{\Delta t}-\langle\nabla
V(x),f(x,u)\rangle}{||f_{2}(x)||_{2}||\nabla V(x)||_{2}}.$ (13)

Since Theorem 5 requires that under the optimal controller $\pi^{*}(x)$,
$\langle\nabla V(x),f(x,\pi^{*}(x))\rangle\leq 0$, this shows that the bound
imposed on $z$ is positive.

######  Remark 1.

Numerical solutions in section 4 show that in harmonic CLBFs, the distinction
between safe and unsafe regions is unclear. To make the distinction more
transparent, we replace the CLBF as a solution to Laplace’s equation with a
solution to Poisson’s equation $\nabla^{2}V=-6$ in this paper, meaning $V$ is
a superharmonic function. However, this method also poses a risk that $V$ has
local minima in the interior of its domain [2]. This is undesirable as local
minima are harder to escape compared to saddle points. We will compare
harmonic CLBFs with superharmonic CLBFs with numerical results.

##  4 Numerial Results

In this section, we will explore two different reach-avoid environments, one
with four small unsafe regions and the other with two big unsafe regions, for
three different systems: Roomba, DiffDrive, and CarRobot [15]. Except that, we
also perform a reach-avoid task for 2D Quadrotor [5]. The dynamics of each of
these systems are provided in Appendix I along with control inputs that
minimize the expression in Equation (5) for any $x\in\mathcal{S}_{safe}$. We
use $c=1$ in the CLBF defined in Definition 1 for all numerical results.

For each environment and system, we compute the harmonic ($\nabla^{2}V=0$) and
superharmonic ($\nabla^{2}V=-6$) CLBF numerically using finite element
methods. For both of the CLBFs, we compute the trajectories of the system with
1,000 different randomly initialized initial conditions and count the number
of time steps (of size $\Delta t=0.1$) it takes for the system to reach
$\mathcal{S}_{goal}$. We test this setting with both deterministic
($\sigma=0$) and stochastic ($\sigma=0.1$, clipped by the upper bound in
Equation (13)) controllers and report the mean time taken to reach the goal
($\mu_{T}$), the standard deviation of the time taken to reach the goal
($\sigma_{T}$), the number of times the system ends in an unsafe region (not
included in $\mu_{T},\sigma_{T}$), and the number of times the system does not
reach $\mathcal{S}_{goal}$ after 1,000 time steps (not included in
$\mu_{T},\sigma_{T}$).

###  4.1 Problem I

In this problem, we explore an environment that contains a goal region near
the origin and four unsafe regions in the interior of the domain. The results
with three car-like dynamical systems, Roomba, DiffDrive, and CarRobot, are
reported.

$\displaystyle S$ $\displaystyle=[-1,1]\times[-1,1]$ (14)
$\displaystyle\mathcal{S}_{goal}$ $\displaystyle=[-0.1,0.1]\times[-0.1,0.1]$
(15) $\displaystyle\mathcal{S}_{unsafe}$ $\displaystyle=\partial S\cup
C_{1}\cup C_{2}\cup C_{3}\cup C_{4},$ (16)

with the subdomains of the unsafe region given by

$\displaystyle C_{1}$ $\displaystyle=[-0.5,-0.3]\times[-0.5,-0.3]$ (17)
$\displaystyle C_{2}$ $\displaystyle=[-0.5,-0.3]\times[0.3,0.5]$ (18)
$\displaystyle C_{3}$ $\displaystyle=[0.3,0.5]\times[-0.5,-0.3]$ (19)
$\displaystyle C_{4}$ $\displaystyle=[0.3,0.5]\times[0.3,0.5].$ (20)

Initial conditions are sampled as:

$\displaystyle x(0),y(0)$ $\displaystyle\sim U[-0.9,-0.6]\cup[0.6,0.9]$ (21)
$\displaystyle\theta(0)$ $\displaystyle\sim U[0,2\pi].$ (22)

We derive a harmonic CLBF using finite element methods using DOLFIN [14]. We
used piecewise linear trial functions and a triangular mesh. The domain has
been divided into 5,000 triangular meshes of equal area. This numerical
solution is plotted in 2D view in Fig. 1(a).

(a) Problem I with $\nabla^{2}V=0$

(b) Problem I with $\nabla^{2}V=-6$

(c) Problem II with $\nabla^{2}V=0$

(d) Problem II with $\nabla^{2}V=-6$

Figure 1: Contour plots of the CLBFs for problems I and II. (a) & (b) The
plots for problem I with the CLBF obtained by solving Laplace equation (left)
and Poisson’s equation (right) respectively; (c) & (d) The plots for problem
II with the CLBF obtained by solving Laplace equation (left) and Poisson’s
equation (right) respectively.

Although in theory, the solution should take values less than 1.0 outside the
unsafe region, the distinction is highly unclear in the numerical solution. To
clarify the distinction, we replace the Laplace equation with the Poisson
equation $\nabla^{2}V=-6$. We plotted the solution in 2D view in Fig. 1(b).
This solution shows a better distinction between safe and unsafe regions, but
it shows more points where $V(x)$ has a local optimum. Local optima are
undesirable as $\langle f(x,u),\nabla V(x)\rangle=0$ in those points so $V(x)$
is unlikely to show a significant decrease.

System | $\nabla^{2}V$ | $\sigma$ | $\mu_{T}$ | $\sigma_{T}$ | unsafe | no reach
---|---|---|---|---|---|---
Roomba | 0 | 0 | 97.278 | 65.607 | 0 | 105
Roomba | -6 | 0 | 58.224 | 50.388 | 0 | 459
Roomba | 0 | 0.1 | 114.706 | 81.527 | 0 | 0
Roomba | -6 | 0.1 | 139.606 | 180.247 | 0 | 185
DiffDrive | 0 | 0 | 301.403 | 216.495 | 0 | 95
DiffDrive | -6 | 0 | 127.045 | 115.348 | 0 | 492
DiffDrive | 0 | 0.1 | 239.759 | 220.689 | 0 | 2
DiffDrive | -6 | 0.1 | 158.837 | 167.421 | 0 | 5
CarRobot | 0 | 0 | 452.125 | 323.942 | 0 | 193
CarRobot | -6 | 0 | 277.398 | 253.148 | 0 | 47
CarRobot | 0 | 0.1 | 21.200 | 69.889 | 0 | 930
CarRobot | -6 | 0.1 | 26.887 | 88.174 | 0 | 920
Table 1: Results for Problem I. Each system was run with 1,000 different
randomly initialized initial conditions under harmonic and superharmonic CLBFs
with deterministic and stochastic controllers.

The numerical results are listed in Table 1. Roomba achieves the best results
under a harmonic CLBF with a stochastic control policy as all the trajectories
converged to $\mathcal{S}_{goal}$. DiffDrive achieves the best results under a
harmonic CLBF with a stochastic control policy with a $99.8\%$ safety rate and
similar results with a superharmonic CLBF with a stochastic control policy
with a $99.5\%$ safety rate. CarRobot achieves the best results under a
superharmonic CLBF with a deterministic policy with a $95.3\%$ safety rate.

###  4.2 Problem II

In this problem, we explore an environment that contains a goal region near
the origin and two unsafe regions in the interior of the domain. In a similar
vein to problem I, the three models are tested with this reach-avoid task.

$\displaystyle S$ $\displaystyle=[-1,1]\times[-1,1]$ (23)
$\displaystyle\mathcal{S}_{goal}$ $\displaystyle=[-0.1,0.1]\times[-0.1,0.1]$
(24) $\displaystyle\mathcal{S}_{unsafe}$ $\displaystyle=\partial S\cup
C_{5}\cup C_{6},$ (25)

with the subdomains of the unsafe region given by

$\displaystyle C_{5}$ $\displaystyle=[-0.5,-0.3]\times[-0.5,0.5]$ (26)
$\displaystyle C_{6}$ $\displaystyle=[0.3,0.5]\times[-0.5,0.5].$ (27)

Initial conditions are sampled as:

$\displaystyle x(0)$ $\displaystyle\sim U[-0.9,-0.6]\cup[0.6,0.9]$ (28)
$\displaystyle y(0)$ $\displaystyle\sim U[-0.3,0.3]$ (29)
$\displaystyle\theta(0)$ $\displaystyle\sim U[0,2\pi].$ (30)

The numerical solution to the harmonic CLBF for this problem is plotted in 2D
view in Fig. 1(c).

To make the distinction between safe and unsafe regions clearer, we plotted
the solution to the Poisson equation $\nabla^{2}V=-6$ in 2D view in Fig. 1(d).

System | $\nabla^{2}V$ | $\sigma$ | $\mu_{T}$ | $\sigma_{T}$ | unsafe | no reach
---|---|---|---|---|---|---
Roomba | 0 | 0 | 85.600 | 54.603 | 12 | 0
Roomba | -6 | 0 | 202.568 | 199.603 | 0 | 91
Roomba | 0 | 0.1 | 86.846 | 55.789 | 9 | 78
Roomba | -6 | 0.1 | 48.033 | 21.363 | 0 | 880
DiffDrive | 0 | 0 | 258.502 | 168.726 | 0 | 55
DiffDrive | -6 | 0 | 96.077 | 28.211 | 0 | 730
DiffDrive | 0 | 0.1 | 265.348 | 164.137 | 7 | 9
DiffDrive | -6 | 0.1 | 209.452 | 171.338 | 0 | 118
CarRobot | 0 | 0 | 380.014 | 230.751 | 30 | 124
CarRobot | -6 | 0 | 415.246 | 257.179 | 0 | 570
CarRobot | 0 | 0.1 | 264.435 | 205.930 | 45 | 942
CarRobot | -6 | 0.1 | 342.060 | 262.481 | 0 | 987
Table 2: Results for Problem II. Each system was run with 1,000 different
randomly initialized initial conditions under harmonic and superharmonic CLBFs
with deterministic and stochastic controllers.

The numerical results are posted in Table 2. Roomba achieves the best results
under a superharmonic CLBF with a deterministic control policy with a $90.9\%$
safety rate. DiffDrive achieves the best results under a harmonic CLBF with a
deterministic control policy with a $94.5\%$ safety rate. And CarRobot
achieves the best results under a harmonic CLBF with a deterministic policy
with a $84.6\%$ safety rate, although it has 30 trajectories converging to
$\mathcal{S}_{unsafe}$ and 124 trajectories not converging to
$\mathcal{S}_{goal}$. Furthermore, this table shows that deterministic control
policies achieve better results than stochastic control policies and have
significantly more trajectories that converge to $\mathcal{S}_{goal}$.

Tables 1 and 2 show that deterministic policies under superharmonic CLBFs
always avoid $\mathcal{S}_{unsafe}$ as expected since the Poisson equation
makes the distinction between safe and unsafe regions more evident. But it
comes with the drawback that there is a high risk that the trajectory will
converge to a local minimum and not reach $\mathcal{S}_{goal}$. Deterministic
policies generally outperform stochastic policies as the noise added to
control inputs comes with the risk that it drives the trajectory away from the
goal or towards an unsafe region. Despite using control inputs that solve the
minimization problem in Equation (5), in many episodes, the trajectory does
not converge to $\mathcal{S}_{goal}$. This is likely due to numerical errors
associated with estimating $\nabla V$ or with solving the system dynamics.

###  4.3 2D Quadrotor

In this problem, we control a quadrotor to land at a safe ground-level region
while avoiding two obstacles. The dynamics of the quadrotor are in Appendix
5.4. We are running the same problem as in [6] and the horizontally flipped
problem in [5]. The safe and unsafe subsets are given as the following.

$\displaystyle S$ $\displaystyle=[-1,2]\times[0,2]$ (31)
$\displaystyle\mathcal{S}_{goal}$ $\displaystyle=[-1,0.5]\times[0.25,0.75]$
(32) $\displaystyle\mathcal{S}_{unsafe}$ $\displaystyle=\partial S\cup
C_{7}\cup C_{8}\cup C_{9},$ (33)

with the subdomains of the unsafe region given by

$\displaystyle C_{7}$ $\displaystyle=[-1,0]\times[1.25,2]$ (34) $\displaystyle
C_{8}$ $\displaystyle=[-1,2]\times[0,0.25]$ (35) $\displaystyle C_{9}$
$\displaystyle=[0.5,1]\times[0,1].$ (36)

Initial conditions are sampled as:

$\displaystyle x(0)$ $\displaystyle\sim U[1.4,1.6]$ (37) $\displaystyle z(0)$
$\displaystyle\sim U[0.3,1.5]$ (38) $\displaystyle
x^{\prime}(0),z^{\prime}(0),\theta(0),\theta^{\prime}(0)$ $\displaystyle\sim
U[-0.1,0.1].$ (39)

The numerical solution to the Harmonic CLBF for this problem is plotted in
Fig. 2.

Figure 2: Countour plot of the harmonic CLBF for the 2D Quadrotor environment

In this case, we use a neural network deterministic controller
$\pi_{\kappa}(\textbf{x})$ parameterized by $\kappa$, where
$\textbf{x}=(x,\dot{x},z,\dot{z},\theta,\dot{\theta})$, and we train the
controller using model-predictive control (MPC). After numerically solving the
ode $\dot{\textbf{x}}=f(\textbf{x},\pi_{\kappa}(\textbf{x}))$ for a horizon of
128 time steps with time step size $\Delta t=0.01$ using RK4, we compute the
loss as the cosine similarity of
$f(\textbf{x}_{t},\pi_{\kappa}(\textbf{x}_{t}))$ and $\nabla
V(\textbf{x}_{t})$ for all time steps. The time step size is 0.01. We repeat
this for 100,000 different initial conditions.

$\mu_{T}$ | $\sigma_{T}$ | unsafe | no reach | safety rate
---|---|---|---|---
70.350 | 16.114 | 43 | 96 | 86.1%
Table 3: Results for 2D Quadrotor

After training, we test the controller on the environment for 1,000 different
initial conditions. The numerical results are posted in Table 3. The safety
rate of the controller is $86.1\%$, which is higher than the safety rates for
rCLBF-QP ($83\%$) and Robust-MPC ($53\%$) reported in [5]. To compare with
[6], we have computed the trajectories of the quadrotor for the initial
conditions $\textbf{x}_{1}=[1.5,0,0.35,0.0.0]$ and
$\textbf{x}_{2}=[1.5,0,1.2,0.0.0]$. In Fig. 3, we overlay the trajectories on
a superhamornic CLBF plot with $\nabla^{2}V=-6$ that we computed to make the
distinction between safe and unsafe regions clearer. The trajectories in the
plot show that the quadrotor can effectively avoid unsafe regions and enter
the goal regions.

Figure 3: Plot of trajectories on the 2D Quadrotor environment with
$z(0)=0.35$ (blue) and $z(0)=1.2$ (red).

##  5 Conclusions

In this paper, we introduced harmonic CLBFs that exploit the maximum principle
that harmonic functions satisfy to encode the properties of CLBFs. This paper
is the first to unify harmonic functions with CLBFs for an application to
control theory. We select control inputs that maximize the inner product of
the system dynamics with the steepest descent direction of the harmonic CLBF.
This has been applied to reach-avoid problems and demonstrated a low risk of
entering unsafe regions while converging to the goal region.

## Appendix I: Dynamics of each system

###  5.1 Roomba

The dynamics of the Roomba are as follows:

$\displaystyle\dot{x}=v\cos\theta,\quad\dot{y}=v\sin\theta,\quad\dot{\theta}=\omega,$

where $v\in[-1,1]$ and $\omega\in[-1,1]$ are the control inputs.

The minimization problem to Equation (5) is:

$\inf_{v\in[-1,1],\omega\in[-1,1]}(v\cos\theta,v\sin\theta)(V_{x}(x,y),V_{y}(x,y))^{T}.$

The infimum is reached when $v$ maximizes the negative of the magnitude of the
inner product, and $w$ is set to the direction that maximizes the magnitude of
the inner product, which yields

$\displaystyle v$
$\displaystyle=-\text{sign}(V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta),$
$\displaystyle w$
$\displaystyle=\text{sign}(-V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta).$

These control input pairs are used in computing the results in Tables 1 and 2.

###  5.2 DiffDrive

The dynamics of the diff-drive robot are as follows:

$\displaystyle\dot{x}$
$\displaystyle=(u_{L}+u_{R})\frac{r}{2}\cos\theta,\dot{y}=(u_{L}+u_{R})\frac{r}{2}\sin\theta,$
$\displaystyle\dot{\theta}$ $\displaystyle=(u_{R}-u_{L})\frac{r}{2d},$

where $u_{L}\in[-1,1]$ and $u_{R}\in[-1,1]$ are the control inputs subject to
the constraint $|u_{L}|+|u_{R}|\leq 1$. $r=0.1$ and $d=0.1$.

The minimization problem to Equation (5) is:

$\displaystyle\inf_{|u_{L}|+|u_{R}|\leq 1}$
$\displaystyle((u_{L}+u_{R})\frac{r}{2}\cos\theta,(u_{L}+u_{R})\frac{r}{2}\sin\theta)\nabla
V(x,y).$

The infimum is reached when $u_{L}+u_{R}$ maximizes the negative of the
magnitude of the inner product, and $u_{R}-u_{L}$ is set to the direction that
maximizes the magnitude of the inner product, which yields

$\displaystyle u_{R}=[$
$\displaystyle\text{sign}(-V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta)$
$\displaystyle-\text{sign}(V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta)]/2,$
$\displaystyle u_{L}=[$
$\displaystyle-\text{sign}(-V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta)$
$\displaystyle-\text{sign}(V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta)]/2.$

###  5.3 CarRobot

The dynamics of the car-like robot are as follows:

$\displaystyle\dot{x}=v\cos\theta,\quad\dot{y}=v\sin\theta,\quad\dot{\theta}=v\tan(\psi)/l,\quad\dot{\psi}=w,$

where $v\in[-1,1]$ and $w\in[-1,1]$ are the control inputs subject to the
constraint $|v|\leq|w|$. $l=0.1$. $\psi(0)$ is initialized to $0$. We set
$l=0.1$.

The minimization problem to Equation (5) is:

$\inf_{|v|\leq|w|\leq 1}(v\cos\theta,v\sin\theta)(V_{x}(x,y),V_{y}(x,y))^{T}.$

The infimum is reached when $v$ maximizes the negative of the magnitude of the
inner product, and $\theta$ is set to the direction that maximizes the
magnitude of the inner product.

$\displaystyle v$
$\displaystyle=-\text{sign}(V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta),$
$\displaystyle w$
$\displaystyle=v\text{sign}(-V_{x}(x,y)\cos\theta+V_{y}(x,y)\sin\theta)-\text{sign}(\psi).$

###  5.4 2D Quadrotor

We use the implementation of the 2D quadrotor environment by [21]. The
dynamics of the quadrotor are as follows:

$\displaystyle\ddot{x}$ $\displaystyle=\sin\theta(T_{1}+T_{2})/m$
$\displaystyle\ddot{z}$ $\displaystyle=\cos\theta(T_{1}+T_{2})/m-g$
$\displaystyle\ddot{\theta}$ $\displaystyle=(T_{2}-T_{1})d/I_{yy},$

where $T_{1},T_{2}\in[g/2-4,g/2+4]$ are the control inputs. $m=0.033$,
$I_{yy}=1.436e-05$, $g=9.81$ and $d=0.00397/\sqrt{2}$ are constants.

## ACKNOWLEDGMENT

The authors would like to thank Professor Yash Pant from the Department of
Electrical and Computer Engineering, University of Waterloo for providing us
with feedback for this paper.

## References

  * [1] Dario Amodei Alex Ray, Joshua Achiam.  Benchmarking safe exploration in deep reinforcement learning.  OpenAI, 2019.
  * [2] Sheldon Axler, Paul Bourdon, and Ramey Wade.  Harmonic function theory, volume 137.  Springer Science & Business Media, 2013.
  * [3] Jason Choi, Fernando Castaneda, Claire J Tomlin, and Koushil Sreenath.  Reinforcement learning for safety-critical control under model uncertainty, using control lyapunov functions and control barrier functions.  arXiv preprint arXiv:2004.07584, 2020.
  * [4] Christopher I Connolly, J Brian Burns, and Rich Weiss.  Path planning using laplace’s equation.  In Proceedings., IEEE International Conference on Robotics and Automation, pages 2102–2106. IEEE, 1990.
  * [5] Charles Dawson, Zengyi Qin, Sicun Gao, and Chuchu Fan.  Safe nonlinear control using robust neural lyapunov-barrier functions.  In 5th Annual Conference on Robot Learning, 2021.
  * [6] Desong Du, Shaohang Han, Naiming Qi, Haitham Bou Ammar, Jun Wang, and Wei Pan.  Reinforcement learning for safe robot control using control Lyapunov barrier functions, 2023.
  * [7] Rong Ge, Furong Huang, Chi Jin, and Yang Yuan.  Escaping from saddle points — online stochastic gradient for tensor decomposition, 2015.
  * [8] David Gilbarg and Neil S. Trudinger.  Elliptic Partial Differential Equations of Second Order, Second Edition.  Springer-Verlag, 1983.
  * [9] Thomas Gurriet, Andrew Singletary, Jacob Reher, Laurent Ciarletta, Eric Feron, and Aaron Ames.  Towards a framework for realizable safety critical control through active set invariance.  In 2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS), pages 98–106. IEEE, 2018.
  * [10] Nicolas Heess, Dhruva Tb, Srinivasan Sriram, Jay Lemmon, Josh Merel, Greg Wayne, Yuval Tassa, Tom Erez, Ziyu Wang, SM Eslami, et al.  Emergence of locomotion behaviours in rich environments.  arXiv preprint arXiv:1707.02286, 2017.
  * [11] Ching-An Cheng Hoai-An Nguyen.  Provable reset-free reinforcement learning by no-regret reduction.  International Conference on Machine Learning, 2023.
  * [12] Jin-Oh Kim and Pradeep Khosla.  Real-time obstacle avoidance using harmonic potential functions.  IEEE Transactions on Robotics and Automation, June 1992.
  * [13] Jason D. Lee, Max Simchowitz, Michael I. Jordan, and Benjamin Recht.  Gradient descent converges to minimizers, 2016.
  * [14] A. Logg and G. N. Wells.  Dolfin: Automated finite element computing.  ACM Transactions on Mathematical Software 37, 2010.
  * [15] Kevin M. Lynch and Frank C. Park.  Modern Robotics: Mechanics, Planning, and Control.  Cambridge University Press, 2017.
  * [16] Yiming Meng, Yinan Li, Maxwell Fitzsimmons, and Jun Liu.  Smooth converse lyapunov-barrier theorems for asymptotic stability with safety constraints and reach-avoid-stay specifications.  Automatica, 144:110478, 2022.
  * [17] Murray H. Protter and Hans F. Weinberger.  Maximum Principles in Differential Equations.  Springer, 1984.
  * [18] S. Kawamura S. Akishita and K. Hayashi.  Laplace potential for moving obstacle avoidance and approach of a mobile robot.  In I990 Japan-USA Symposium on Flexible Automation, A Pacific Rim Conference, 1990.
  * [19] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [20] A Yanushauskas.  The zeros of the gradient and the hessian of an harmonic function.  Siberian Mathematical Journal, 10(3):497–501, 1969.
  * [21] Zhaocong Yuan, Adam W. Hall, Siqi Zhou, Lukas Brunke, Melissa Greeff, Jacopo Panerati, and Angela P. Schoellig.  Safe-control-gym: A unified benchmark suite for safe learning-based control and reinforcement learning in robotics.  IEEE Robotics and Automation Letters, 7(4):11142–11149, 2022.


# Conformalized Multimodal Uncertainty Regression and Reasoning

###### Abstract

This paper introduces a lightweight uncertainty estimator capable of
predicting multimodal (disjoint) uncertainty bounds by integrating conformal
prediction with a deep-learning regressor. We specifically discuss its
application for visual odometry (VO), where environmental features such as
flying domain symmetries and sensor measurements under ambiguities and
occlusion can result in multimodal uncertainties. Our simulation results show
that uncertainty estimates in our framework adapt sample-wise against
challenging operating conditions such as pronounced noise, limited training
data, and limited parametric size of the prediction model. We also develop a
reasoning framework that leverages these robust uncertainty estimates and
incorporates optical flow-based reasoning to improve prediction prediction
accuracy. Thus, by appropriately accounting for predictive uncertainties of
data-driven learning and closing their estimation loop via rule-based
reasoning, our methodology consistently surpasses conventional deep learning
approaches on all these challenging scenarios–pronounced noise, limited
training data, and limited model size–reducing the prediction error by
2–3$\times$.

Index Terms—  Conformal inference. Visual odometry.

##  1 Introduction

Cutting-edge deep learning frameworks are striving not only for accurate point
predictions but also to express predictive uncertainties. This is achieved by
incorporating uncertainty-aware learning and prediction techniques such as
Bayesian neural networks, Gaussian processes, Monte Carlo dropout, variational
inference, etc. [1, 2, 3, 4]. However, the robustness of learning uncertainty-
aware prediction comes at a considerable computational cost. For example,
methods like Monte Carlo dropout require multiple model samplings to capture
the predictive model’s distribution, where each sampling entails running model
inference. For applications such as autonomous drones, where predictions must
be made in real-time and with limited onboard resources, the computational
challenges for uncertainty-aware predictions thus become excessive [5, 6].

To address these limitations of computationally efficient uncertainty-aware
deep learning, there has been a notable surge of interest in conformal
inference [7, 8, 9, 10, 11]. This approach offers a systematic framework for
estimating predictive uncertainty intervals by incorporating calibration
measures into training procedures. The process entails training a learning
model on a labeled dataset and utilizing calibration data to construct a
predictive uncertainty region around each test instance. In the case of
classification problems, conformal inference adapts them into a minimal set
prediction task, where there is a high confidence that the true class lies
within the prediction set. Likewise, regression problems transform into
minimal interval prediction tasks, ensuring that the true value falls within
the interval with high confidence. Consequently, by employing conformal
prediction, learning models can provide point predictions as well as a
reliable measure of confidence associated with each prediction. A key
computational advantage of conformal inference is its explicit focus on the
uncertainty bounds rather than computing the distribution, which enhances its
computational efficiency.

Fig. 1: Multimodal Uncertainties: Environmental features such as flying domain
symmetries and sensor measurements under ambiguities and occlusion can result
in multimodal (disjoint) uncertainty bands in practical applications such as
visual odometry. We develop lightweight multimodal uncertainty regression and
reasoning to address these practical challenges. Fig. 2: Multimodal
Conformalized Regression: (a) Our framework for embedding conformalization in
a deep learning regressor. (b) Example demonstration of our scheme by
conformalization of softmax scores of multiple interval classes resulting in
multimodal (disjoint) uncertainty interval prediction. Sample demonstration of
our scheme for (c) 10-class and (d) 50-class multi-head configuration.

Yet, a notable drawback of current conformal inference methods for regression
is that they can only predict contiguous (single mode) uncertainty intervals
[10, 12, 13, 14]. Meanwhile, predictive uncertainties can be multi-modal in
many cases. For example, in Fig. 1, consider an autonomous drone navigating a
hallway with multiple similar rooms, doors, and intersections. Determining the
vehicle’s position in such surroundings will likely involve multimodal
uncertainties due to environmental symmetries. Additionally, sensor
measurements can be noisy or prone to occlusions, resulting in multiple
plausible interpretations. Hence, computationally efficient multimodal
uncertainty extraction is crucial.

Addressing these limitations, our work makes the following key contributions:
Using visual odometry (VO) as a driving application, we present a novel
conformalized multimodal uncertainty regression by transforming regression
into a set prediction problem. Results show that the uncertainty estimates
adapt sample-wise against varying operating conditions such as input noise,
limited training data, and the model’s limited parametric size. We also
present multimodal uncertainty reasoning that leverages the robust uncertainty
estimates and closes the estimation loop via rule-based reasoning.
Specifically, under challenging scenarios such as extreme noise, limited
training data, and limited computational constraints, our integrated
methodology consistently outperforms traditional deep learning by a
2–3$\times$ improvement in regression accuracy.

##  2 Multimodal Uncertainty Regression

Overcoming the limitations of current conformalized regression methods [10,
12, 13] and using VO as a driving application, we present a novel
conformalized multimodal uncertainty regression. VO is a predominant computer
vision technique in robotics to estimate the pose of a mounted camera on a
drone/robot [15, 16, 17]. Deep learning-based VO employs neural networks to
automatically learn and extract features from consecutive image frames to
estimate the camera’s pose.

Our method reformulates the learning-based pose regression task as a set
prediction problem. In this context, the predicted set can represent
uncertainty intervals that are not confined to being contiguous. To achieve
this transformation, we extract a calibration set from the initial training
data, enabling the conversion of the pose classification challenge into a set
prediction task through conformalization, as depicted in Fig. 2(b). The first
step involves segmenting the drone’s navigational space into $K$ unique sets,
spanning both position and orientation (or pose) dimensions. This segmentation
is achieved using a non-uniform space discretization, which is informed by the
training set trajectories within each dimension. These trajectories are split
into $K$ quantiles to establish class boundaries. The quantile-based non-
uniform discretization groups frequently visited poses into more confined
spatial intervals. This allows for their identification with enhanced
precision. The outcome of this space encoding process is a one-hot encoded
matrix, contingent on the $K$ classes. In the subsequent stage, a neural
network extracts features from input images. These features are then relayed
to a multi-head classifier, as illustrated in Fig. 2(a). Each classifier head
generates the softmax scores for interval classes along its respective pose
dimension.

Based on the softmax scores, our procedure then targets the sample-adaptive
minimum prediction set $C(X)\subset\\{1,...,K\\}$ such that the correct class
resides within the set with $1-\alpha$ probability, where $\alpha$ denotes an
arbitrary miscoverage rate (e.g., 10%). This marginal coverage can be
expressed as $1-\alpha\leq\mathbb{P}\\{Y\in C(X)\mid X=x\\}\leq
1-\alpha+\frac{1}{n+1}$, where $n$ denotes the size of the calibration set.
Conformal scores arise by deducting the softmax output of the appropriate
class for every input from one. Finally, $\hat{q}$ is computed as the
$\left\lceil{\frac{(n+1)(1-\alpha)}{n}}\right\rceil$ empirical quantile of
these conformal scores, and the conformalized prediction set $C(X)$ is formed
by incorporating classes with softmax scores greater than $(1-\hat{q})$. The
product of the predicted set of classes along each pose dimension produces the
net uncertainty regions, which are not necessarily contiguous.

Fig. 3: Optical flow-based correspondence between consecutive frames.

Figs. 2(c-d) shows the example multimodal uncertainty regions predicted using
the above procedure on the KITTI dataset [18]. The uncertainty regions are
cuboid-shaped due to the product of intervals along respective pose
dimensions. In Fig. 2(d), the uncertainty cuboids contract as the number of
classes $K$ increases from ten to fifty due to higher precision interval
classification. The multimodality and non-contiguity of uncertainty intervals
are evident in the figure. The uncertainty estimates also adapt to drone’s
motion intricacies. For example, the predictive uncertainty increases at sharp
turns.

Fig. 4: Sample Efficiency: In the right panel, the proposed framework
maintains high accuracy despite a significant reduction in training data.
Predictive uncertainties adaptively increase in our framework with lower
training data. [Testcase: RGB-D dataset [19], first test sequence, 50-class
multihead]

##  3 Multimodal Uncertainty Reasoning

###  3.1 Relative Motion Estimation

For reasoning among multimodal uncertainty estimates, our framework determines
the relative motion between two frames using the Harris corner detector [20]
as shown in Fig. 3. The process identifies key feature points, often located
at angles or regions with intensity discontinuities in the first image.
Subsequently, the Lucas-Kanade algorithm [21] pinpoints the corresponding key
points in the second frame. The algorithm presumes minimal movement of
interest points between consecutive frames – which is a reasonable assumption
under a sufficiently high sampling rate and consistent brightness. Under this
assumption, starting from the optical flow equation for each point $(x,y)$:
$I(x+dx,y+dy,t+dt)\approx I(x,y,t)$, a linear approximation of the optical
flux in a neighborhood of the points of interest is computed using a Taylor
series.

The algorithm assumes that the change in brightness of a pixel of the scene is
totally compensated by the gradient of the scene itself, i.e.,
$I_{x}u+I_{y}v+I_{t}=0$. Therefore, the displacement vector $[dx,dy]$ can be
estimated by minimizing the squared error between points in the initial image
and their counterparts in the subsequent frame. Based on the computed optical
flow, points of interest are continuously updated, discarding those with
significant tracking discrepancies, which also yields the key feature points
for the second frame.

Based on the feature correspondence, the relative motion between two poses can
be extracted by exploiting the essential matrix $E$, given by $x_{1}^{T}\cdot
E\cdot x_{0}=0$, where $x_{0}$ is a corresponding point in the first image and
$x_{1}$ is the corresponding point in the destination image. By factoring the
matrix using singular value decomposition (SVD), rotation matrix $R$ and
translation unit vector $t$ can be extracted as in [22]. Notably, essential
matrix-based rotation and translation estimates often prove challenging in
practice due to their sensitivity to small errors in point correspondences,
leading to significant inaccuracies in the derived transformations. Despite
these challenges, we leverage them solely for discerning uncertainty
intervals, i.e., not as the main predictor. As elaborated in the subsequent
section, such integrated prediction combining uncertainty-aware deep learning
and optical flow-based reasoning yields impressive predictive accuracy across
benchmarks, including noise robustness, sample efficiency, and parametric
efficiency.

Fig. 5: Parametric Efficiency: In the left panel, the proposed uncertainty-
aware prediction framework maintains high accuracy despite a smaller network
size of four million parameters. In the right panel, the proposed framework
reduces overfitting risks by uncertainty-oriented predictions. [Testcase:
RGB-D dataset [19], first test sequence, 50-class multihead]

###  3.2 Uncertainty Discrimination

After determining the relative rotation and translation between two frames,
the next step involves discerning among multiple uncertainty intervals to
ascertain the optimal value. We initiate this by computing the mean value
within each uncertainty interval to approximate the mean directions of
displacements. Given that distinct predicted bounds exist for every dimension,
we evaluate all possible combinations to identify the most fitting value. The
algorithm takes the first pose with the highest softmax score from the
prediction model and then iteratively searches for the best pose of the next
frame. To find the best prediction of the orientation corresponding to the
next frame, first, we calculate the corresponding rotation matrix as
$R_{next}=R_{relative}\cdot R_{previous}$, then we transform it into a
quaternion $q_{next}$ and at this point we look for the best prediction by
solving:

$\displaystyle\underset{q_{\text{predicted}}\in Q}{\text{minimize}}\quad\lVert
q_{\text{predicted}}-q_{\text{next}}\rVert$ (1)

where $Q$ is the set of all quaternions describing the orientation of the
camera that took the frame. Likewise, to find the best prediction of the
spatial position corresponding to the next frame, we use the translation
direction $t_{relative}$ calculated in the relative motion estimation phase
and solve:

$\displaystyle\underset{t_{\text{predicted}}\in
T}{\text{minimize}}\quad\left\lVert
t_{\text{relative}}-\frac{t_{\text{predicted}}-t_{\text{previous}}}{\left|t_{\text{predicted}}-t_{\text{previous}}\right|}\right\rVert$
(2)

where $T$ is the set of all translation vectors describing the position of the
camera that took the frame. Solving these optimization steps selects the
uncertainty interval with the highest likelihood for mean prediction and
uncertainty discrimination on successive predictions.

##  4 Simulation Results

In Table 1 and Figs. 4-6, we compare the proposed framework to a conventional
regression model on sample efficiency (i.e., at varying sizes of labeled
training data), parametric efficiency (i.e., the necessary deep learning model
size), and noise robustness. For the comparative study, we only modify the
training procedure from conventional to conformalized uncertainty-based
reasoning while keeping all other settings, such as feature extraction,
training data, etc. the same.

Fig. 6: Noise Robustness: In the right panel, the proposed framework shows
exquisite resilience to extreme noise. Predictive uncertainty levels
appropriately adapt to input image noise levels. [Testcase: RGB-D dataset
[19], first test sequence, 50-class multihead]

Sample Efficiency: Gathering sufficient labeled data for applications like VO
poses significant challenges [23]. Under these practical challenges, Fig. 4
compares the training sample efficiency of the proposed framework against the
conventional. In the right panel, the proposed framework maintains high
accuracy even when the training data reduces to 40% of the original. Notably,
under limited training data, our framework appropriately expresses higher
predictive uncertainties, as seen in left vs. right panels.

Table 1: Comparison of RMSE (Root Mean Square Error) of Pose Prediction Trajectories in Classical Inference vs. Our Uncertainty-Aware Inference. Metric | Conformal | Classical | Improvement
---|---|---|---
40% total data | 0.0534 | 0.136 |  2.5$\times$
80% total data | 0.026 | 0.078 |  3$\times$
Efficientnet-b0 (4M) | 0.029 | 0.083 |  2.9$\times$
Efficientnet-b3 (10.7M) | 0.031 | 0.077 |  2.5$\times$
Efficientnet-b7 (63.8M) | 0.028 | 0.088 |  3.1$\times$
No noise | 0.050 | 0.136 |  2.7$\times$
Gaussian ($\sigma=0.05$) | 0.055 | 0.126 |  2.3$\times$
Gaussian ($\sigma=0.1$) | 0.073 | 0.132 |  1.8$\times$
Gaussian ($\sigma=0.2$) | 0.099 | 0.201 |  2$\times$

Parametric Efficiency: Smaller deep networks are essential for edge computing
of VO to minimize computational and memory demands, enabling real-time
processing on resource-constrained devices. In Fig. 5, we compare the
performance of our framework against conventional deep learning regression at
two variants of feature extractor–EfficientNet-b0 (with four million
parameters) and EfficientNet-b7 (with $\sim$sixty-four million parameters)
[24]. Across these diverse configurations, the performance of the conformal
models remained strikingly consistent despite the substantial variation in
model parameters. While deep learning-based regressors underfit or overfit at
too small or too large a model, in the proposed framework, the parametric
resources are used under the guidance of uncertainty calibration. Our
framework’s mean prediction curves in the left and right panels are
consistent. In contrast, the uncertainty bounds reduce with higher parameters
in the right panel for EfficientNet-b7.

Noise Robustness: Input noise tolerance is crucial for VO to ensure accurate
pose estimation under challenging environmental conditions such as low light
scenarios. In Fig. 6, we compare predictions from both procedures while
injecting Gaussian noise to input image pixels under varying noise variance
($\sigma$). Notably, the proposed framework also shows exquisite tolerance to
input image noise and responds to higher noise levels in input images by
expanding predictive uncertainty intervals.

Table 1 compares the RMSE (Root Mean Square Error) of pose prediction
trajectories on the first testing sequence of RGB-D Scenes dataset [19] in
classical deep learning-based regression against our procedure. By
systematically acknowledging the variabilities in training data, model size,
and input noise, our multimodal uncertainty-aware procedure improves
prediction accuracy by $\sim$2-3$\times$ across all these diverse variations.
When faced with multiple modes of uncertainty (i.e., disjoint uncertainty
intervals), the reasoning framework is able to evaluate each mode’s
implications separately, thus capturing these complexities more accurately
than an unimodal distribution, which otherwise oversimplifies.

##  5 Conclusions

We presented a novel conformalized multimodal uncertainty regressor and
reasoning framework for VO. The proposed framework demonstrated
$\sim$2-3$\times$ lower prediction error than conventional deep learning under
challenging scenarios such as lower training data, model size, and extreme
input noise.

## References

  * [1] E. Goan and C. Fookes, “Bayesian neural networks: An introduction and survey,” _Case Studies in Applied Bayesian Data Science: CIRM Jean-Morlet Chair, Fall 2018_ , pp. 45–87, 2020.
  * [2] E. Schulz, M. Speekenbrink, and A. Krause, “A tutorial on gaussian process regression: Modelling, exploring, and exploiting functions,” _Journal of Mathematical Psychology_ , vol. 85, pp. 1–16, 2018.
  * [3] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_. PMLR, 2016, pp. 1050–1059.
  * [4] 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.
  * [5] P. Shukla, S. Nasrin, N. Darabi, W. Gomes, and A. R. Trivedi, “Mc-cim: Compute-in-memory with monte-carlo dropouts for bayesian edge intelligence,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 70, no. 2, pp. 884–896, 2022.
  * [6] P. Shukla, A. Muralidhar, N. Iliev, T. Tulabandhula, S. B. Fuller, and A. R. Trivedi, “Ultralow-power localization of insect-scale drones: Interplay of probabilistic filtering and compute-in-memory,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 30, no. 1, pp. 68–80, 2021.
  * [7] B.-S. Einbinder, Y. Romano, M. Sesia, and Y. Zhou, “Training uncertainty-aware classifiers with conformalized deep learning,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 22 380–22 395, 2022.
  * [8] H. Karimi and R. Samavi, “Quantifying deep learning model uncertainty in conformal prediction,” _arXiv preprint arXiv:2306.00876_ , 2023.
  * [9] A. M. Alaa, Z. Hussain, and D. Sontag, “Conformalized unconditional quantile regression,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2023, pp. 10 690–10 702.
  * [10] Y. Romano, E. Patterson, and E. Candes, “Conformalized quantile regression,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [11] A. C. Stutts, D. Erricolo, T. Tulabandhula, and A. R. Trivedi, “Lightweight, uncertainty-aware conformalized visual odometry,” _arXiv preprint arXiv:2303.02207_ , 2023.
  * [12] S. Feldman, S. Bates, and Y. Romano, “Improving conditional coverage via orthogonal quantile regression,” _Advances in neural information processing systems_ , vol. 34, pp. 2060–2071, 2021.
  * [13] W. Tang, G. Shen, Y. Lin, and J. Huang, “Nonparametric quantile regression: Non-crossing constraints and conformal prediction,” _arXiv preprint arXiv:2210.10161_ , 2022.
  * [14] K. Stankeviciute, A. M Alaa, and M. van der Schaar, “Conformal time-series forecasting,” _Advances in neural information processing systems_ , vol. 34, pp. 6216–6228, 2021.
  * [15] D. Nistér, O. Naroditsky, and J. Bergen, “Visual odometry,” in _Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004._ , vol. 1. Ieee, 2004, pp. I–I.
  * [16] D. Scaramuzza and F. Fraundorfer, “Visual odometry [tutorial],” _IEEE robotics & automation magazine_, vol. 18, no. 4, pp. 80–92, 2011.
  * [17] M. He, C. Zhu, Q. Huang, B. Ren, and J. Liu, “A review of monocular visual odometry,” _The Visual Computer_ , vol. 36, no. 5, pp. 1053–1065, 2020.
  * [18] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the kitti vision benchmark suite,” _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2012.
  * [19] 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.
  * [20] K. G. Derpanis, “The harris corner detector,” _York University_ , vol. 2, pp. 1–2, 2004.
  * [21] S. Baker and I. Matthews, “Lucas-kanade 20 years on: A unifying framework,” _International journal of computer vision_ , vol. 56, pp. 221–255, 2004.
  * [22] R. I. Hartley and F. Kahl, “Global optimization through searching rotation space and optimal estimation of the essential matrix,” in _2007 IEEE 11th International Conference on Computer Vision_. IEEE, 2007, pp. 1–8.
  * [23] S. Li, X. Wang, Y. Cao, F. Xue, Z. Yan, and H. Zha, “Self-supervised deep visual odometry with online adaptation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 6339–6348.
  * [24] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International conference on machine learning_. PMLR, 2019, pp. 6105–6114.