# Multi-Axis and Multi-Vector Gradient Estimations: Using Multi-Sampled
Complex Unit Vectors to Estimate Gradients of Real Functions

Ergun Akleman
Visual Computing & Computational Media,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
& Alan D. Freed
Department of Mechanical Engineering
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
Joint with Computer Science and Engineering Department.

(October 5, 2023)

###### Abstract

In this preliminary study, we provide two methods for estimating the gradients
of functions of real value. Both methods are built on derivative estimations
that are calculated using the standard method or the Squire-Trapp method
Squire and Trapp (1998) for any given direction. Gradients are computed as the
average of derivatives in uniformly sampled directions. The first method uses
a uniformly distributed set of axes that consists of orthogonal unit vectors
that span the space. The second method only uses a uniformly distributed set
of unit vectors. Both methods essentially minimize the error through an
average of estimations to cancel error terms. Both methods are essentially a
conceptual generalization of the method used to estimate normal fractal
surfaces Hart et al. (1989).

##  1 Introduction

In this work, we study the gradients of general functions of the real value of
$n$ dimensional $F:\Re^{n}\rightarrow\Re$. We will simply write these
functions as $F(\mathbf{p})$, where $\mathbf{p}$ is a point in a
$n$-dimensional space, $\Re^{n}$. In other words, $\mathbf{p}=x$ in 1D,
$\mathbf{p}=(x,y)$ in 2D, $\mathbf{p}=(x,y,z)$ in 3D, etc. Our goal is to
estimate the gradient of $F(\mathbf{p})$ numerically at any given point
$\mathbf{p}_{0}$ in our $n$ dimensional space. We do not expect that the
function has a well-defined derivative. In this short treatment of the
problem, we focus only on the 2D and 3D gradients, but the methods are general
and can be extended to higher dimensions.

We essentially provide two types of estimation: (1) Multi-Axis and (2) Multi-
Vector. Both are conceptual generalizations of the multi-sample method Hart et
al. (1989), which was originally developed to compute normals on fractal
surfaces. These are also closely related to the classical method. Note that
the classical method is based on taking and combining the $n$ number of 1D
derivatives in orthogonal directions for the gradients in the $n$ dimensional
space. In multi-axis gradient estimation, we calculate the gradient in
multiple orthonormal directions and use their average as an estimate. In
multi-vector gradient estimation, we calculate the gradient in multiple
directions and use their average as an estimate.

Both of these methods require a good estimation of the derivatives in a given
direction. To compute derivatives in the given direction, we consider using
both the standard numerical approximation of derivatives, that is, finite
differences, and Squire and Trapp’s Squire and Trapp (1998) method of using
imaginary numbers. These give us six possible ways to compute a gradient: (1)
single-axis sampling using finite differences, (2) single-axis sampling using
imaginary values, (3) multi-axis sampling using finite differences, (4) multi-
axis sampling using imaginary values, (5) multi-vector sampling using finite
differences, and (6) multi-vector sampling using imaginary values. Our initial
qualitative assessment suggests that multi-axis using imaginary values
provides the best results. Unlike the case of a single axis, in the multi-
axis, many equally sampled orthogonal directions can provide a better estimate
of the actual gradient as an average of multiple estimations of gradient
vectors. Multi-vector method also provides an estimate of the average gradient
along a unit circle or a unit sphere (or a unit hypersphere for higher
dimensions) around the given point Hart et al. (1989), however, since these
directions may not include the orthogonal sets of vectors that represent
randomly rotated axis, the error can be higher.

##  2 Preliminaries and Previous Work

In this section, we will briefly discuss existing methods for estimating
derivatives and gradients of functions whose derivatives may or may not exist.
In these methods, it is sufficient that the function can work as a black box
that gives us a single real number for any given point in the domain of the
function. We are only interested in n-dimensional Euclidean domains that are
defined as $\Re^{n}$, where $n$ is the dimension of the domain, which we call
nD such as 1D, 2D, and 3D. In this study, we only discuss the estimation of
gradients in 2D and 3D. Note that 1D gradients are simply derivatives.

###  2.1 Estimation of Derivatives

The classical method for estimating the derivative of a 1D function $F(x)$ at
a point $x_{0}$ is based on the standard central difference formula to define
a derivative

$\left.\frac{\mathrm{d}F(x)}{\mathrm{d}x}\right|_{x_{0}}=\lim_{h\rightarrow
0}\frac{F(x_{0}+h)-F(x_{0}-h)}{2h}$

is as follows:

$\left.\frac{\mathrm{d}F(x)}{\mathrm{d}x}\right|_{x_{0}}\approx\frac{F(x_{0}+h)-F(x_{0}-h)}{2h}$

In theory, the error in this numerical calculation must be $O(h^{2})$.
However, the round-off error that occurs when subtracting the two terms
prevents such accuracy in practice for small values of $h$. In 1998, Squire
and Trapp Squire and Trapp (1998) introduced a new way to improve the
precision of estimation of derivatives of 1D functions with real value using
complex algebra. Squire and Trapp used a Taylor series to derive an
alternative formula involving complex numbers and obtained the following
approximation:

$\left.\frac{\mathrm{d}F(x)}{\mathrm{d}x}\right|_{x_{0}}\approx\frac{\Im(F(x_{0}+ih))}{h}.$

This method symbolically produces the same result as the standard derivative.
To demonstrate this, consider that we have a general quadratic function
$F(x)=ax^{2}+bx+c$, we can ignore the constant term since it will not appear
in the imaginary term, and therefore:

$\displaystyle\frac{dF}{dx}$ $\displaystyle=$
$\displaystyle\frac{1}{h}\Im\left(F(x+ih)\right)$ $\displaystyle=$
$\displaystyle\Im\left(a(x+ih)^{2}+b(x+ih)\right)$ $\displaystyle=$
$\displaystyle\frac{1}{h}\Im\left(ax^{2}+2axih-ah+bx+ibh)\right)$
$\displaystyle=$ $\displaystyle\frac{1}{h}\left(2axh+bh)\right)$
$\displaystyle=$ $\displaystyle 2ax+b$

It is possible to demonstrate that for any function this method symbolically
provides the same derivative as the original derivative. The advantage of this
formula is that it not only gives higher precision digits for small values of
$h$, but the precision stays the same for smaller values of $h$ Greenwell
(1998).

Despite this advantage, this method is useful only for 1D cases. There is a
need for a method that computes the gradient for cases of higher dimensions.
In this study, we generalize Squire and Trapp’s method to handle higher-
dimensional cases.

It is important to note that neither the standard method nor Squire and
Trapp’s imaginary method requires $F(x)$ to have a well-defined derivative.
Even if $F(x)$ does not have a well-defined derivative at a given point, these
methods provide an answer. This is very useful in Computer Graphics.

###  2.2 Multisample Method for Gradient Computation for Fractal Functions

Hart et al. introduced the multi-sample method to render 3D fractals. To
compute diffuse illumination and specular reflection, one needs to estimate
the surface normal in regions around the boundaries of 3D fractals Hart et al.
(1989). The fractals are not represented by nice analytical functions, and the
boundaries have a very high curvature; the estimation is even harder. We
observe that Hart’s solution is to view the fractal description algorithms in
a black-box function that can provide an inside-outside test. Figure 1(a)
shows an example of the boundaries defined by such black-box functions as an
interface between inside and outside. Such a black-box function
$F(\mathbf{p})$ can be given in the following form:

$F(\mathbf{p})=\left\\{\begin{array}[]{ll}a&\mbox{if $\mathbf{p}$ is
outside}\\\ b&\mbox{if $\mathbf{p}$ is inside}\\\ \end{array}\right.$

where $a$ and $b$ are real numbers with $a>b$ to be used to identify inside
and outside.

Note that this is a terrible case. The function is not only discontinuous, but
its boundary is so complicated that it is hard to assign a normal vector. To
compute the gradient in such cases, the first step in Hart’s method is to
create a set of unit vectors $\vec{\mathbf{n}}_{k}$ that are uniformly
distributed on the surface of a unit sphere of dimensions $n$, where
$k=0,1,\cdots,K-1$ (see Figure 1(b)). We then classify these vectors as
inside- and outside-vectors for a given position $\mathbf{p}_{0}$. The vector
$\vec{\mathbf{n}}_{k}$ is an inside vector if
$F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{k})=a$, otherwise it is an outside vector
(see Figure 1(c)). The normal vector is then calculated by adding all inside
vectors (see Figure 1(d)). Note that the main idea in this method is that the
tangent components cancel out each other, while the perpendicular directions
reinforce each other. If we choose $a=1$ and $b=0$, this can be formulated as
follows.

$\left.\nabla
F(\mathbf{p})\right|_{\mathbf{p}=\mathbf{p}_{0}}\approx\sum_{k=0}^{K-1}F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{k})\vec{\mathbf{n}}_{k}.$

Note that the values of $a$ and $b$ can be anything (as far as the $a>b$
condition is satisfied) to obtain a somewhat acceptable normal vector. Since
this uses all vectors, the result is slightly more robust. A good choice for
the values of $a$ and $b$ is $a=1$ and $b=-1$.

(a) A black-box function that defines an interface.

(b) Uniformly distributed 2D unit vectors.

(c) Classifications of unit vectors as inside and outside.

(d) An estimated normal by adding outside vectors.

Figure 1: Text.

###  2.3 Antialising with Multisample Method

One problem with the multisampling method is that it can suffer from aliasing.
To solve it, we need to randomize the unit vectors while keeping their
distribution intact Cook et al. (1984); Cook (1986). The simple solution is to
rotate all vectors at random. For each new position to compute the gradient
vector, we can create a random rotation matrix $R$ in $SO(3)$, then update all
unit vectors as follows:

$\vec{\mathbf{n}}_{k}\leftarrow R\vec{\mathbf{n}}_{k}.$

Note that in this case, we cannot randomize vectors since that can change the
distribution. The random rotation of all vectors creates a new randomization
while keeping the original distribution of vectors. This idea is nothing more
than a reformulation of the original stochastic sampling Cook et al. (1984);
Cook (1986) that is adapted for a uniformly distributed set of unit vectors.

##  3 Theoretical Framework

Now, we are ready to extend the multi-sample method to be used beyond to
estimate the gradients of Fractal-type black-box functions. First, we need to
define the derivative in a given direction.

###  3.1 Derivative in a Given Direction

Numerical derivatives in 1D can be easily extended to compute the derivative
of any high-dimensional function in any given direction. Let
$\vec{\mathbf{n}}$ denote a unit vector, and $\mathbf{p}_{0}$ denote a
position in a $n$-dimensional space. The $F(\mathbf{p})$ becomes a one-
dimensional function along the direction of a ray $\vec{n}$ emanating from
$\mathbf{p}_{0}$. The equation for this ray can be written simply using an
auxiliary parameter $t$ as

$\mathbf{p}=\mathbf{p}_{0}+\vec{\mathbf{n}}t$

and a one-dimensional function along this ray is given in terms of the
parameter $t$ as

$F(\mathbf{p}_{0},\vec{\mathbf{n}}(t))=F(\mathbf{p}_{0}+\vec{\mathbf{n}}t)$

whose derivative along this ray can be simply estimated by the classical
method as follows

$\left.\frac{dF(\mathbf{p}+\vec{\mathbf{n}}t)}{dt}\right|_{\mathbf{p}=\mathbf{p}_{0}}\approx\frac{F(\mathbf{p}_{0}+\vec{\mathbf{n}}h)-F(\mathbf{p}_{0}-\vec{\mathbf{n}}h)}{2h}$
(2)

We can also use Squire and Trapp’s method by replacing $t$ with $ih$ Al-Mohy
and Higham (2010); Squire and Trapp (1998) as follows

$\left.\frac{dF(\mathbf{p}+\vec{\mathbf{n}}t)}{dt}\right|_{\mathbf{p}=\mathbf{p}_{0}}\approx\frac{\Im(F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{k}ih))}{h}.$
(3)

###  3.2 Derivative in a Given Line

Note that a line passing through a point consists of two rays defined by two
opposing vectors. The derivative has to be the same in both directions once we
consider the vector directions. In other words, we expect the following.

$F(\mathbf{p}_{0},\vec{\mathbf{n}}(t))\vec{\mathbf{n}}=F(\mathbf{p}_{0},\vec{\mathbf{-n}}(t))\vec{\mathbf{-n}}$

where $\vec{\mathbf{n}}+(\vec{\mathbf{-n})}=0$. This is a useful property,
since if we estimate the derivative in both directions, their average can
provide us with a better estimate by averaging errors if we assume the errors
to be random. Note that Equation 2 is written as an average of the two
derivatives. However, the negative creates a problem in estimation. Equation 3
does not have this problem, and the result can be improved by computing the
average as follows:

$\left.\frac{dF(\mathbf{p}+\vec{\mathbf{n}}t)}{dt}\right|_{\mathbf{p}=\mathbf{p}_{0}}\approx\frac{\Im(F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{k}ih))+\Im(F(\mathbf{p}_{0}+(\vec{\mathbf{-n}}_{k}ih)))}{2h}.$
(4)

##  4 Multi-Axis Gradient Estimation

Now, we are ready to develop the first method, which we call the multi-axis
gradient. In this case, we simply compute the gradient for multiple axes and
find the average. The following subsection presents the multi-axis gradient
method by estimating the gradient in multiple sets of orthonormal vectors.

###  4.1 Gradient with Orthonormal Vectors

Let $\vec{\mathbf{n}}_{j}$, where $j=0,1,\cdots,n-1$ denote any set of
orthonormal vectors in the $n$ dimensional space, $\Re^{n}$. The orthonormal
means that all these unit vectors are perpendicular to each other.

$\nabla
F(\mathbf{p})|_{\mathbf{p}=\mathbf{p}_{0}}=\sum_{j=0}^{n-1}\vec{\mathbf{n}}_{j}\left.\frac{dF(\mathbf{p}+\vec{\mathbf{n}_{j}}t)}{dt}\right|_{\mathbf{p}=\mathbf{p}_{0}}.$

Our intuition observation is that this equation must be correct regardless of
how orthogonal vectors are chosen. Now, consider a 2D case to test this
intuitive observation.

$\nabla
F(\mathbf{p})\approx\vec{\mathbf{n}_{0}}\frac{dF(\mathbf{p}+\vec{\mathbf{n}_{0}}t)}{dt}+\vec{\mathbf{n}_{1}}\frac{dF(\mathbf{p}+\vec{\mathbf{n}_{1}}t)}{dt}.$

where two perpendicular unit vectors can be given either
$\vec{\mathbf{n}}_{0}=(\cos\theta,\sin\theta)$ and
$\vec{\mathbf{n}}_{1}=(-\sin\theta,\cos\theta)$, or
$\vec{\mathbf{n}}_{0}=(\cos\theta,\sin\theta)$ and
$\vec{\mathbf{n}}_{1}=(\sin\theta,-\cos\theta)$. Now, let us check if this
gradient gives us the actual formula of the gradient:

$\displaystyle\nabla F(\mathbf{p})$ $\displaystyle=$
$\displaystyle\frac{dF(\mathbf{p}+\vec{\mathbf{n_{0}}}t)}{dt}\vec{\mathbf{n_{0}}}+\frac{dF(\mathbf{p}+\vec{\mathbf{n_{1}}}t)}{dt}\vec{\mathbf{n_{1}}}$
$\displaystyle=$
$\displaystyle\frac{dF(x+t\cos\theta,y+t\sin\theta)}{dt}(\cos\theta,\sin\theta)+\frac{dF(x-t\sin\theta,y+t\cos\theta)}{dt}(-\sin\theta,\cos\theta)$
$\displaystyle=$ $\displaystyle\left(\frac{\partial F(x,y)}{\partial
x}\cos\theta+\frac{\partial F(x,y)}{\partial
y}\sin\theta\right)(\cos\theta,\sin\theta)+\left(\frac{-\partial
F(x,y)}{\partial x}\sin\theta+\frac{\partial F(x,y)}{\partial
y}\cos\theta\right)(-\sin\theta,\cos\theta)$ $\displaystyle=$
$\displaystyle\frac{\partial F(x,y)}{\partial x}(1,0)+\frac{\partial
F(x,y)}{\partial y}(0,1)$

Note that this is not complete proof since it just holds for 2D cases. On the
other hand, it provides additional support to our intuitive observation.

Let us assume that our intuitive observation supported for 2D cases is correct
for higher dimensions. This means that we can compute multiple estimations of
the gradient using different sets of orthonormal vectors. For calculation, we
can use either the classical method or the imaginary method such as the one
given below:

$\nabla
F(\mathbf{p})|_{\mathbf{p}=\mathbf{p}_{0}}\approx\frac{1}{h}\sum_{j=0}^{n-1}\Im\bigl{(}F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{j}ih)\bigr{)}\vec{\mathbf{n}}_{j}.$

Each of these calculations provides an estimate based on the precision of the
computation method. Now, let us also assume that in these estimates the errors
are symmetric. we can calculate the gradient as the average of multiple
estimates. If the symmetry condition holds, this average should give us a
better estimate of the gradient. An additional property is that, even when
derivatives do not exist, it gives us a good estimate of the gradient vector.

It is important to note that any orthonormal set of vectors in n-dimension can
be obtained using a transformation matrix $M\in O(n)$. Here, $O(n)$ is called
the orthogonal group in n-dimension, and it is the group of all distance-
preserving transformations of an n-dimensional Euclidean that preserve a fixed
point. Note that these operations include mirror operations along with
rotation operations. Now, let $M_{k}\in O(n)$ where $k=0,1,\cdots,K-1$ and let
$\vec{\mathbf{n}}_{j}$ denote a set of initial orthonormal vectors where
$j=0,1,\cdots,n-1$, then we can compute the average gradient starting from the
following formula:

$\nabla
F(\mathbf{p})|_{\mathbf{p}=\mathbf{p}_{0}}=\frac{1}{K}\sum_{k=0}^{K-1}\sum_{j=0}^{n-1}M_{k}\vec{\mathbf{n}}_{j}\left.\frac{dF(\mathbf{p}+M_{k}\vec{\mathbf{n}_{j}}t)}{dt}\right|_{\mathbf{p}=\mathbf{p}_{0}}.$

To ensure that the errors are symmetric, we need to get $M$’s by equally
sampling $O(n)$. If the samples are skewed in one region, the average may not
necessarily be good. Higher dimensions require special care to obtain the
symmetric error.

We want to point out that the term

$\left.\frac{dF(\mathbf{p}+M_{k}\vec{\mathbf{n}_{j}}t)}{dt}\right|_{\mathbf{p}=\mathbf{p}_{0}}$

can be calculated by either the standard method or the imaginary method to
obtain an estimate. This gives us two ways to compute the gradient.

###  4.2 Demonstration of Multi-Axis Method for 2D

It is easy to demonstrate this work for 2D, as it is straightforward to sample
equally $O(2)$. Consider a set of unit vectors that are equally distributed in
a circle. These vectors can be created using the functions $\cos$ and $\sin$
as follows:

$\vec{\mathbf{n}}_{k}=\left(\cos\frac{\pi k}{K},\sin\frac{\pi
k}{K}\right),\qquad k=0,1,2,\ldots,K-1.$

Here, $M$’s are implicit, which are given as rotation matrices that are
equally spaced $S(2)$. If $K$ is divisible by four, then this set of vectors
provides us with the $K/2$ number of the axis. For any given point
$\mathbf{p}=(x,y)$ in any given direction $\vec{\mathbf{n}}_{k}$, we
approximate the derivative along this ray as:

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle x+t\cos\frac{\pi k}{K}$
$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle y+t\sin\frac{\pi k}{K}$

so that for any given 2D function $F(x,y)$, the derivative along this ray is
computed as follows by replacing $t$ with $ih$:

$\nabla F(\mathbf{p})\approx\frac{\Im(F(x_{0}+\cos\frac{\pi
k}{K}ih,y_{0}+\sin\frac{\pi k}{K}ih))}{h}.$

To see that this method symbolically produces the same result as the standard
gradient, now consider the simple circle equation $F(x,y)=x^{2}+y^{2}$, we can
ignore the constant term since it will not appear in the imaginary term, and
therefore:

$\frac{1}{h}\Im\left(\left(x+\cos\frac{\pi
k}{K}ih\right)^{2}+\left(y+\sin\frac{\pi
k}{K}ih\right)^{2}\right)=2x\cos\frac{\pi k}{K}+2y\sin\frac{\pi k}{K}$ (6)

whose gradient is computed as:

$\nabla F(\mathbf{p})\approx\sum_{k=0}^{K-1}\left(2x\cos\frac{\pi
k}{K}+2y\sin\frac{\pi k}{K}\right)\left(\cos\frac{\pi k}{K},\sin\frac{\pi
k}{K}\right).$ (7)

For $K=4$, the formula consists of the addition of the four vectors as
follows:

$\displaystyle\nabla F(\mathbf{p})$ $\displaystyle\approx$
$\displaystyle\sum_{k=0}^{3}\left(2x\cos\frac{\pi k}{4}+2y\sin\frac{\pi
k}{4}\right)\left(\cos\frac{\pi k}{4},\sin\frac{\pi k}{2}\right)$ (8)
$\displaystyle=$
$\displaystyle\left(4x\right)\left(1,0\right)+\left(4y\right)\left(0,1\right)$
$\displaystyle=$ $\displaystyle(4x,4y)$

In this case, there are two sets of axes, and when we divide this result by
the number of axes, the computation gives us an actual gradient $(2x,2y)$.
This is the simplest version of the multi-axis. For any value of $K$ divisible
by $4$, we obtain a multi-axis method. In this case, the result becomes

$\nabla F(\mathbf{p})\approx\left(2x\sum_{k=0}^{K-1}\cos^{2}\frac{\pi
k}{K}+2y\sum_{k=0}^{K-1}\sin\frac{\pi k}{K}\cos\frac{\pi
k}{K},2y\sum_{k=0}^{K-1}\sin^{2}\frac{\pi
k}{K}+2x\sum_{k=0}^{K-1}\sin\frac{\pi k}{K}\cos\frac{\pi k}{K}\right)$ (9)

In Equation 9, if $K$ is even

$\displaystyle\sum_{k=0}^{K-1}\cos^{2}\frac{\pi k}{K}$ $\displaystyle=$
$\displaystyle\frac{K}{2}$ (10) $\displaystyle\sum_{k=0}^{K-1}\sin\frac{\pi
k}{K}\cos\frac{\pi k}{K}$ $\displaystyle=$ $\displaystyle 0$ (11)

As a result, the Equation 9 is calculated as

$\nabla F(\mathbf{p})\approx(Kx,Ky)$ (12)

Note that to obtain the actual gradient, we need to take the average by
dividing the resulting vector with $K/2$, which corresponds to averaging of
the computation. Therefore, the multi-axis method can potentially give better
numerical results since the average removes the errors.

We also point out that in this example $K$ must be divisible by four since we
use uniformly spaced vectors. For single-axis computation, we will need to use
one of the two following axes:

$\displaystyle\vec{\mathbf{n}}_{0}=\left(\cos\theta,\sin\theta\right)\;\;\;\vec{\mathbf{n}}_{1}=\left(\cos(\theta+\pi/2),\sin(\theta+\pi/2)\right)$
$\displaystyle\vec{\mathbf{n}}_{0}=\left(\cos\theta,\sin\theta\right)\;\;\;\vec{\mathbf{n}}_{1}=\
\left(\cos(\theta-\pi/2),\sin(\theta-\pi/2)\right)$ (13)

For higher dimensions, $K$ must be at least a multiple of $n$ to have a multi-
axis. In higher dimensions, the normalization term will be $K/n$ instead of
$K/2$. It is also important to note that this symbolic computation is possible
only for the $2D$ case. For higher dimensions, we do not always get such nice
closed formulas. Another important point is that even in the 2D case, there
will be round-off errors that can accumulate in numerical computations.

Another tricky problem is the number of samples. Note that to estimate
gradients with the single-axis method, we only need the $n$ number of
orthonormal vectors in $n$ dimensions. The standard finite-difference formula
effectively uses its mirror. In practice, we use the $2N$ number of vectors.
In other words, the axis consists of $2n$ vectors. This is the reason in a 2D
space that we need $4$ vectors instead of $2$. If we populate the space with
these axes, we need to have $K=2nM$ number of vectors, where any $M$ is any
positive integer. This is why we need $K$ divisible by $4$ in 2D. In higher
dimensions, populating the space with a uniformly distributed set of $M$
number of axes may not always be straightforward.

Now, assume that $K$ is not divisible by $2n$ or even $n$. We also assume that
the vectors may not include orthonormal sets. It is still possible to
calculate a sum and divide the sum by $K/n$ to obtain an average that gives us
an estimation of the gradient. The main problem is that we will not have the
theoretical guarantee that there will be a $K/n$ number of axes, and each of
them gives us a good estimate of the actual gradient. On the other hand, from
the multi-sample method discussed in Section 2.2 we know that it could still
provide an approximate estimate. This observation gives us the next method.

###  4.3 Multi-Vector Gradient Estimation

Now, we create a set of unit vectors $\vec{\mathbf{n}}_{k}$ that are uniformly
distributed on the surface of a unit sphere of dimensions $n$, where
$k=0,1,\cdots,K-1$. This is an interesting problem in itself. We assume that
we can do this in the general case. This problem is very simple in 2D as we
discussed in the previous section. The problem is still reasonably easy in 3D
regular and semi-regular polyhedra. However, for higher-dimensional cases, we
need to use regular and semi-regular polytops Coxeter (1973, 1985), which can
be more difficult to conceptualize. Now, let us assume that it is possible to
obtain uniform distribution in any dimension.

For every unit vector $\vec{\mathbf{n}}_{k}$, we compute a derivative of
$F(\mathbf{p})$ along the direction of a ray $\vec{\mathbf{n}}_{k}$ emanating
from $\mathbf{p}_{0}$. The equation for this ray can be written simply using
an auxiliary parameter $t$ as

$\mathbf{p}=\mathbf{p}_{0}+\vec{\mathbf{n}}_{k}t$

whose derivative along this ray we estimate by replacing $t$ with $ih$ Al-Mohy
and Higham (2010); Squire and Trapp (1998)

$\frac{\Im(F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{k}ih))}{h}.$

The gradient of $F(\mathbf{p})$ is then approximated as the sum of all
derivative estimates multiplied with their corresponding unit vectors as
follows:

$\nabla
F(\mathbf{p})|_{\mathbf{p}=\mathbf{p}_{0}}=\frac{1}{h}\sum_{k=0}^{K-1}\Im\bigl{(}F(\mathbf{p}_{0}+\vec{\mathbf{n}}_{k}ih)\bigr{)}\vec{\mathbf{n}}_{k}.$

Note that this formula needs to be resized by dividing $K/n$ to compute the
average over all the sample points.

##  5 Uniformly Distributed Vector Creation

In the 2D case, the creation of equally spaced unit vectors can be done using
the vertex positions of the regular polygons that are circumscribed inside a
unit circle whose center is located in the origin such as a regular triangle,
square, or regular pentagon. As discussed in Section 4.2, for the multi-axis
method we need regular polygons with $4n$ sides such as square, regular
octagon, regular dodecagon, and regular hexadecagon. Note that for the multi-
axis method, it is also possible to use only half of the circle by eliminating
the inverse vectors. In other words, if $\vec{\mathbf{n}}_{k}$ is used,
$-\vec{\mathbf{n}}_{k}$ will not be used. This removes half of the vectors,
and using only even numbers of vectors, we can estimate the gradient with the
multi-axis method. Note that polygons with an even number of corners but not
divisible by $4$, such as regular hexagons and regular decagons, cannot still
be used for the multi-axis method since they do not include perpendicular
vectors. Of course, these and regular polygons with odd numbers of sides can
be used with the multi-vector method. This discussion provides a glimpse of
how the problem can be complicated in higher dimensions.

In the 3D case, the creation of equally spaced unit vectors can be done using
the vertex positions of the regular or semiregular polyhedra that are
circumscribed inside of a unit sphere whose center is located in origin. Some
examples are shown in Figure 2. However, most of them are useful only for
multi-vector gradient estimation. For example, the cube cannot be used for
multi-axis gradient estimation since its number of vertices is not divisible
by $6$. The octahedron is the simplest that can be used for multi-axis
gradient estimation (see Figure 2(b)). The icosahedron with its $12$ vertices
can also be used for multi-axis gradient estimation (see Figure 2(d)). Note
that the number of vertices of the icosahedron is not only divisible by $6$,
but its vertex positions consist of two distinct sets of perpendicular
directions. Similarly, the truncated octahedron shown in Figure 2(f) consists
of four distinct sets of perpendicular directions and can be used for multi-
axis gradient estimation. Any of these polyhedra. i.e. the octahedron, the
icosahedron, and the truncated octahedron can be used to create new vectors by
randomly rotating them. Since the initial distribution is uniform, the random
rotations create vectors that can cover the surface of the sphere in a
reasonably uniform way. We can approach this problem in higher dimensions in a
similar way.

(a) Tetrahedron.

(b) Octahedron.

(c) Cube.

(d) Icosahedron.

(e) Dodecahedron.

(f) Truncated Octahedron.

(g) Soccer Ball.

(h) Geodesic Dome (Not Semi-Regular).

Figure 2: Examples of regular and semiregular polyhedra. To use them for
multisampling, their centers must be the origin, and every vertex position
must be given by a unit vector. That unit vector provides the direction. Note
that not all of them can still be used for 3D multisampling. In these cases,
the octahedron is the shape that corresponds to the classical method by
providing three orthogonal directions. Note that the four vertices of the
tetrahedron provide the same distribution as the vertices of the cube.
Therefore, the tetrahedron is practically useless. The icosahedron provides 10
direction vectors, and the dodecahedron provides 6 direction vectors. For
semiregular polyhedra, all vertices must be congruent for multisampling. The
two examples here, the truncated octahedron and the soccer ball, are such
examples. If we remove the same distribution criteria, we can use a wide
variety of shapes such as the Geodesic Domes.

###  5.1 Discussion

In this section, we provide a brief discussion on difficulty in comparison of
the numerical computational methods with ground truths. Let the shapes $S$ be
defined through implicit functions $F(\mathbf{p})$’s as

$S=\\{\mathbf{p}|F(\mathbf{p})\leq 0\\}.$

One tricky problem for implicit functions, there exist infinitely many
functions that can describe the same shape. Concerning gradients, we are
interested in accurately describing tangent hyper-planes on our shape
boundaries whose unit vectors are

$\vec{\mathbf{n}}=\frac{\nabla F(\mathbf{p})}{\|\nabla F(\mathbf{p})\|}.$

In other words, we are interested only in the direction of these normal
vectors, independent of how the functions are formed. A simple example that
demonstrates this idea is a circle.

Notice that a circle of radius $r$ may be described using any of the following
functions:

$\displaystyle F_{0}(\mathbf{p})$ $\displaystyle=$ $\displaystyle
x^{2}+y^{2}-r^{2}$ $\displaystyle F_{1}(\mathbf{p})$ $\displaystyle=$
$\displaystyle\frac{x}{r}^{2}+\frac{y}{r}^{2}-1$ $\displaystyle
F_{2}(\mathbf{p})$ $\displaystyle=$ $\displaystyle\sqrt{x^{2}+y^{2}}-r$
$\displaystyle F_{3}(\mathbf{p})$ $\displaystyle=$
$\displaystyle(x^{2}+y^{2})^{a}-r^{2}a$

where $a$ is an arbitrary positive real number. The lengths of the gradient
vectors for these functions are vastly different, as can be easily seen in the
following list:

$\displaystyle\nabla F_{0}(\mathbf{p})$ $\displaystyle=$ $\displaystyle
2(x,y)$ $\displaystyle\nabla F_{1}(\mathbf{p})$ $\displaystyle=$
$\displaystyle\frac{2(x,y)}{r^{2}}$ $\displaystyle\nabla F_{2}(\mathbf{p})$
$\displaystyle=$ $\displaystyle\frac{(x,y)}{\sqrt{x^{2}+y^{2}}}$
$\displaystyle\nabla F_{3}(\mathbf{p})$ $\displaystyle=$ $\displaystyle
2a(x,y)(x^{2}+y^{2})^{a-1}$

On the other hand, it is also easy to check that their unit gradients are the
same on the boundary; they being precisely $\nabla F_{2}(\mathbf{p})$.

Therefore, it is possible to claim that it is important to check the accuracy
against the unit gradient instead of against its length. Note that if we agree
with this argument, there is no need to divide $K/n$. We can simply normalize
the gradient. One problem is that this normalization will also include some
errors. Let us assume that we can ignore this and we want to compare numerical
results with ground truth. In 2D, it can be easy to visualize the precision of
the unit numerical gradient $\vec{\mathbf{n}}_{n}$ with that of the ground
truth $\vec{\mathbf{n}}_{gt}$ using $\cos\theta$, where $\theta$ is simply the
angle between the two vectors, which can be computed using the dot product as
follows:

$\cos\theta=\vec{\mathbf{n}}_{n}\cdot\vec{\mathbf{n}}_{gt}.$

A nice property of this error measure is that it is bounded between $[-1,1]$,
where $1$ means no error, and $-1$ means the maximum possible. This error can
be visually represented through gray-level colors by mapping $[-1,1]$ to
$[0,1]$ using the following transformation:

$\epsilon=\frac{\cos\theta+1}{2}$

Based upon these observations, we initially conjecture that it can be possible
to compare the numerical gradient computations with the classical numerical
gradient computations and with the ground truth. One simply creates an
ancillary image of some image whose gradients are sought. The points
comprising this ancillary image are dot products that describe an error,
thereby producing a visual representation of the error over its spatial
domain. Images can be constructed using classical methods to estimate
numerical gradients. If a method turns out to be a significantly brighter
image, then we can visually claim that that method is superior for the problem
considered.

Unfortunately, in these cases, we quickly realized that even ground truths are
approximate since they should be calculated from symbolic results. The
additional normalization term makes the ground-truth computation not that
reliable. We also realized that low-dynamic range images are not suitable for
comparison. All methods provide relatively bright results that can be
difficult to visually compare results in a large domain. Therefore, we provide
a qualitative assessment of these methods.

  1. 1.

Single-Axis/Finite-Difference: This approach potentially suffers from
inaccuracies that come from a low number of sampling directions due to the use
of only a single set of orthogonal directions, and from the roundoff error
caused by the difference in the finite-difference formulæthat is used to
estimate the derivative.

  2. 2.

Single-Axis/Imaginary: The imaginary step derivative no longer succumbs to the
roundoff error, but the results still suffer from inaccuracies due to the low
number of sampling directions, because only a single orthogonal direction is
used.

  3. 3.

Multi-Axis/Finite-Difference: The integral method eliminates the problem of
low numbers of sampling directions, but the roundoff error from the standard
derivative approximation can still be significant.

  4. 4.

Multi-Axis/Imaginary: The imaginary step derivative coupled with integral
sampling significantly improves the results.

  5. 5.

Single-Axis/Imaginary: The imaginary step derivative no longer succumbs to the
roundoff error, but the results still suffer from inaccuracies due to the low
number of sampling directions, because only a single orthogonal direction is
used.

  6. 6.

Multi-Vector methods: These are not expected to be as good as multi-axis
methods for estimating gradients. However, they have the significant advantage
that we can use any semi-regular polyhedra or semi-regular polytops with
congruent vertices to obtain uniformly distributed vectors on the surfaces of
unit spheres or unit hyperspheres.

##  6 Conclusion and Future Work

This is a preliminary study that shows the potential of the derivative
estimation of Squire and Trapp for estimating the gradient in higher
dimensions. The extension of these gradient estimation methods to higher
dimensions is theoretically straightforward. Unfortunately, it is hard to
compute the quality of the results, since ground truths for the gradients may
not necessarily be reliable. We are still able to make some observations on
the nature of these operations. We believe that there is strong merit in
combining Squire and Trapp with the multi-axis extension of multi-sample
methodHart et al. (1989). For example, we suspect that these methods can be
used to obtain more accurate and fast solutions for optimization problems.
Therefore, we believe that these methods could be better tested to solve
optimization problems by demonstrating faster and more accurate results.

## References

  * Squire and Trapp [1998] William Squire and George Trapp.  Using complex variables to estimate derivatives of real functions.  _SIAM Review_ , 40(1):110–112, 1998.
  * Hart et al. [1989] John C Hart, Daniel J Sandin, and Louis H Kauffman.  Ray tracing deterministic 3-d fractals.  _ACM SIGGRAPH Computer Graphics_ , 23(3):289–296, 1989.
  * Greenwell [1998] Raymond N Greenwell.  Using complex variables to estimate derivatives of real functions.  _The College Mathematics Journal_ , 29(5):447, 1998.
  * Cook et al. [1984] Robert L Cook, Thomas Porter, and Loren Carpenter.  Distributed ray tracing.  In _Proceedings of the 11th annual conference on Computer graphics and interactive techniques_ , pages 137–145, 1984.
  * Cook [1986] Robert L Cook.  Stochastic sampling in computer graphics.  _ACM Transactions on Graphics (TOG)_ , 5(1):51–72, 1986.
  * Al-Mohy and Higham [2010] A. H. Al-Mohy and N. J. Higham.  The complex step approximation to the Fréchet derivative of a matrix function.  _Numerical Algorithms_ , 53:133–148, 2010.
  * Coxeter [1973] Harold Scott Macdonald Coxeter.  _Regular polytopes_.  Courier Corporation, 1973.
  * Coxeter [1985] HSM Coxeter.  Regular and semi-regular polytopes. ii.  _Mathematische Zeitschrift_ , 188(4):559–591, 1985.


# Informative Planning in the Presence of Outliers

Weizhe Chen and Lantao Liu The authors are with the Luddy School of
Informatics, Computing, and Engineering, Indiana University, Bloomington, IN,
47408, USA. {chenweiz<EMAIL_ADDRESS>
###### Abstract

Informative planning seeks a sequence of actions that guide the robot to
collect the most informative data to map a large environment or learn a
dynamical system. Existing work in informative planning mainly focus on
proposing new planners, and applying them to various robotic applications such
as environmental monitoring, autonomous exploration, and system
identification. The informative planners optimize an objective given by a
probabilistic model, e.g. Gaussian process regression. In practice, the model
can be easily affected by the ubiquitous sensing outliers, resulting in a
misleading objective. A straightforward solution is to filter out the outliers
in the sensing data stream using an off-the-shelf outlier detector. However,
informative samples are also scarce by definition, so they might be falsely
filtered out. In this paper, we propose a method to enable the robot to re-
visit the locations where outliers were sampled besides optimizing the
informative planning objective. By doing so, the robot can collect more
samples in the vicinity of outliers and update the outlier detector to reduce
the number of false alarms. This is achieved by designing a new objective on
top of a Pareto variant of Monte Carlo tree search. We demonstrate that the
proposed framework achieves better performance than simply applying an outlier
detector.

##  I INTRODUCTION

Oftentimes intelligent robots need to sense the environment with onboard
sensors and use collected data to understand their surroundings before
performing other subsequent tasks. _Informative planning_ [1] seeks a sequence
of actions that allow the robot to obtain the most informative data, i.e. the
data that contributes most to learning the environment model, while minimizing
the cost of traveling and sampling. Usually the level of informativeness for
sampling data from a certain location can be quantified by the reduction of
predictive uncertainty in a probabilistic model. For example, when dealing
with spatial phenomena, a common choice for the probabilistic model is
_Gaussian process regression_ (GPR), and we typically use _mutual information_
or _variance/entropy reduction_ to quantify the information contained in the
data [2, 3].

|
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
Figure 1: (a) Sensing data of Chlorophyll level over time. Most of the orange
and red points are considered as outliers since the variation of Chlorophyll
cannot take place this fast [4]. (b) Ground-truth environmental model and the
dataset. Green points are normal data points and red points are outliers. Red
and blue colors on the surface represent high and low values, respectively.
(c) Prediction of the GPR trained on the normal data. The upper and lower
black grids indicate the upper and lower confidence bounds. (d) Two-down view
of the predictive standard deviation of the GPR in (c). The black tree shows
the searching process of an informative planner while green line is the
resulting sampling trajectory. (e-f) are the same as (c-d) but the GPR is
contaminated by outliers. The small amount of outliers not only affect the
prediction and uncertainty estimate of the GPR but also the informative
planning result.

Many informative planners have been proposed during the last two decades [5,
6, 7, 8, 9, 10]. The majority of the literature assumes the sensing data is
the ground-truth value corrupted by adding some Gaussian noise. In practice,
however, sensing data typically contains outliers due to various reasons. For
instance, ranging sensors such as sonars can contain failed specular
reflection measurements, leading to returned false values of maximum range
even there are obstacles blocking the beams; When there are multiple
exteroceptive sensors in use, cross talks among them also cause “random
measurements” [11]; For outdoor environments, unexpected fast transient
obstacles such as insects, fish, or large dust particles in air/water can also
cause unavoidable outliers. In addition, sensor imperfect conditions such as
worn probes, occasional hardware malfunction, and many other unfathomable
reasons can also lead to outliers. Gaussian model can be easily skewed by
outliers, which makes the “optimal” informative sampling decisions computed
from such a model become sub-optimal. Fig. 1 illustrates a scenario with
outliers that mislead the subsequent modeling/learning and planning
performances.

A straightforward idea to mitigate the influence of outliers is to filter them
out from the sensing data stream using an off-the-shelf outlier detector.
However, this leads to high false-alarm rate because informative samples and
outliers are similar in nature: they both change the model significantly and
are observed less frequently.  Otherwise, the informative samples would not
bring much additional information. To reduce the high false-alarm rate, we
propose to optimize two objectives simultaneously. In addition to the original
information-seeking objective, the other objective encourages re-visiting the
locations where outliers were sampled. In this way, the robot can collect more
samples in the region where outliers are detected frequently and update the
outlier detector to reduce the number of false alarms. We implement this idea
by designing a new objective in the recently proposed Pareto Monte Carlo tree
search [12]. Combined with a state-of-the-art outlier detector, we demonstrate
that the proposed framework achieves better performance than applying only the
outlier detector.

Although dealing with outliers is essential in practical applications of
informative planning algorithms, to the best of our knowledge, it has not been
well discussed in the literature. The contribution of this paper is to point
out the interesting relationship between informative samples and outliers:
outliers can also be viewed as “extremely informative samples”. With this
discovery, we develop an effective approach to significantly reduce the number
of false alarms brought by this relationship.

##  II Related Work

Informative planning has attracted growing interest especially in the scenario
of autonomous environmental monitoring [13]. Informative planning extends the
optimal sensor placement problem [14] by considering the traversal costs and
motion constraints of the mobile robots. Starting from the seminal work [15],
various recursive greedy algorithms have been proposed [1, 16, 17]. These
methods are based on the _submodularity_ of the objective function and provide
performance guarantee. The submodularity requirement has been relaxed to
monotonicity [18]. Dynamic programming (DP) based methods do not require the
objective function to have these special properties. In [2], a sequence of
informative waypoints are selected via DP by assuming the underlying map to be
in a rectangular and sliced shape. This framework was extended to arbitrary
continuous space by connecting the informative waypoints via traveling
salesman problem solver [19]. Integrated with online sparse Gaussian
processes, this framework is lifted to online planning [3]. To achieve
efficient informative planning, Hollinger et. al. [8] extends sampling-based
motion planner, e.g., rapidly-exploring random tree, to robotic information
gathering algorithms. Recently, this framework has been developed to online
variants [5, 20]. Monte-Carlo tree search (MCTS) based methods are
conceptually similar to sampling-based informative planners and have recently
garnered great attention [21, 7, 22, 12, 23]. Instead of randomly growing the
search tree, the MCTS expands tree nodes in a best-first search manner [24].
Sine Gaussian processes have been the de facto standard for modeling spatio-
temporal phenomena in many environmental monitoring applications, Bayesian
optimization becomes a natural choice for informative planning [9, 25, 26].
New frameworks for informative planning are constantly emerging with
attractive mechanisms such as evolutionary methods [6, 27, 28] and imitation
learning [10].

Many of the aforementioned works employ GPR as the probabilistic model.
Compared to the large body of work in informative planners, investigations on
the probabilistic model in the informative planning framework are relatively
sparse. Methods using online sparse GPs [29] have been proposed to tackle the
computational bottleneck. Mixture of GPs [30, 31] has also been used to
capture environmental non-stationarity.

Different from all existing work, we investigate the effects of sensing
outliers on model learning and informative planning. To the best of our
knowledge, this is the first time that outliers in informative planning has
been discussed.

##  III Problem Formulation

Let $\mathcal{T}\triangleq\\{1,2,\dots,T\\}$ be the set of decision epochs. At
time $t\in\mathcal{T}$, a robot with a fully observed state
$\mathbf{s}_{t}\in\mathcal{S}$ collects scalar measurement $y_{t}$ via its
sensor and takes an action $a_{t}$, arriving at the next state
$\mathbf{s}_{t+1}$ with transition probability
$p(\mathbf{s}_{t+1}|\mathbf{s}_{t},a_{t})$. In this work, we assume the
transition to be known and deterministic and focus on inferring the noise-free
environment state $f_{\text{env}}$ from its noisy observation collected so far
$\mathbf{y}=[y_{1},\dots,y_{N}]^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}$.
One example is to estimate an elevation map from range sensor measurements. In
this case, the elevation map is a function $f_{\text{env}}(\mathbf{x})$ of the
input spatial locations $\mathbf{x}\in\mathbb{R}^{D}$. The robot maintains an
internal model $f$ of the environment learned from the noisy observations
$\mathbf{y}$. The model is updated after taking new observations into account.

Ideally, informative planning should use error reduction as the reward signal
for the planner. Let $f_{e}(\bm{a})$ be a function that quantifies the
modeling error reduction after incorporating the data obtained by taking an
action sqeuence $\bm{a}$, and $f_{c}(\bm{a})$ be a function that returns the
cost of executing $\bm{a}$. Informative planning requires solving the
following optimization problem given an available budget $B$,

$\bm{a}^{\star}=\operatorname*{arg\,max}_{\bm{a}\in\bm{\mathcal{A}}}f_{e}(\bm{a}),~{}s.t.~{}f_{c}(\bm{a})\leq{B},$
(1)

where $\mathcal{A}$ is the space of all possible action sequences. However,
the ground-truth environment state $f_{\text{env}}$ is unknown, hindering the
computation of error reduction. Informative planning bypasses this problem by
utilizing a surrogate reward function $f_{i}(\mathbf{a})$ given by a
probabilistic model instead of $f_{e}(\mathbf{a})$. Typically, $f_{i}(\bm{a})$
quantifies the informativeness of the data obtained by taking actions $\bm{a}$
by measuring the reduction in predictive uncertainty of the probabilistic
model.

##  IV Method

Figure 2: Illustration of each module of the system.

Fig. 2 illustrates the system with different modules and their relationship at
a high level. The outlier detector takes raw sensing data, filters out the
outliers, and outputs an outlier matrix indicating the density of outliers at
each location (Sec. IV-A). The clean data is fed into the GPR model which
produces the predictive mean and standard deviation at each location (Sec.
IV-B). The standard deviation matrix and the outlier matrix together form the
multi-objective reward for the Pareto Monte Carlo tree search (Sec. IV-C).
Finally, the robot collects samples along the optimal trajectory given by the
planner.

###  IV-A Copula-Based Outlier Detection

COPula-based Outlier Detector (COPOD) is one of the top performing outlier
detector that is efficient and tunning-free [32]. It is available in the open-
source PyOD library [33] and can be readily used in our problem. Firstly,
COPOD fits empirical left tail cumulative distribution functions (CDFs)
$\mathrm{F}^{l}_{1}(y)\dots,\mathrm{F}^{l}_{O}(y)$ using

$\mathrm{F}^{l}_{o}(y)=\mathbb{P}((-\infty,y])=\frac{1}{N}\sum_{n=1}^{N}\mathbb{I}(y_{n,o}\leq{y}),$
(2)

and empirical right tail CDFs
$\mathrm{F}^{r}_{1}(y)\dots,\mathrm{F}^{r}_{O}(y)$ by replacing $y$ with $-y$,
where $O$ is the observation dimension. The skewness vector
$\mathbf{b}=[b_{1},\dots,b_{O}]$ is also computed using the standard
estimation formula

$b_{o}=\frac{\frac{1}{N}\sum_{n=1}^{N}{(y_{n,o}-\overline{y_{o}})}^{3}}{\sqrt{\frac{1}{N-1}\sum_{n=1}^{N}{(y_{n,o}-\overline{y_{o}})}^{2}}^{3}},$
(3)

where the overline $\overline{y_{o}}$ denotes the mean. Then, we apply
$F^{l}_{o}(y)$ and $F^{r}_{o}(y)$ on each sample to get the empirical copula
observations $l_{n,o}$and $r_{n,o}$, respectively. The skewness corrected
empirical copula observations are given by

$s_{n,o}=\begin{cases}l_{n,o}&\text{ if }b_{o}<0,\\\ r_{n,o}&\text{
otherwise}.\end{cases}$ (4)

Lastly, for each dimension, we compute the probability of observing a point at
least as extreme as each $y_{n,o}$. The outlier score is the maximum negative
log probability given by the left-tail, right-tail, and sknewness-corrected
empirical copula. Intuitively, outliers should occur less frequently, which
means that the tail probabilities should be small, or equivalently, its
negative log probabilities should be large. We refer the reader to [32] for
more details.

###  IV-B Gaussian Process Regression

A Gaussian process (GP) is a collection of random variables, any finite number
of which have a joint Gaussian distribution [34]. We place a Gaussian process
prior over the function

$f(\mathbf{x})\sim\mathcal{GP}(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime}))$
(5)

which is specified by a mean function $m(\mathbf{x})$ and a covariance
function $k(\mathbf{x},\mathbf{x}^{\prime})$ (a.k.a. kernel). Popular choices
of kernel functions include Matérn and Gaussian kernels. In this paper, we use
a zero mean and the anisotropic Gaussian kernel

$k(\mathbf{x},\mathbf{x}^{\prime})=\alpha^{2}\exp{\left(-\frac{1}{2}\sum_{d=1}^{D}\frac{(x_{d}-x_{d}^{\prime})^{2}}{\ell_{d}^{2}}\right)}.$
(6)

The amplitude $\alpha^{2}$ controls the variance of the random functions and
$\ell_{d}$ defines the lengthscale, which informally can be thought of as the
“radius of a neighborhood” in the $d$-th dimension. GP Regression (GPR)
augments the model by assuming observations are corrupted by an additive
Gaussian white noise

$p(y|\mathbf{x})=\mathcal{N}(y|f(\mathbf{x}),\sigma^{2}).$ (7)

We will collect the parameters into
$\bm{\theta}\triangleq\\{\alpha,\ell_{d},\sigma\\}$ which are termed as
hyperparameters.

According to the definition of GP and the likelihood in Eq. (7), the joint
distribution of the observations $\mathbf{y}$ and a set of latent function
values $\mathbf{f}_{\star}$ at arbitrary test input locations
$\mathbf{X}_{\star}$ is a multivariate Gaussian distribution

$\begin{bmatrix}\mathbf{y}\\\
\mathbf{f}_{\star}\end{bmatrix}\sim\mathcal{N}\left(\mathbf{0},\begin{bmatrix}\mathbf{K}_{y}&\mathbf{K}_{\star}\\\
\mathbf{K}_{\star}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}&\mathbf{K}_{\star\star}\end{bmatrix}\right),$
(8)

where $\mathbf{K}_{y}=\mathbf{K}_{\mathbf{x}}+\sigma^{2}\mathbf{I}$,
$\mathbf{K}_{\mathbf{x}}$ is the covariance matrix given by the covariance
function evaluated at each pair of observations, $\mathbf{K}_{\star}$ is the
covariance matrix between observations and the test function values
$\mathbf{f}_{\star}$, and $\mathbf{K}_{\star\star}$ is the covariance matrix
of the test function values. The predictive distribution for GPR is a
conditional Gaussian distribution

$\displaystyle p(\mathbf{f}_{\star}\rvert\mathbf{y})$
$\displaystyle=\mathcal{N}(\mathbf{f}_{\star}\rvert{\bm{\mu}},\bm{\Sigma}),\text{
where}$ (9) $\displaystyle\bm{\mu}$
$\displaystyle=\mathbf{K}_{\star}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\mathbf{K}_{y}^{-1}\mathbf{y},$
(10) $\displaystyle\mathbf{V}$
$\displaystyle=\mathbf{K}_{\star\star}-\mathbf{K}_{\star}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\mathbf{K}_{y}^{-1}\mathbf{K}_{\star}.$
(11)

Learning in Gaussian process regression refers to determining appropriate
values of the hyperparameters, and the common approach is to maximize the log
marginal likelihood of the observations

$\ln{p(\mathbf{y}|\bm{\theta})}=-\frac{1}{2}\left(\mathbf{y}^{{\mathchoice{\raisebox{0.0pt}{$\displaystyle\intercal$}}{\raisebox{0.0pt}{$\textstyle\intercal$}}{\raisebox{0.0pt}{$\scriptstyle\intercal$}}{\raisebox{0.0pt}{$\scriptscriptstyle\intercal$}}}}\mathbf{K}_{y}^{-1}\mathbf{y}+\ln{\det(\mathbf{K}_{y})}+N\ln(2\pi)\right),$
(12)

where $\det(\cdot)$ denotes the matrix determinant.

Figure 3: (a) A matrix storing the number outliers occur at each location. (b)
A smoothed version of (a) using a Gaussian kernel.

###  IV-C Pareto Monte Carlo Tree Search

We treat the informative planning problem as a sequential decision making
problem and approximate the solution using Monte Carlo tree search (MCTS)
[24]. A node in the tree contains a state/pose and some statistics. The MCTS
iterates in four steps: _node selection, expansion, simulation_ , and _back-
propagation_. We now describe the reward function and each step of the MCTS in
detail.

Reward Functions: The reward of an action should be the amount of the
“informativeness” of the collected data by taking this action. A natural
choice for measuring informativeness is the mutual information between the
visited locations and the remainder of the space [14]. However, calculating
entropy involves the determinant computation of the predictive covariance
matrix of the GPR in Eq. (11), which is computationally costly. A cheaper
approximation is the reduction in the trace of the predictive covariance
matrix, namely, variance reduction. Observing that the predictive uncertainty
at a location becomes very small after sampling at that location, we further
simplify the maximization of variance reduction to simply maximizing the sum
of predictive standard deviation along the sampling trajectory. To encourage
re-visiting the locations where outliers were detected, another reward is the
number of outlier occurance. For better planning efficiency, we pre-compute
the predictive standard deviation and the number of outliers on a _query grid_
representing all discretized sampling locations before the tree search
process. When evaluating the reward function values, we simply access the
corresponding matrix elements. To better guide the robot, we also filter the
outlier occurance matrix via a Gaussian kernel. Fig. 3 shows the outlier
occurance matrix and its smoothed version.

####  IV-C1 Selection

As a best-first search algorithm, MCTS selects nodes that are expected to have
a high reward but at the same time try other nodes sufficiently to avoid
greedy choices. We trade-off _exploration_ and _exploitation_ using the tree
policy. A well known tree policy is the upper confidence bound (UCB) [35, 36].
We recursively select child nodes with the highest UCB value until a node with
unexpanded children is found. For node $j$, the UCB value adapted to our
problem is defined by

$\operatorname{UCB}_{j}=\overline{f_{i}}_{j}+C\sqrt{\frac{2\log
N_{p}}{N_{j}}},$ (13)

where $N_{p}$ is the number of times that the parent of node $j$ has been
selected and $N_{j}$ is the number of times that node $j$ has been selected.
The average reward $\overline{f_{i}}_{j}$ encourages the selection of nodes
that currently look the most promising and thus stimulates exploitation. The
value of the second term is large when $N_{i}$ is small, which encourages
exploring node $i$ if it is not selected enough times. Constant $C$ balances
between the exploration and the exploitation which should be set to a similar
scale as the rewards.

For vector reward given by the multiple reward functions, we use the Pareto
variant of UCB proposed in [12]:

$\operatorname{PUCB}_{j}=\overline{\mathbf{f}_{i}}_{j}+C\sqrt{\frac{4\log
N_{p}+\ln{D_{r}}}{2*N_{j}}},$ (14)

where $\overline{\mathbf{f}_{i}}_{j}$ is the average reward vector and $D_{r}$
is the dimensionality of this vector. When selecting the best child node, we
first compute the Pareto optimal set from the Pareto UCB vectors of all child
nodes. The vectors in the Pareto optimal set cannot be improved for any
objective without hurting other objectives. We then randomly select a child
node from this set because the Pareto optimal solutions are considered equally
optimal if no preference information is given. We refer the reader to [12] for
further details of how Pareto MCTS balances multiple objectives. When the UCB
tree policy is applied to a MCTS, the algorithm is commonly referred to as
UCT. Similarly, we use PUCT to represent Pareto UCB applied to MCTS.

####  IV-C2 Expansion

We randomly select an action from the selected node and delete this action
from the list of available actions. We then take this action and calculate the
corresponding reward. A new node is created based on the new state and reward.
We append this new node to the children list of the selected node.

####  IV-C3 Simulation

The goal of this step is to estimate the expected reward of the action of the
newly expanded node by executing a default/rollout policy. Here we define the
rollout policy to be “moving forward” for some steps. This policy is designed
to estimate the expected reward along the direction of executing an action.

####  IV-C4 Back-propagation

We add the average reward received in simulation to each node along the
selection path and increase the number of visits of these nodes.

These four steps are repeated until the computational budget for the robot has
reached. We recursively select the child node with the highest number of
visits to get the final informative action sequence. Fig. 4 demonstrates the
planning results of PUCT and UCT. We can observe that PUCT can balance the two
rewards.

|
---|---
(a) PUCT | (b) UCT
Figure 4: Illustrations of the searching results of PUCT and UCT. In the
occupancy map of (a), black and white colors represent occupied space and free
space, respectively. In a reward map, darker color indicates higher reward
value. The goal is to search a trajectory that maximizes the cumulative
reward. The green tree shows the searching process which outputs the best
trajectory colored in red and best action shown in blue. The first reward map
tells the PUCT to search towards the upper right while the second reward map
encourages going up. The reward map of the UCT is the same as the first reward
map in (a). The PUCT balances the two rewards and gives different trajectory
at the end.

###  IV-D Overall Framework

Collect the initial training data

Initialize and optimize gpr with the initial data

Initialize the outlier matrix outlier to be zeros

while  _has sampling budget_ do

mean,std=gpr.predict(x_test)

trajectory=PUCT(pose,std,outlier)

Follow trajectory and sample new data

Train the outlier detector on all collected data

Filter out the outliers in the new data

gpr appends the new data and optimize for several iterations

end while

Algorithm 1 Main Method |  |
---|---|---
(a) Ground-truth environment | (b) PUCT-COPOD final predictive mean | (c) UCT-COPOD final predictive mean
Figure 5: Ground-truth environment and the final predictive means of PUCT-
COPOD and UCT-COPOD.

Alg. 1 summarized main steps of the proposed framework. First, we collect the
_pilot_ data following a path that does not depend on the informative planning
framework. For example, we will use Bezier curve in Fig. 6 (a). The pilot data
is used for an initial optimization of the GPR hyperparameters through Eq.
(12) and computing the statistics for pre/post-processing of the data. In a
decision epoch, the robot searches for an optimal trajectory using PUCT based
on the predictive standard deviation given by the GPR and the outlier matrix
given by the outlier detector. The robot then collects samples along the
trajectory. Finally, we detect outliers in this batch of new samples, get rid
of them, feed the clean data to the GPR, and optimize the hyperparameters for
several iterations.

|
---|---
(a) | (b)
Figure 6: (a) The Bezier curve to generate _a priori_ sampling path for
collecting initial training data for the initial optimization of the GPR and
computing data pre/post-processing statistics. (b) Final predictive standard
deviation of UCT-NONE.

##  V Experimental Results

We would like to answer the following questions.

  * •

How do outliers affect an informative planning system if they are not filtered
out? To answer this question, we use the standard GPR and UCT without any
outlier detector as the worst-case baseline (UCT-NONE).

  * •

What is the ideal performance if we have a perfect outlier detector? We use
the ground-truth outlier labels as the prediction of a perfect outlier
detector. We refer to this best-case competitor as UCT-BEST.

  * •

What if we simply apply an off-the-shelf detector to deal with the outliers?
In this method, we use COPOD to get rid of outliers (UCT-COPOD).

  * •

Is the proposed method able to mitigate the negative effect of outliers? We
denote our proposed method as PUCT-COPOD.

|
---|---
(a) $\rho=0.05$ | (b) $\rho=0.1$
|
(c) $\rho=0.15$ | (d) $\rho=0.1$
Figure 7: Root mean squared errors of the proposed PUCT-COPOD and other three
baselines. UCT-BEST filters outliers with the ground-truth labels, so it
serves as a best-case baseline; UCT-NONE, which is the worst case baseline,
does not filter outliers at all. UCT-COPOD directly applies the outlier
detector to the sensing data stream. Parameter $\rho$ controls the number of
outliers: higher value indicates more outliers (MO) while lower value means
less outliers (LO). In (a)-(c), the outlier detector only filters out the
outliers in the newly acquired batch without changing the historical data. In
(d), the outlier detector re-examine _all_ the historical data at each
decision epoch. (a) Performance difference of all the methods is insignificant
when the number of outliers is very small. (b) PUCT-COPOD and UCT-COPOD are
both in between the best-case baseline and the worst-case baseline after
collecting $1000$ samples while PUCT-COPOD performs better. At the early
stage, UCT-COPOD is even worse than UCT-NONE because it filters out some
informative samples. (c) When the number of outliers increases, UCT-COPOD
loses advantage over UCT-NONE. The performance of PUCT-COPOD is also affected
but its still better than UCT-COPOD and UCT-NONE. (d) Both PUCT-COPOD and UCT-
COPOD are ineffective if the detector filters outliers in all historical data.

Considering that the second question requires the ground-truth outlier labels
and to better understand the results, we inject some “spike” outliers to the
sensing data which are similar to those in Fig. 1 (a). The robot collects a
batch of data along the planned trajectory (red arrows in Fig. 4) in each
decision epoch. We randomly select $\rho$ portion of the batch for outlier
injection. Specifically, we first sample a random amplitude uniformly from the
range $[1,2]$ with a random sign. Then, the amplitude is multiplied with the
data range of the batch computed by the $0.05$ and $0.95$ quantiles.

We use the elevation map of the Mount Saint Helens, shown in Fig. 5(a), as the
environment. The planar robot follows the Dubins car dynamics:
$\begin{bmatrix}\dot{x}_{1},\dot{x}_{2},\dot{\theta}\end{bmatrix}=\begin{bmatrix}v\cos{\theta},v\sin{\theta},u\end{bmatrix},u\in\mathcal{U},$
where $\mathbf{x}=[x_{1},x_{2}]^{\intercal}\in\mathbb{R}^{2}$ and
$\theta\in[0,2\pi)$ are the position and orientation of the vehicle, $u$ is
the control input that represents the robot’s steering angle. We set the
linear velocity $v$ to be a constant and only control the steering angle $u$.
The possible steering angles $[-0.15,0.15]$ (in radians) are discretized into
$5$ choices.

The GPR is initialized using $100$ training samples generated by a Bezier
curve shown in Fig. 6 (a) and optimize the hyperparameters for $500$
iterations. This _a priori_ Bezier sampling path can be applied to any
environment for collecting the _pilot_ data for the initial GPR optimization
and computing data pre/post-processing statistics. The path circles around the
whole environment to help us get a sense of the workspace extent and a rough
estimate of the observation range. Using the statistics computed from the
initial training data, all inputs $\mathbf{X}$ will be scaled to the range of
$[-1,1]$ and all observations $\mathbf{y}$ will be standardized to have zero
mean and unit variance. We use $500$ tree search iterations and $5$ rollout
iterations. The sampling budget is set to $2000$ samples.

|
---|---
(a) UCT-COPOD | (b) PUCT-COPOD
Figure 8: (a) and (b) shows the absolute error maps and training samples after
cleaning the outliers of UCT-COPOD and PUCT-COPOD, respectively. Many
informative samples at the high-variability region are falsely filtered out in
(a), leading to some high-error “holes”. As a comparison, (b) shows that the
proposed method can mitigate this effect to some extent.

Fig. 7 shows the root mean squared error (RMSE) of the four methods versus the
number of training samples. When the number of outliers is very small
($\rho=0.05$ in Fig. 7 (a)), the difference among all the methods is
negligible. After increasing $\rho$ to $0.1$, as expected, the error of UCT-
BEST drops drastically and reaches the lowest error. After the decrease at the
beginning, the error of UCT-None is actually increasing. The inclusion of
outliers not only makes the prediction deteriorated but also messes up the
standard deviation (see Fig. 6 (b)), which further affects the downstream
informative planning. The result of UCT-COPOD is interesting — although the
final error is lower than that of UCT-NONE, it is the worst at the beginning.
The reason for this can be seen from Fig. 8 (a). Many informative samples at
the high-variability region are falsely filtered out, leading to some high-
error “holes”. As a comparison, Fig. 8 (b) shows that the proposed method can
mitigate this effect to some extent. This brings better RMSE curve to PUCT-
COPOD in Fig. 7 (b). When there are more outliers ($\rho=0.15$), UCT-COPOD
loses the advantage over UCT-NONE. The performance of PUCT-COPOD is also
affected but its still better than UCT-COPOD and UCT-NONE. Finally, we also
study the case where the detector filters outliers in all the historical data
instead of the newly acquired batch. Both PUCT-COPOD and UCT-COPOD are
ineffective in this case. One way to understand the result is that, when the
robot has not sampled the volcanic area, the detector can easily detect
extremely high or low values because the detector is trained on the historical
data sampled from the flatland. However, after collecting samples from the
volcanic area and retaining the detector, it is difficult for the detector to
distinguish the aforementioned extreme values and the normal elevation
measurements around the volcanic area during re-examination. In addition to
the quantitative difference in the error reduction, the final predictions of
the proposed method and the baseline are also qualitatively different. Fig. 5
shows that the prediction of PUCT-COPOD is closer to the ground-truth
environment than that of UCT-COPOD.

##  VI CONCLUSION

We present a framework to enable the robot to re-visit the locations where
outliers were sampled besides optimizing the conventional informative planning
objective. This is very different from existing informative planning
approaches which mainly focus on planners that optimize over uncertainty or
confidence based informativeness, which can be easily deteriorated by various
sensing outliers, resulting in a misleading objective and sub-optimal sampling
behaviors. We propose an approach to filter out the outliers from the sensing
data stream using an off-the-shelf outlier detector. By designing a new
planning objective with a Pareto variant of Monte Carlo tree search, our new
framework allows the robot to collect more samples in the vicinity of the
outliers and update the outlier detector to reduce the number of false alarms.
Results show that the proposed framework performs much better than only
applying the outlier detector.

## References

  * [1] A. Singh, A. Krause, C. Guestrin, W. J. Kaiser, and M. A. Batalin, “Efficient planning of informative paths for multiple robots,” in _IJCAI_ , vol. 7, 2007, pp. 2204–2211.
  * [2] N. Cao, K. H. Low, and J. M. Dolan, “Multi-robot informative path planning for active sensing of environmental phenomena: A tale of two algorithms,” in _Proceedings of the 2013 international conference on Autonomous agents and multi-agent systems_. International Foundation for Autonomous Agents and Multiagent Systems, 2013, pp. 7–14.
  * [3] K.-C. Ma, L. Liu, H. K. Heidarsson, and G. S. Sukhatme, “Data-driven learning and planning for environmental sampling,” _Journal of Field Robotics_ , vol. 35, no. 5, pp. 643–661, 2018.
  * [4] B. Espinola, “Statistical analysis of water data from an online exo2 monitoring sonde: Modelling chlorophyll-a using an exo2 sonde,” Ph.D. dissertation, Metropolia University of Applied Sciences, Helsinki, Finland, 2017\.
  * [5] L. M. Schmid, M. Pantic, R. Khanna, L. Ott, R. Siegwart, and J. Nieto, “An efficient sampling-based method for online informative path planning in unknown environments,” _IEEE Robotics and Automation Letters_ , 2020.
  * [6] M. Popović, G. Hitz, J. Nieto, I. Sa, R. Siegwart, and E. Galceran, “Online informative path planning for active classification using uavs,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 5753–5758.
  * [7] G. Best, O. M. Cliff, T. Patten, R. R. Mettu, and R. Fitch, “Dec-mcts: Decentralized planning for multi-robot active perception,” _The International Journal of Robotics Research_ , vol. 38, no. 2-3, pp. 316–337, 2019\.
  * [8] G. A. Hollinger and G. S. Sukhatme, “Sampling-based robotic information gathering algorithms,” _The International Journal of Robotics Research_ , vol. 33, no. 9, pp. 1271–1287, 2014.
  * [9] R. Marchant and F. Ramos, “Bayesian optimisation for intelligent environmental monitoring,” in _2012 IEEE/RSJ international conference on intelligent robots and systems_. IEEE, 2012, pp. 2242–2249.
  * [10] S. Choudhury, M. Bhardwaj, S. Arora, A. Kapoor, G. Ranade, S. Scherer, and D. Dey, “Data-driven planning via imitation learning,” _The International Journal of Robotics Research_ , vol. 37, no. 13-14, pp. 1632–1672, 2018.
  * [11] S. Thrun, “Probabilistic robotics,” _Communications of the ACM_ , vol. 45, no. 3, pp. 52–57, 2002.
  * [12] W. Chen and L. Liu, “Pareto monte carlo tree search for multi-objective informative planning,” in _Proceedings of Robotics: Science and Systems_ , FreiburgimBreisgau, Germany, June 2019.
  * [13] M. Dunbabin and L. Marques, “Robots for environmental monitoring: Significant advancements and applications,” _IEEE Robotics & Automation Magazine_, vol. 19, no. 1, pp. 24–39, 2012.
  * [14] A. Krause, A. Singh, and C. Guestrin, “Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies,” _Journal of Machine Learning Research_ , vol. 9, no. Feb, pp. 235–284, 2008\.
  * [15] C. Chekuri and M. Pal, “A recursive greedy algorithm for walks in directed graphs,” in _46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05)_. IEEE, 2005, pp. 245–253.
  * [16] A. Meliou, A. Krause, C. Guestrin, and J. M. Hellerstein, “Nonmyopic informative path planning in spatio-temporal models,” in _AAAI_ , vol. 10, 2007, pp. 16–7.
  * [17] G. Hollinger, S. Singh, J. Djugash, and A. Kehagias, “Efficient multi-robot search for a moving target,” _The International Journal of Robotics Research_ , vol. 28, no. 2, pp. 201–219, 2009.
  * [18] J. Binney, A. Krause, and G. S. Sukhatme, “Optimizing waypoints for monitoring spatiotemporal phenomena,” _The International Journal of Robotics Research_ , vol. 32, no. 8, pp. 873–888, 2013.
  * [19] K.-C. Ma, L. Liu, and G. S. Sukhatme, “An information-driven and disturbance-aware planning method for long-term ocean monitoring,” in _Intelligent Robots and Systems (IROS), 2016 IEEE/RSJ International Conference on_. IEEE, 2016, pp. 2102–2108.
  * [20] M. Ghaffari Jadidi, J. Valls Miro, and G. Dissanayake, “Sampling-based incremental information gathering with applications to robotic exploration and environmental monitoring,” _The International Journal of Robotics Research_ , vol. 38, no. 6, pp. 658–685, 2019.
  * [21] A. Arora, P. M. Furlong, R. Fitch, S. Sukkarieh, and T. Fong, “Multi-modal active perception for information gathering in science missions,” _Autonomous Robots_ , pp. 1–27, 2018.
  * [22] P. Morere, R. Marchant, and F. Ramos, “Sequential bayesian optimization as a pomdp for environment monitoring with uavs,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 6381–6388.
  * [23] G. Flaspohler, V. Preston, A. P. Michel, Y. Girdhar, and N. Roy, “Information-guided robotic maximum seek-and-sample in partially observable continuous environments,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3782–3789, 2019.
  * [24] C. B. Browne, E. Powley, D. Whitehouse, S. M. Lucas, P. I. Cowling, P. Rohlfshagen, S. Tavener, D. Perez, S. Samothrakis, and S. Colton, “A survey of monte carlo tree search methods,” _IEEE Transactions on Computational Intelligence and AI in games_ , vol. 4, no. 1, pp. 1–43, 2012.
  * [25] S. Bai, J. Wang, F. Chen, and B. Englot, “Information-theoretic exploration with bayesian optimization,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2016, pp. 1816–1822.
  * [26] C. K. Ling, K. H. Low, and P. Jaillet, “Gaussian process planning with lipschitz continuous reward functions: Towards unifying bayesian optimization, active learning, and beyond,” in _Thirtieth AAAI Conference on Artificial Intelligence_ , 2016.
  * [27] M. Popović, T. Vidal-Calleja, G. Hitz, J. J. Chung, I. Sa, R. Siegwart, and J. Nieto, “An informative path planning framework for uav-based terrain monitoring,” _Autonomous Robots_ , vol. 44, no. 6, pp. 889–911, 2020.
  * [28] G. Hitz, E. Galceran, M.-È. Garneau, F. Pomerleau, and R. Siegwart, “Adaptive continuous-space informative path planning for online environmental monitoring,” _Journal of Field Robotics_ , vol. 34, no. 8, pp. 1427–1449, 2017.
  * [29] K.-C. Ma, L. Liu, and G. S. Sukhatme, “Informative planning and online learning with sparse gaussian processes,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 4292–4298.
  * [30] W. Luo and K. Sycara, “Adaptive sampling and online learning in multi-robot sensor coverage with mixture of gaussian processes,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 6359–6364.
  * [31] R. Ouyang, K. H. Low, J. Chen, and P. Jaillet, “Multi-robot active sensing of non-stationary gaussian process-based environmental phenomena,” in _Proceedings of the 2014 international conference on Autonomous agents and multi-agent systems_ , 2014, pp. 573–580.
  * [32] Z. Li, Y. Zhao, N. Botta, C. Ionescu, and X. Hu, “Copod: copula-based outlier detection,” in _2020 IEEE International Conference on Data Mining (ICDM)_. IEEE, 2020, pp. 1118–1123.
  * [33] Y. Zhao, Z. Nasrullah, and Z. Li, “Pyod: A python toolbox for scalable outlier detection,” _Journal of Machine Learning Research_ , vol. 20, pp. 1–7, 2019\.
  * [34] C. E. Rasmussen and C. K. I. Williams, _Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)_. The MIT Press, 2005.
  * [35] P. Auer, “Using confidence bounds for exploitation-exploration trade-offs,” _Journal of Machine Learning Research_ , vol. 3, no. Nov, pp. 397–422, 2002\.
  * [36] L. Kocsis and C. Szepesvári, “Bandit based monte-carlo planning,” in _European conference on machine learning_. Springer, 2006, pp. 282–293.