# An Open-Source Gazebo Plugin for GNSS Multipath Signal Emulation in Virtual
Urban Canyons

Kartik Anand Pant111Graduate Student, School of Aeronautics and Astronautics.
<EMAIL_ADDRESS>Purdue University, West Lafayette, IN, 47906  Zhanpeng
Yang222Graduate Student, School of Aeronautics and Astronautics.
<EMAIL_ADDRESS>Purdue University, West Lafayette, IN, 47906  James M
Goppert  Purdue University, West Lafayette, IN, 47906  333Lecturer, School of
Aeronautics and Astronautics<EMAIL_ADDRESS>and Inseok Hwang 444
Professor, School of Aeronautics and Astronautics<EMAIL_ADDRESS>and AIAA
Associate Fellow  Purdue University, West Lafayette, IN, 47906

###### Abstract

One of the major errors affecting GNSS signals in urban canyons is GNSS
multipath error. In this work, we develop a Gazebo plugin which utilizes a ray
tracing technique to account for multipath effects in a virtual urban canyon
environment using virtual satellites. This software plugin balances accuracy
and computational complexity to run the simulation in real-time for both
software-in-the-loop (SITL) and hardware-in-the-loop (HITL) testing. We also
construct a 3D virtual environment of Hong Kong and compare the results from
our plugin with the GNSS data in the publicly available Urban-Nav dataset, to
validate the efficacy of the proposed Gazebo Plugin. The plugin is openly
available to all the researchers in the robotics community.
https://github.com/kpant14/multipath_sim

##  1 Nomenclature

$m$ = | Multipath offset
---|---
$d$ = | Additional distance covered by the reflected GNSS signal before reflection
$l$ = | Additional distance covered by the reflected GNSS signal after reflection
$\theta_{i}$ = | Elevation angle of the $i^{th}$ satellite
$\theta_{0}$ = | Elevation angle threshold
$\alpha$ = | Angle subtended by the normal of the reflected ray with the reflected ray
$P_{i}$ = | Pseudorange of the $i^{th}$ satellite
$c_{0}$ = | speed of light
$\Delta t$ = | clock bias of the receiver with respect to the satellite
$e$ = | measurement error
$\sigma$ = | standard deviation of the measurement error
$\rho_{i}$ = | True range of the $i^{th}$ satellite
$(x_{i},y_{i},z_{i})$ = | Position coordinate of the $i^{th}$ satellite
$(x,y,z)$ = | Position coordinate of the receiver
$\sigma_{G}$ = | Geometric dilution-of-precision
$r^{LOS}$ = | Laser range for LOS ray
$r^{REF}$ = | Laser range for Mirror ray

##  2 Introduction

The Global Navigation Satellite System (GNSS) is being widely used for
numerous modern transportation and logistic applications, from school buses to
military convoys. A robust GNSS positioning fix is essential for safe and
secure navigation in a multitude of civil and military applications. With the
advent of urban air mobility (UAM), such accurate positioning capability
becomes paramount to the safe operation of aerial vehicles, whether unmanned
or passenger-carrying. The accuracy of GNSS receivers is affected by several
factors, for example, ionospheric errors, satellite clock errors, and
multipath errors. By using differential techniques, ionospheric effects can be
mitigated [1]. However, the effects of multipath errors are site and time
dependent and urban environments are particularly impacted by GNSS multipath.
There are two major challenges for GNSS operations in an urban canyon.
Firstly, the densely packed high-rise buildings significantly obstruct the
line-of-sight (LOS) of satellites. Secondly, the high-rise buildings also
provide ideal surfaces for the multipath of GNSS signals to occur, which
further reduces the accuracy of the position fix provided by the GNSS
constellation. Many studies have been conducted for multipath detection and
mitigation in the past. A 3D model of the environment is constructed in [2, 3]
to detect multipath errors by using the ray tracing method in an urban canyon
environment. However, both approaches require computationally expensive ray-
tracing algorithms to check every 3D surface in the virtual environment. A
similar approach using ray-tracing is developed in [4], however, it requires
precise knowledge of the satellite–reflector–antenna geometry, the reflector
material, and antenna characteristics. Statistical methods have been proposed
in [5, 6] to detect multipath using the double difference of the pseudo-ranges
and carrier phase measurements. The multipath errors have been modeled as
statistical distributions in [7, 8, 9]. In [10, 11], a simulated satellite
using a Zeppelin is utilized to compute the precise measurements of
statistical distributions of multipath delays. Although statistical methods
provide good estimates of the multipath properties, they require large-scale
raw GNSS data collection to fit the error profile of a specific region and
they do not accurately represent the position and time dependent offsets
inherent in multipath errors.

In this work, we present an open-source GNSS sensor emulation plugin for the
Gazebo simulator that can be easily integrated with popular flight control
firmware such as PX4 and Ardupilot. The proposed GNSS plugin utilizes a ray-
tracing algorithm on 3D building models. The building models in the Gazebo
simulator are generated using data obtained from Open Street Map (OSM), to
estimate the GNSS signal paths at each location in the simulated urban
environment. This framework can be utilized for:

  1. 1.

Software-in-the-Loop (SITL) testing of UASs in a virtual urban environment
constructed in Gazebo.

  2. 2.

Mixed Reality testing on real UAVs flying in the tracking volume of the motion
capture (MoCap) system while a virtual copy moves in the simulated urban
environment.

  3. 3.

Computing heat maps of GNSS multipath error for particular times and locations
in an urban environment.

The main contribution of our work is to provide an open-source GNSS emulation
plugin in the ROS-Gazebo environment that is accessible to the wider UAS
community. The rapid adaptation of Cyber-Physical Systems (CPS) in urban
environments in recent years have led to a strong emphasis on the safety and
security of autonomous systems. The proposed GNSS plugin can be used as an
effective tool to model and prototype the effects of GNSS multipath errors on
large-scale CPS systems. It can be further utilized to identify state and time
dependent vulnerabilities in the urban environment to synthesize stealthy GNSS
spoofing attacks and develop effective detection and mitigation strategies
against them.

The rest of the paper is organized as follows: Section 3 describes the
construction of the virtual 3D environment and the mathematical model used for
the computation of pseudo-ranges in the presence of multipath errors. Section
4 highlights the architecture followed for designing the proposed GNSS plugin
in Gazebo. A comparison of the performance of the proposed with a publicly
available GNSS dataset of Hong Kong city is presented in Section 5. Finally,
Section 6 summarizes the research and lays out future directions.

##  3 Background and Preliminaries

This section describes mathematical formulations used for modeling the
multipath effects on the GNSS signals and discusses the theory of operation in
estimating the receiver’s position. It also describes the construction of the
virtual urban 3D environment and virtual satellites in the Gazebo simulator.

The GNSS receivers periodically receive radio signals from various satellites
in a constellation and compute the ranges by measuring the signal delay. These
ranges are then used to estimate the position of the receiver. As the
satellite and the receiver have different clocks, the receiver also estimates
its own clock bias with respect to the satellite’s clock. In order to uniquely
calculate the coordinates of the receiver, the receiver requires a minimum of
four satellites to estimate its position and the clock bias [12]. In this
work, the GNSS signals are modeled as a straight line segment between the
transmitter at each satellite and the receiver on the Earth. In reality,
however, these signals do not follow a straight line path due to the
atmospheric effects. The straight-line assumption is valid as the effects of
atmospheric errors (ionospheric and tropospheric) on GNSS signals are well
modeled [13, 14]. The corrections for these errors are available on most
commercially available receivers through the wide-area augmentation system
(WAAS) [15].

This work focuses on the effects of multipath errors on GNSS localization
which are dependent on the geometry of the environment. Multipath errors are
difficult to model as they occur due to reflective surfaces such as tall
buildings, dense vegetation, the ground, etc. In the absence of LOS to a
specific satellite, the GNSS receiver may only receive reflected multipath
signals. The multipath signal travels a longer distance and thus creates a
synchronization bias with respect to LOS signals from other satellites,
resulting in errors in code and carrier pseudo-range measurements. Ultimately,
multipath errors lead to significant errors in positioning for the GNSS
receiver, especially in urban canyons.

To model the behavior of the multipath signals, first, a 3D virtual
environment of a part of the city is constructed in the Gazebo simulator,
followed by a ray tracing algorithm to compute the pseudo-range offset
resulting from the multipath effects. The modified pseudo-ranges are then
utilized to compute the position fix and the estimation accuracy, using the
least squares estimation algorithm. The details of each of these steps are
described below:

###  3.1 Virtual Urban Canyon and Virtual Satellites

To develop the 3D virtual environment in Gazebo, the building information for
the desired part of the city is obtained from the publicly available dataset
using Open Street Maps (OSM) [16]. The data from OSM is pre-processed to
extract the 3D meshes and geodetic coordinates of the buildings. The geodetic
positions are then converted into Universal Transverse Mercator (UTM)
coordinates to align with the coordinate system of Gazebo. The virtual
satellites are set up using the raw ephemeris data either retrieved from the
publicly available GNSS repositories [17] or acquiring it with an actual GNSS
receiver over a period of time. The raw ephemeris data contains the Keplerian
elements of the satellite’s orbit, for example, orbital inclination and
eccentricity. The orbital parameters are used to calculate the position of the
satellites in the Earth Centered Earth Fixed (ECEF) coordinate system. The
construction of the virtual 3D model of the city of Hong Kong is shown in Fig.
1.

Figure 1: Virtual Hong Kong Environment in Gazebo

###  3.2 Ray Tracing Algorithm

For each satellite, the ray-tracing algorithm simulates two rays originating
from the GNSS receiver, one connecting the satellite and receiver directly and
another in the complementary direction, i.e., the direction of the reflected
signal. The LOS signal is simulated as a straight line segment connecting the
GNSS receiver with the satellite. Similarly, the reflected signal is simulated
as a straight line segment connecting the GNSS receiver and the obstruction
(if there exists one). The lengths of these line segments are utilized to (1)
detect the presence of multipath error and (2) calculate the multipath offset,
i.e., additional distance covered by the reflected signal before reaching the
receiver.

Figure 2: Multipath reflection in a typical urban setting

In this work, only a single reflection of the ray from a building’s surface is
considered. This is because attenuation and degradation in the GNSS signal
strength due to multiple reflections causes less significant positioning error
as stated in [18]. In addition, it is assumed that the reflected surface is
normal to the ray. This is a simplification, but it captures many behaviors of
GPS multipath error without detailed information on the reflection geometry
and reflection surface composition. In practice, multipath reflections with
high incidence angles are often ignored due to outlier rejection. Therefore,
normal reflections are often the most detrimental to the performance of the
system, and with the assumption of a single reflection, normal incidence
angle, and clear LOS of the reflected signal we derive the 2D geometric
relationship shown in Fig. 2.

We use the following trigonometric identity to derive the additional path
length covered by the GNSS signals after reflection.

$\begin{gathered}m=d+l\\\ m=d(1+sin\alpha)\\\
m=d(1+sin(\pi/2-2\theta))\end{gathered}$ (1)

###  3.3 Measurement Noise

The pseudo ranges obtained from the GNSS constellations contain inherent
randomness which can be modeled as a combination of white noise, pink noise,
and random walk. Various noise models have been studied in the literature to
model GNSS measurement noise. [19, 20] highlight that correlated noise models
fit the GNSS measurement noise better than uncorrelated noise models. It also
highlighted that Ornstein Ulenback (OU) process can be utilized to model
correlated GNSS noise. An OU process has the following attractive properties:
mean-reverting, continuous and has few parameters; making it suitable for fast
embedded computation and sensor fusion applications.

The OU process is defined by the stochastic differential equation:

$dx_{t}=\theta_{ou}(\mu_{ou}-x_{t})+\sigma_{ou}dW_{t}$ (2)

where $\theta_{ou}$ and $\sigma_{ou}$ are parameters that can be tuned to
match the randomness corresponding to a specific model of GNSS receiver and
$\mu_{ou}$ is a parameter that can be adjusted to introduce drifting behavior
of a position fix when the receiver velocity is slow or losing LOS to too many
satellites. In this work, we use the OU process to model the measurement noise
by the GNSS receiver.

###  3.4 Localization

The pseudoranges calculated by the GNSS receiver describe the distance between
the satellites and the receiver as well as the errors from external sources,
for example, atmospheric delays and multipath interference. To imitate this
behavior in simulation, the true ranges are first computed using the true
positions of the virtual satellites and the true position of the receiver
using Gazebo/MoCap depending on the simulation/mixed reality test setup. For
each satellite $i$, the pseudorange $P_{i}$ is calculated as:

$P_{i}=[(x_{i}-x)^{2}+(y_{i}-y)^{2}+(z_{i}-z)^{2}]^{\frac{1}{2}}+m_{i}+c_{0}\Delta
t+e$ (3)

The right side of (3) is the sum of (i) mean-squared difference between the
$i^{th}$ satellite’s position $(x_{i},y_{i},z_{i})$ and the receiver’s
position $(x,y,z)$, (ii) multipath offset $m_{i}$ calculated by the ray
tracing algorithm, and (iii) additional delay terms $c_{0}\Delta t$ and
measurement noise, where $c_{0}$ represents the speed of light, $e$ accounts
for measurement noise and $\Delta t$ is the delay associated with the clock
offset between satellites and the receiver.

To calculate the position of the receiver using the pseudorange measurements,
a minimum of four independent pseudorange measurements are required. This is
because the system of equations contains four unknown variables $(x,y,z,\Delta
t)$. When the number of measurements is more than four, the system of
equations becomes over-determined. To decode the position of the receiver,
first (3) is linearized about some initial estimate or guess for the
receiver’s position (the linearization point) using Taylor’s theorem. Then,
the corrections to these initial estimates are recursively adjusted to obtain
the receiver’s actual position and clock offset [12]. The compact
representation of the set of the linearized equation in matrix form can be
described as:

$\Delta P=A\Delta x+e$ (4)

in which $\Delta P$ is a vector of difference between the corrected
pseudorange and modeled pseudorange of size $n$, where $n$ is the number of
pseudorange measurements, $\Delta x$ designates the four-element vector of
unknowns, the receiver’s position, and the clock offset, from the
linearization point. A is an $n\times 4$ matrix of the partial derivatives of
the pseudoranges $wrt$ the unknowns and is defined as:

$A=\begin{bmatrix}\frac{x_{1}-x}{\rho_{1}}&\frac{y_{1}-y}{\rho_{1}}&\frac{z_{1}-z}{\rho_{1}}&-1\\\
\vdots&\vdots&\vdots&\vdots\\\
\frac{x_{i}-x}{\rho_{i}}&\frac{y_{i}-y}{\rho_{i}}&\frac{z_{i}-z}{\rho_{i}}&-1\\\
\vdots&\vdots&\vdots&\vdots\\\
\frac{x_{n}-x}{\rho_{n}}&\frac{y_{n}-y}{\rho_{n}}&\frac{z_{n}-z}{\rho_{n}}&-1\\\
\end{bmatrix}$ (5)

where $\rho_{i}=[(x_{i}-x)^{2}+(y_{i}-y)^{2}+(z_{i}-z)^{2}]^{\frac{1}{2}}$ is
the actual distance between $i^{th}$ satellite and the receiver. The receiver
then solves the matrix equation (4) by minimizing the estimation error using
the least squares algorithm to get,

$\Delta x=(A^{T}A)^{-1}A^{T}\Delta P$ (6)

In general, when the initial estimates are incorrect, the solution of the non-
linear problem requires multiple iterations to converge. The least squares
algorithm is iterated until $\Delta x$ becomes close to a small value
$\epsilon$, which in practice is kept typically around $10^{-6}$.

###  3.5 Localization Accuracy

As the measurements obtained from the satellites are not perfect, due to noise
and other unmodeled errors, the position estimates given by (6) are also not
accurate. Assuming (a priori) standard deviation of modeling error in
pseudorange measurements, $\sigma$, the covariance matrix $Q$ of the expected
error in unknown parameters $(x,y,z,\Delta t)$ can be formulated as:

$Q=(A^{T}A)^{-1}\sigma^{2}\equiv D\sigma^{2}$ (7)

With a given standard deviation of measurement errors, $\sigma$, the quality
of least squares solutions can be computed by using the diagonal elements of
the $Q$ matrix described as:

$Q=\begin{bmatrix}{\sigma_{x}}^{2}&\sigma_{xy}&\sigma_{xz}&\sigma_{xt}\\\
\sigma_{yx}&{\sigma_{y}}^{2}&\sigma_{yz}&\sigma_{yt}\\\
\sigma_{zx}&\sigma_{zy}&{\sigma_{z}}^{2}&\sigma_{zt}\\\
\sigma_{tx}&\sigma_{ty}&\sigma_{tz}&{\sigma_{t}}^{2}\\\ \end{bmatrix}$ (8)

where $\sigma_{x},\sigma_{y},\sigma_{z}$ are the estimated rms errors in the
user position $(x,y,z)$ and $\sigma_{t}$ is the rms value of the error in the
user clock bias estimate. $\sigma_{xy},\sigma_{yz},\sigma_{zx}$ are covariance
in the position estimates in $x$, $y$ and $z$ directions, and
$\sigma_{xt},\sigma_{yt},\sigma_{zt}$ are the covariance of position estimates
$(x,y,z)$ with the clock bias $\Delta t$. Observe that the elements of matrix
$Q$ are dependent only on the receiver-satellite geometry. The overall quality
of the least-squares solution can be obtained by taking the square root of the
diagonal elements of the $Q$ matrix.

$\displaystyle\sigma_{G}$
$\displaystyle=\sqrt{{\sigma_{x}}^{2}+{\sigma_{y}}^{2}+{\sigma_{z}}^{2}+{\sigma_{t}}^{2}}$
(9)
$\displaystyle=\sqrt{{D_{11}}^{2}+{D_{22}}^{2}+{D_{33}}^{2}+{D_{44}}^{2}}\sigma$
(10)

The standard deviation, $\sigma_{G}$, is a scaled version of the standard
deviation of the measurement errors. As observed by authors in [21], the
scaling factor is typically greater than 1 which results in the standard
deviation of the least squares solution being greater than $\sigma$. This
inflation of variance of the estimated unknown variables due to receiver-
satellite geometry is known as Dilution of Precision (DOP) and $\sigma_{G}$ is
called geometric dilution of precision (GDOP). To understand this effect
better, let us consider two satellites transmitting to a ground receiver in
two different configurations as depicted in Fig. 3. The receiver’s true
position can be uniquely identified by intersecting the two range circles.
However, due to uncertainty in the measurements, the range circles are
inexact, resulting in an error in positioning. This error depends on the
geometry related to the receiver and the transmitters. When the transmitters
are located far apart shown in Fig. 3(a), there is a small region of
uncertainty where the receiver must be located, whereas, the region of
uncertainty is larger when the transmitters are closer to each other as shown
in Fig. 3(b).

(a) Small uncertainty when transmitters are far away

(b) Large uncertainty when transmitters are closer

Figure 3: Influence of satellite-receiver geometry on GNSS position accuracy
adapted from [22]

Instead of assessing the quality of the overall solution, typically the
specific components of DOPs, i.e., the three-dimensional position component
(PDOP), the horizontal component (HDOP), the vertical component (VDOP), and
the time component (TDOP) are evaluated. Each of these components is computed
as:

PDOP
$\displaystyle=\frac{\sqrt{{\sigma_{x}}^{2}+{\sigma_{y}}^{2}+{\sigma_{z}}^{2}}}{\sigma}$
(11) HDOP
$\displaystyle=\frac{\sqrt{{\sigma_{x}}^{2}+{\sigma_{y}}^{2}}}{\sigma}$ (12)
VDOP $\displaystyle=\frac{\sigma_{z}}{\sigma}$ (13) TDOP
$\displaystyle=\frac{\sigma_{t}}{\sigma}$ (14)

Note that $\text{PDOP}^{2}=\text{HDOP}^{2}+\text{VDOP}^{2}$ and
$\text{GDOP}^{2}=\text{PDOP}^{2}+\text{TDOP}^{2}$. The DOP values can be pre-
computed as they are only dependent on the receiver and satellite coordinates.
The higher the DOP values the poorer the accuracy of the GNSS position fix.
Although, there is no common agreement on the classification of the quality of
the GNSS solution based on DOP coefficient values. Based on the study
conducted in [21], Table. 2 highlights the rating of the resulting solution as
per the DOP values.

Table 2: Characterization based on DOP coefficient DOP value |  <1 |  1-2 |  2-5 |  5-10 |  10-20 |  >20
---|---|---|---|---|---|---
Rating |  Ideal |  Exellent |  Good |  Moderate |  Fair |  Poor

##  4 Architectural Design

Figure 4: GNSS multipath plugin flowchart

Before describing the architecture of the proposed plugin, it is important to
understand some useful terminologies related to ROS and Gazebo. This section
explains the basic tools and terminologies used for designing Gazebo plugins
and presents the software architecture of the proposed GNSS multipath plugin.

###  4.1 ROS and Gazebo

The Robot Operating System(ROS) is a framework for the development and
deployment of robotic systems [23]. It contains collections of drivers,
libraries, and tools for interprocess communication between various robot
hardware and software. ROS follows a computational graph structure in which
ROS processes are graph nodes. The nodes transfer information among each other
using edges called ROS topics. Each ROS topic broadcasts information in the
form of ROS messages. This inter-process communication allows easy development
of user libraries and executables.

Gazebo is a 3D simulator built for both indoor and outdoor robotic
applications [24]. Together with ROS, Gazebo provides an efficient physics
platform for high fidelity 3D simulation of various robotic systems. The
interface between ROS and Gazebo is controlled by a set of shared libraries
called gazebo_plugins.

###  4.2 Gazebo Plugins

A plugin in Gazebo is a shared library that can access all the functionalities
of the simulator from the physics engine to the simulated environment. The
plugins are flexible pieces of code that can be easily added and removed from
the existing simulation system. There are four different types of Gazebo
plugins to modify and control properties of various simulated entities:

  1. 1.

World Plugin: for altering physics settings.

  2. 2.

Sensor Plugin: for changing sensor parameters.

  3. 3.

Visual Plugin: for modifying visual properties of the 3D models.

  4. 4.

System Plugin: for system-level settings.

###  4.3 Simulation Description File (SDF)

Gazebo uses an XML-based configuration file for saving and loading the
parameters related to the simulated models. It captures the necessary
information about the scene, physics, lighting, and plugins of the Gazebo
simulator.

###  4.4 GNSS Multipath Plugin

The architecture of the GNSS plugin is similar to the different sensor plugins
already present in the Gazebo simulator. The sensor plugins can be attached to
any 3D object in the simulator using Gazebo’s SDF file. The plugin is
developed as a C++ class, which inherits from ray sensor plugin of the Gazebo.
In Gazebo, the ray sensor is the backbone for LiDAR and laser-based range
sensors. The line segments mentioned in the previous section are constructed
by orienting rays from the receiver in the direction of the satellites
(azimuth and elevation angles). Using the native C++ APIs of the ray sensor
plugin, the additional pseudorange offset is calculated by using algorithm 2.
Our major contribution is to modify and optimize the existing ray sensor
plugin APIs to allow real-time updating of the azimuth and elevation angles of
the satellites. The noise levels and satellite constellations can be modified
using the SDF, allowing the proposed GNSS plugin to simulate a wide range of
GNSS constellations and a variety of commercially available GNSS receivers.

Input: True Ranges $\rho_{i}$, satellites position
$(x_{sat},y_{sat},z_{sat})$, elevation angles $\theta_{i}$, multipath offset
$m_{i}$, measurement noise $e_{i}$

Output: Receiver’s position, ($x,y,z$)

for  _$i\in(1,n)$_ do

// Checking for low elevation angle

if  _$\theta_{i} >\theta_{0}$_ then

$m_{i}\leftarrow ComputeMultipathOffset(r_{i}^{LOS},r_{i}^{REF})$ // multipath
offset, $m_{i}$, is computed by the ray tracing algorithm

$P_{i}\leftarrow\rho_{i}+m_{i}+e_{i}$ // $P_{i}$ is the computed pseudorange,
$e_{i}$ is the measurement error

$vis\\_sat\leftarrow vis\\_sat+1$ // visible satellites are incremented

// Checking for number of visible satellites

if  _$vis\\_sat >4$_ then

$(x,y,z)\leftarrow RecursiveLeastSquares(P,(x_{0},y_{0},z_{0}))$ // receiver’s
position, $(x,y,z)$, calculated using RLS algorithm, $(x_{0},y_{0},z_{0})$ is
the initial guess of the receiver’s position

$DOP\\_coeff\leftarrow ComputeDOP((x_{sat},y_{sat},z_{sat}),(x,y,z))$ //
compute DOP coefficient using satellites and receiver’s estimated position

Algorithm 1 GNSS Localization

Input: Ranges $r_{i}^{LOS}$ and $r_{i}^{REF}$, $i\in(1,n)$ from the ray
sensor. where n is the number of satellites, $r_{i}^{LOS}$ is the laser range
for LOS ray and $r_{i}^{REF}$ is the laser range for mirror ray.

Output: Multipath offset, $m_{i}$

for  _$i\in(1,n)$_ do

// Checking for LOS obstruction

if  _$r_{i}^{LOS} <MAX\\_RANGE$_ then

// Checking whether the mirror ray exists

if  _$r_{i}^{REF} <MAX\\_RANGE$_ then

$m_{i}\leftarrow CalculateOffset(r_{i}^{REF})$// using trigonometric
identities to calculate pseudorange offset

// No GPS multipath

Algorithm 2 Multipath Offset Calculation

Algorithm. 1 describes the estimation of the receiver’s position and
calculation of the estimation accuracy. For calculating the input ranges,
$\rho$, the true position of the receiver is extracted from the
simulator/MoCap system based on the simulation/mixed reality experiment. The
plugin obtains the position of each satellite using the ephemeris data
provided. The true range for each satellite is computed using the extracted
position of the receiver and the satellite. Using Algorithm. 2, the plugin
computes the multipath offset, $m_{i}$, for each satellite. The low elevation
satellites are excluded from the estimation by checking against a threshold
angle, $\theta_{0}$, as those are most prone to atmospheric effects and
multipath errors [25]. Finally, the coordinates of the receiver are computed
only if the minimum number of visible satellites is more than 4. The recursive
least square estimation algorithm is used to estimate the receiver’s
coordinates. The various components of the DOP (PDOP, HDOP, VDOP, and TDOP)
are computed using the receiver and the satellite’s positions.

The multipath pseudorange offset is computed using the ray tracing algorithm
defined in Algorithm. 2. The ray sensor first calculates the ranges of line-
of-sight and reflected rays, $r_{i}^{LOS}$ and $r_{i}^{REF}$, to infer the
presence of multipath errors. Multipath errors are detected when the line-of-
sight ray is obstructed and there exists a obstructed reflected ray
(indicating the presence of a building or a reflecting surface). If the range
criterion is satisfied, the function calculateOffset computes the additional
pseudorange offsets using the derived trigonometric identities in (1). In the
absence of direct line of sight to the satellite (LOS) and lack of a building
for reflection, the satellite is declared as not visible (no signal received).

The plugin logs the computed coordinates of the receiver, components of the
DOP, visible satellites, noise levels, and the pseudorange offset of the
satellites. All this information is wrapped into a custom ROS message and
published as a ROS topic for downstream applications such as path planning and
navigation.

###  4.5 GNSS Multipath Plugin’s message and SDF

As mentioned previously, the plugin combines all the information into a ROS
message and publishes it using a ROS topic. The details of each entry of the
ROS message are tabulated in Table. 3.

Field |  Datatype |  Purpose
---|---|---
header |  std_msgs/Header |  Stores the timestamp and frame information
rec_pos |  sensor_msgs/NavSatFix |  Stores the receiver’s geodetic location
dop |  float32[] |  Stores the components of Dilution-of-Precision
range_offset |  float32[] |  Stores the multipath offset for each satellite[m]
sats_blocked |  int32[] |  Stores the PRN of blocked satellites
num_vis_sat |  int32 |  Stores the number of visible satellite
num_block_sat |  int32 |  Stores the number of blocked satellite
noise |  float32[] |  Stores the measurement noise [m]
Table 3: ROS message for GNSS multipath plugin

Below is the snippet of the SDF for adding the GNSS multipath plugin. To add
the plugin to an existing model these lines must be included. More details
about the SDF can be found at
https://github.com/kpant14/multipath_sim/tree/main/worlds

⬇

<plugin name="laserscan" filename="libMultipathSensorPlugin.so">

<ros>

<namespace>/multipath</namespace>

<remapping>~/out:=scan</remapping>

</ros>

<output_type>sensor_msgs/LaserScan</output_type>

<disableNoise>0</disableNoise>

<satNum>8</satNum>

<satAzimuth>2.93 3.39 5.56 0.23 4.02 1.38 2.23 0.65</satAzimuth>

<satElevation>0.78 0.32 1.02 0.90 0.47 0.26 1.09 0.38</satElevation>

</plugin>

##  5 Experimental Results

In this section, we examine the performance of the proposed plugin by
comparing its performance with the publicly available Urban-Nav dataset,
containing the GNSS records of the city of Hong Kong.

###  5.1 Urban-Nav Dataset

We utilize the open-source Urban-Nav Dataset collected and prepared by Hsu et.
al. [26] to test the simulation of multipath error from our plugin. The
dataset contains multiple GNSS navigation records of the city of Hong Kong and
Tokyo in the form of RINEX files measured using various GNSS receivers,
namely, u-blox ZED-F9P, EVK-M8T and NovAtel Flexpak6, etc. The ground truth
position of the receiver is computed using sensor fusion from IMU and NovAtel
SPAN-CPT receiver. Fig. 5(a) shows the map of Hong Kong where the data is
acquired which we used in our study. The recorded data captures the essential
components of an urban environment, for example, medium and tall buildings and
wide streets. Fig. 5(b) illustrates the satellite constellation diagram at the
begining of the data acquisition.

(a) Map of Hong Kong where the data is recorded.

(b) Satellite constellation diagram for u-blox EVK-M8T

Figure 5: Urban-Nav dataset for Hong Kong city: Map and Satellite
Constellation (adapted from Hsu et. al. [26])

###  5.2 Results

In Fig. 7, the GNSS position solution obtained from the plugin is depicted. We
observe a close match between the noise properties of the simulated position
fix by the proposed plugin and the measured GNSS position solution by a
commercially available receiver, u-blox EKV-M8T. The exact position fix differ
as the commercially available receivers performs additional post-processing on
the raw position fix generated by the recursive least squares estimation
algorithm by using an Extended Kalman Filter (EKF). Although there is a
mismatch in the position fix, the proposed plugin is able to effectively
capture the multipath prone regions in the city. The plugin is also able to
capture the loss of position fix during the phase of the trajectory when
vehicle passes near tall buildings. This is attributed to the obstruction of
the satellites due to the buildings.

Fig. 6 shows a typical instance of the plugin, the left hand side shows the
sky plot using a simulated fish-eye camera pointing up toward the sky, and, on
the right side a drone is equipped with our proposed plugin. The spherical
bubble around the drone shows the intensity of the multipath error at the
given location. Additionally, we have created a heat map highlighting the
intensity of the GNSS position error throughout the entire city at an altitude
of 15m as shown in Fig. 8 to understand the effects of the multipath errors.

Figure 6: An instance of the GNSS multipath plugin being attached to a drone
model in the Gazebo simulator Figure 7: Comparison of GNSS position fix of
u-blox EVK-M8T and the proposed plugin Figure 8: Heat Map depicting multipath
error in virtual Hong Kong environment

##  6 Conclusion and Future Work

In this paper, we presented the design of a Gazebo plugin for GNSS sensor in a
virtual 3D environment. Our design enables the construction of a virtual 3D
environment of a part of any city and presents a realistic high-fidelity
simulation of various commercially available GNSS receivers. We chose the ROS-
Gazebo framework because of its ubiquitous application and support in the
robotics community. The proposed plugin has a wide range of applications,
especially in the growing field of high-assurance robotics and autonomy. In
the future, we plan to extend the functionality of the proposed GNSS plugin to
incorporate the orbital equations of the satellites to ease the process of
retrieving the coordinates of the satellites in a constellation. Finally, we
plan to officially release our plugin inside the Gazebo community.

## Acknowledgments

This research is funded by the Secure Systems Research Center(SSRC) at the
Technology Innovation Institute (TII), UAE. The authors are grateful to Dr.
Shreekant (Ticky) Thakkar and his team members at the SSRC for their valuable
comments and support.

## References

  * Monteiro et al. [2005] Monteiro, L. S., Moore, T., and Hill, C., “What is the accuracy of DGPS?” _Journal of Navigation_ , Vol. 58, No. 2, 2005, p. 207–225.  10.1017/S037346330500322X.
  * Bauer et al. [2012] Bauer, S., Obst, M., and Wanielik, G., “3D environment modeling for GPS multipath detection in urban areas,” _International Multi-Conference on Systems, Signals Devices_ , 2012, pp. 1–5.  10.1109/SSD.2012.6198054.
  * Zhang et al. [2018] Zhang, S., Lo, S., Chen, Y.-H., Walter, T., and Enge, P., “GNSS multipath detection in urban environment using 3D building model,” _2018 IEEE/ION Position, Location and Navigation Symposium (PLANS)_ , 2018, pp. 1053–1058.  10.1109/PLANS.2018.8373486.
  * Lau and Cross [2007] Lau, L., and Cross, P., “Development and testing of a new ray-tracing approach to GNSS carrier-phase multipath modelling,” _Journal of Geodesy_ , Vol. 81, No. 11, 2007, pp. 713–732.  10.1007/s00190-007-0139-z.
  * Lee et al. [2004] Lee, H. K., Lee, J. G., and Jee, G.-I., “GPS multipath detection based on sequence of successive-time double-differences,” _IEEE Signal Processing Letters_ , Vol. 11, No. 3, 2004, pp. 316–319.  10.1109/LSP.2003.821744.
  * Wu and Hsieh [2010] Wu, J., and Hsieh, C.-H., “Statistical modeling for the mitigation of GPS multipath delays from day-to-day range measurements,” _Journal of Geodesy_ , Vol. 84, No. 4, 2010, pp. 223–232.
  * Brenner et al. [1998] Brenner, M., Reuter, R., and Schipper, B., “GPS landing system multipath evaluation techniques and results,” _Proceedings of the 11th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 1998)_ , 1998, pp. 999–1008.
  * Wang et al. [2018] Wang, Y., Chen, X., and Liu, P., “Statistical Multipath Model Based on Experimental GNSS Data in Static Urban Canyon Environment,” _Sensors_ , Vol. 18, No. 4, 2018.  10.3390/s18041149, URL https://www.mdpi.com/1424-8220/18/4/1149.
  * Kong [2011] Kong, S.-H., “Statistical Analysis of Urban GPS Multipaths and Pseudo-Range Measurement Errors,” _IEEE Transactions on Aerospace and Electronic Systems_ , Vol. 47, No. 2, 2011, pp. 1101–1113.  10.1109/TAES.2011.5751245.
  * Steingass et al. [2004] Steingass, E., Aerospace, G., and German, A. L., “Measuring the navigation multipath channel - a statistical analysis,” _in Proceedings of the ION GNSS 2004_ , 2004.
  * Lehner and Steingass [2005] Lehner, A., and Steingass, A., “A Novel Channel Model for Land Mobile Satellite Navigation,” _ION GNSS, Long Beach, California, 13-16 September, 2005_ , Vol. Seesio, 2005.  URL https://elib.dlr.de/19171/, lIDO-Berichtsjahr=2005, monograph_id=Paper #4,.
  * Parkinson and Spilker [1996] Parkinson, B., and Spilker, J., _Global Positioning System: Theory and Applications_ , No. v. 1 in Ciencia militar y naval, American Institute of Aeronautics & Astronautics, 1996.  URL https://books.google.com/books?id=lvI1a5J_4ewC.
  * Klobuchar [1987] Klobuchar, J. A., “Ionospheric Time-Delay Algorithm for Single-Frequency GPS Users,” _IEEE Transactions on Aerospace and Electronic Systems_ , Vol. AES-23, No. 3, 1987, pp. 325–331.  10.1109/TAES.1987.310829.
  * Bevis et al. [1994] Bevis, M., Businger, S., Chiswell, S., Herring, T. A., Anthes, R. A., Rocken, C., and Ware, R. H., “GPS Meteorology: Mapping Zenith Wet Delays onto Precipitable Water,” _Journal of Applied Meteorology and Climatology_ , Vol. 33, No. 3, 1994, pp. 379 – 386.  10.1175/1520-0450(1994)033<0379:GMMZWD>2.0.CO;2.
  * Enge and Van Dierendonck [1996] Enge, P. K., and Van Dierendonck, A., “Wide area augmentation system,” _Progress in Astronautics and Aeronautics_ , Vol. 164, 1996, pp. 117–142.
  * ope [2022] “OpenStreetMap Foundation: Openstreetmap,” , 2022.  URL https://www.openstreetmap.org/.
  * Noll [2015] Noll, C., “CDDIS Global Data Center Technical Report 2014,” _IGS Central Bureau_ , 2015, p. 121.
  * Ray et al. [1999] Ray, J. K., Cannon, M. E., and Fenton, P. C., “Mitigation of static carrier-phase multipath effects using multiple closely spaced antennas,” _Navigation_ , Vol. 46, No. 3, 1999, pp. 193–201.
  * Martin et al. [2020] Martin, A., Parry, M., Soundy, A. W., Panckhurst, B. J., Brown, P., Molteno, T. C., and Schumayer, D., “Improving real-time position estimation using correlated noise models,” _Sensors_ , Vol. 20, No. 20, 2020, p. 5913\.
  * Soundy et al. [2020] Soundy, A. W., Panckhurst, B. J., Brown, P., Martin, A., Molteno, T. C., and Schumayer, D., “Comparison of Enhanced Noise Model Performance Based on Analysis of Civilian GPS Data,” _Sensors_ , Vol. 20, No. 21, 2020, p. 6050\.
  * Isik et al. [2020] Isik, O. K., Hong, J., Petrunin, I., and Tsourdos, A., “Integrity Analysis for GPS-Based Navigation of UAVs in Urban Environment,” _Robotics_ , Vol. 9, No. 3, 2020.  10.3390/robotics9030066, URL https://www.mdpi.com/2218-6581/9/3/66.
  * Langley [1999] Langley, R. B., “Dilution of Precision,” GPS World, 10(5):52–59, 1999\.
  * Quigley et al. [2009] Quigley, M., Gerkey, B., Conley, K., Faust, J., Foote, T., Leibs, J., Berger, E., Wheeler, R., and Ng, A., “ROS: an open-source Robot Operating System,” _Proc. of the IEEE Intl. Conf. on Robotics and Automation (ICRA) Workshop on Open Source Robotics_ , 2009.
  * Koenig and Howard [2004] Koenig, N., and Howard, A., “Design and use paradigms for Gazebo, an open-source multi-robot simulator,” _2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (IEEE Cat. No.04CH37566)_ , Vol. 3, 2004, pp. 2149–2154 vol.3.  10.1109/IROS.2004.1389727.
  * Möller and Landskron [2019] Möller, G., and Landskron, D., “Atmospheric bending effects in GNSS tomography,” _Atmospheric Measurement Techniques_ , Vol. 12, No. 1, 2019, pp. 23–34.
  * Hsu et al. [2021] Hsu, L.-T., Kubo, N., Wen, W., Chen, W., Liu, Z., Suzuki, T., and Meguro, J., “UrbanNav:an open-sourced multisensory dataset for benchmarking positioning algorithms designed for urban areas,” _Proceedings of the 34th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2021)_ , 2021.  10.33012/2021.17895.


$\displaystyle\left(\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\left(\widetilde{D}_{t}-\mathbb{E}(\widetilde{D}_{t})\right)\right)$
$\displaystyle=\mathbb{E}\left(\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\widetilde{D}_{t}\right)-\mathbb{E}\left(\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\right)\mathbb{E}\left(\widetilde{D}_{t}\right)$
$\displaystyle=\mathbb{E}\left(\widetilde{D}_{t}|\,d\leq\widetilde{D}_{t}\right)\mathbb{E}\left(\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\right)-\mathbb{E}\left(\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\right)\mathbb{E}\left(\widetilde{D}_{t}\right)$
$\displaystyle=\left(\mathbb{E}\left(\widetilde{D}_{t}|\,d\leq\widetilde{D}_{t}\right)-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\mathbb{P}\left(d\leq\widetilde{D}_{t}\right)\geq
0,$

since
$\mathbb{E}(\widetilde{D}_{t}|\,d\leq\widetilde{D}_{t})-\mathbb{E}(\widetilde{D}_{t})\geq
0$ for $d\in[\underline{d},\,\overline{d}]$. To see that the weights integrate
to one, note that

$\displaystyle\widetilde{D}_{t}$
$\displaystyle=\underline{d}+\int_{\underline{d}}^{\widetilde{D}_{t}}\mathrm{d}\widetilde{d}=\underline{d}+\int_{\underline{d}}^{\overline{d}}\mathbf{1}\left\\{\widetilde{d}\leq\widetilde{D}_{t}\right\\}\mathrm{d}\widetilde{d}.$

Using this we have,

$\displaystyle 1=\mathrm{Var}\left(\widetilde{D}_{t}\right)$
$\displaystyle=\mathbb{E}\left[\left(\widetilde{D}_{t}-\underline{d}\right)\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\right]$
$\displaystyle\quad-\mathbb{E}\left[\mathbb{E}\left(\int_{\underline{d}}^{\overline{d}}\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\mathrm{d}d\right)\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\right]$
$\displaystyle=\mathbb{E}\left[\left(\widetilde{D}_{t}-\underline{d}\right)\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\right]$
$\displaystyle\quad-\mathbb{E}\left(\int_{\underline{d}}^{\overline{d}}\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\mathrm{d}\widetilde{d}\right)\mathbb{E}\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)$
$\displaystyle=\mathbb{E}\left[\left(\widetilde{D}_{t}-\underline{d}\right)\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\right]$
$\displaystyle=\mathbb{E}\left[\int_{\underline{d}}^{\overline{d}}\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\mathrm{d}\widetilde{d}\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\right]$
$\displaystyle=\int_{\underline{d}}^{\overline{d}}\mathbb{E}\left(\mathbf{1}\left\\{d\leq\widetilde{D}_{t}\right\\}\left(\widetilde{D}_{t}-\mathbb{E}\left(\widetilde{D}_{t}\right)\right)\right)\mathrm{d}d.$

Finally, the fact that $\Delta_{t}$ is decreasing in $\sigma_{e,t}^{2}$
follows directly from (A.3) which continues to hold when Assumption 4 is
replaced by Assumption 7. $\square$

###  A.2 Proof of Theorem 2

The result follows directly from Theorem 1 and Lemma A.1. $\square$

###  A.3 Proof of Theorem 3

We have

$\displaystyle\widehat{\beta}_{\mathrm{ES}}$
$\displaystyle=\frac{\sigma_{D}^{-2}T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}-\overline{D}\right)\left(Y_{t}-\overline{Y}\right)}{\sigma_{D}^{-2}T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}-\overline{D}\right)^{2}}$
$\displaystyle=\frac{T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}^{*}-\overline{D}^{*}\right)\left(Y_{t}^{*}-\overline{Y}^{*}\right)}{T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}^{*}-\overline{D}^{*}\right)^{2}}$
$\displaystyle\overset{\mathbb{P}}{\rightarrow}\int_{0}^{1}c\left(D^{*},\,Y^{*},\,s\right)\mathrm{d}s$

by Assumption 8. Given the definition of $c\left(D^{*},\,Y^{*},\,s\right)$,
the statements of the theorem then follow from the same arguments as in the
proof of Theorem 1. $\square$

###  A.4 Proof of Theorem 4

Note that Theorem 3 is obtained under Assumption 7 while the current theorem
uses Assumption 4. The difference between the two assumptions is that
$\mathbb{E}(g_{D,e}\left(e_{t},\,t\right)^{2})<\infty$ in Assumption 7 while
$\mathbb{E}(g_{D,e}\left(e_{t},\,t\right)^{2})\rightarrow\infty$ as
$\sigma_{e,t}^{2}\rightarrow\infty$ in Assumption 4. Given that we take the
limits sequentially, we can use the limiting result from Theorem 3 before
taking the limit as $\sigma_{e,t}^{2}\rightarrow\infty$. Thus, the theorem
follows directly from Theorem 3 and Lemma A.1. $\square$

###  A.5 Proof of Theorem 5

We have

$\displaystyle\widehat{\beta}_{\mathrm{ES}}-\beta_{\mathrm{ES}}$
$\displaystyle=\frac{T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}^{*}-\overline{D}^{*}\right)\left(Y_{t}^{*}-\overline{Y}^{*}-\beta_{\mathrm{ES}}\left(D_{t}^{*}-\overline{D}^{*}\right)\right)}{T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}^{*}-\overline{D}^{*}\right)^{2}}=\frac{T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}^{*}-\overline{D}^{*}\right)\varepsilon_{t}}{T_{P}^{-1}\sum_{t=1}^{T_{P}}\left(D_{t}^{*}-\overline{D}^{*}\right)^{2}}.$

Using Assumption 8(ii) and Assumption 9, we obtain the statement of the
theorem. $\square$

###  A.6 Proof of Corollary 1

The result follows directly from Theorem 5 and Lemma A.1. $\square$


Continuum-Armed Bandits: A Function Space Perspective

Shashank Singh

<EMAIL_ADDRESS>Carnegie Mellon University

###### Abstract

Continuum-armed bandits (a.k.a., black-box or $0^{th}$-order optimization)
involves optimizing an unknown objective function given an oracle that
evaluates the function at a query point, with the goal of using as few query
points as possible. In the most well-studied case, the objective function is
assumed to be Lipschitz continuous and minimax rates of simple and cumulative
regrets are known in both noiseless and noisy settings. This paper studies
continuum-armed bandits under more general smoothness conditions, namely Besov
smoothness conditions, on the objective function. In both noiseless and noisy
conditions, we derive minimax rates under simple and cumulative regrets. Our
results show that minimax rates over objective functions in a Besov space are
identical to minimax rates over objective functions in the smallest Hölder
space into which the Besov space embeds.

##  1 Introduction

Many computational problems involve optimizing an unknown function of
continuous inputs given only the ability to evaluate the function at a point.
In settings where the function is expensive (in terms of computation, time,
money, or another resource) to evaluate, it is desirable to estimate the
global optimum of the function using as few evaluations as possible. This
problem, which we refer to as _continuum-armed bandits_ (following Agrawal
(1995)), is also known as _black-box_ (Golovin et al., 2017), _derivative-
free_ (Larson et al., 2019), or _$0^{th}$ -order_ (Golovin et al., 2019)
optimization, due to the unavailability of analytic information, including
gradients, about the function.

Recent interest in continuum-armed bandits has been fueled by a need to tune
hyperparameters of machine learning models and other complex algorithms (Li et
al., 2017; Kandasamy, 2018; Feurer and Hutter, 2019). However, the problem is
a natural framework for improving many scientific and engineering processes as
diverse as scientific experimentation (Nakamura et al., 2017), taxonomy
matching (Pandey et al., 2007), and baking cookies (Kochanski et al., 2017).

To make continuum-armed bandits tractable, most work assumes the “arms” (i.e.,
values being optimized over) lie in a compact set
$\mathcal{X}\subset\mathbb{R}^{d}$ and makes smoothness assumptions on the
objective function. For example, much work has focused on the case where the
objective function lies in a Lipschitz or Hölder space over $\mathcal{X}$.
However, compared to our understanding of closely related estimation problems
such as nonparametric regression, our understanding of the relationship
between the assumptions made on the objective function and the difficulty of
optimization is limited.

This paper considers a large family of function spaces, called Besov spaces,
which include previously considered function spaces, such as Hölder and
Sobolev spaces, as special cases. Classical results for nonparametric
regression shown that convergence rates over Besov spaces can sometimes be
improved by used special “spatially adaptive” methods, such as certain wavelet
or spline methods, instead of commonly used “linear” methods, such as most
kernel methods (Donoho et al., 1995; Donoho and Johnstone, 1998). Our primary
motivation is to investigate whether such an improvement is possible for
continuum-armed bandits. To understand this, we prove new lower bounds on
minimax convergence rates for continuum-armed bandits over Besov spaces under
a variety of conditions. Our lower bounds imply that minimax rates over Besov
spaces are identical to those obtained by embedding the Besov space in a
Hölder space, and we conclude that, unlike in regression, Besov smoothness
provides no advantage over Hölder smoothness, clarifying the nature of the
smoothness properties of a function that determine the difficulty of its
optimization. Additionally, as we discuss in Section 3, our results
generalize, unify, and strengthen several prior results.

##### Organization

The remainder of this paper is organized as follows: Section 2 provides
necessary notation, defines relevant function spaces, and formalizes the
optimization problem we study. Section 3 reviews related work, primarily work
in the Hölder and reproducing kernel Hilbert space (RKHS) cases. Sections 4
and 5 state and prove our main results, for the noiseless and noisy cases,
respectively, and Section 6 concludes with some discussion.

##  2 Notation, Problem Setup, and Technical Background

##### Notation

For non-negative real sequences $\\{a_{n}\\}_{n\in\mathbb{N}}$,
$\\{b_{n}\\}_{n\in\mathbb{N}}$, $a_{n}\lesssim b_{n}$ indicates
$\limsup_{n\to\infty}\frac{a_{n}}{b_{n}}<\infty$, and $a_{n}\asymp b_{n}$
indicates $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$. For a real number
$\sigma$, $\lfloor\sigma\rfloor:=\max\\{s\in\mathcal{Z}:s<\sigma\\}$ denotes
the greatest integer strictly less than $\sigma$. For non-negative integers
$k\in\mathbb{N}$, $[k]:=\\{1,...,k\\}$ denotes the set of positive integers
$\leq k$. For the remainder of this paper, we fix the domain of the
optimization problem to be the $d$-dimensional unit cube, denoted
$\mathcal{X}:=[0,1]^{d}$.

###  2.1 Hölder and Besov Spaces

We now define the main function spaces and norms discussed in this paper, and
review basic relationships between these spaces. We begin with Hölder spaces:

######  Definition 1 (Hölder Seminorm, Space, and Ball).

Fix $s\in[0,\infty)$. When $s$ is an integer, $\mathcal{C}^{s}(\mathcal{X})$
denotes the class of $s$-times differentiable functions in
$\mathcal{L}_{\infty}(\mathcal{X})$. Otherwise, the _Hölder seminorm_
$\|\cdot\|_{\mathcal{C}^{s}(\mathcal{X})}:\mathcal{C}^{\lfloor
s\rfloor}\to[0,\infty]$ is defined by

$\|f\|_{\mathcal{C}^{s}(\mathcal{X})}=\sup_{\beta\in\mathbb{N}^{d}:\|\beta\|=\lfloor
s\rfloor}\sup_{x\neq
y\in\mathcal{X}}\frac{\left|f^{\beta}(x)-f^{\beta}(y)\right|}{\|x-y\|_{2}^{s-\lfloor
s\rfloor}},$

for all $f\in\mathcal{C}^{\lfloor s\rfloor}(\mathcal{X})$, where $f^{\beta}$
denotes the $\beta^{th}$ mixed derivative of $f$. The _Hölder space_
$\mathcal{C}^{s}(\mathcal{X})$ is defined by

$\mathcal{C}^{s}(\mathcal{X}):=\\{f\in\mathcal{L}_{\infty}:\|f\|_{\mathcal{C}^{s}(\mathcal{X})}<\infty\\}$

and, for $L>0$, the _Hölder ball_ $\mathcal{C}^{s}(\mathcal{X},L)$ is defined
by

$\mathcal{C}(\mathcal{X},L):=\\{f\in\mathcal{L}_{\infty}:\|f\|_{\mathcal{C}^{s}(\mathcal{X})}\leq
L\\}.$

We now turn to defining Besov spaces. While many equivalent definitions of
Besov spaces exist, we use a definition in terms of orthonormal wavelet bases
that has classically been used in the context of nonparametric statistics
(Donoho et al., 1995). This wavelet basis can be constructed from a “father”
wavelet $\phi\in\mathcal{L}_{\infty}(\mathcal{X})$ by defining, for each
$j\in\mathbb{N}$ and $\lambda\in[2^{j}]^{d}$ (i.e., $\lambda$ is a $d$-tuple
of integers between $1$ and $2^{j}$),

$\phi_{j,\lambda}(x):=\left\\{\begin{array}[]{ccl}2^{dj/2}\phi(2^{dj}x-\lambda)&:&\text{if
}2^{dj}x-\lambda\in\mathcal{X}\\\ 0&:&\text{ else}.\end{array}\right.$

for all $x\in\mathcal{X}$. In particular, the set
$\\{\phi_{j,\lambda}\\}_{j\in\mathbb{N},\lambda\in[2^{j}]^{d}}$ is a basis of
$\mathcal{L}^{2}(\mathcal{X})$. We will additionally require the basis to be
orthonormal and to satisfy a certain regularity condition called
“$r$-regularity”, for which the father wavelet must be carefully constructed;
we refer the reader to Meyer (1992) or Daubechies (1992) for standard
constructions. For this paper, the exact wavelet construction is not crucial
as it affects only constant factors in our rates, and it suffices to note that
(a) such an orthonormal basis exists, (b) it can be indexed as
$\left\\{\phi_{j,\lambda}\right\\}_{j\in\mathbb{N},\lambda\in[2^{j}]^{d}}$,
and (c), for each $j\in\mathbb{N}$,
$\lambda\neq\lambda^{\prime}\in[2^{j}]^{d}$, $\phi_{j,\lambda}$ and
$\phi_{j,\lambda^{\prime}}$ have disjoint supports. Given such a wavelet
basis, we can define Besov spaces and norms as follows:

######  Definition 2 (Besov Norm, Space, and Ball).

Fix $\sigma\in(0,\infty)$ and $p,q\in[1,\infty]$. Given a wavelet basis
$\\{\phi_{j,\lambda}\\}_{j\in\mathbb{N},\lambda\in[2^{j}]^{d}}$ of
$\mathcal{L}_{2}(\mathcal{X})$, the _Besov norm_
$\|\cdot\|_{\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})}:\mathcal{L}_{2}(\mathcal{X})\to[0,\infty]$
is defined by

$\|f\|_{\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})}=\left(\sum_{j\in\mathbb{N}}2^{jq(\sigma+d\frac{p-2}{2p})}\left(\sum_{\lambda\in[2^{j}]^{d}}f_{j,\lambda}^{p}\right)^{\frac{q}{p}}\right)^{\frac{1}{q}},$
(1)

for all $f\in\mathcal{L}_{2}(\mathcal{X})$, where $f_{j,\lambda}:=\langle
f,\phi_{j,\lambda}\rangle_{\mathcal{L}_{2}(\mathcal{X})}$ denotes the
$(j,\lambda)^{th}$ coefficient in the wavelet expansion of $f$. The _Besov
space_ $\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})$ is defined by

$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X}):=\\{f\in\mathcal{L}_{2}:\|f\|_{\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})}<\infty\\}$

and, for $L>0$, the _Besov ball_ $\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)$
is defined by

$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L):=\\{f\in\mathcal{L}_{2}:\|f\|_{\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})}\leq
L\\}.$

Besov spaces are appealing to work with in part because they provide a unified
representation of many commonly used function spaces. Examples include Hölder
spaces
$\mathcal{C}^{s}(\mathcal{X})=\mathcal{B}_{\infty,\infty}^{s}(\mathcal{X})$,
Hilbert-Sobolev spaces $\mathcal{B}_{2,2}^{s}(\mathcal{X})$ (which also
correspond to the Matérn kernel RKHS of order $s-d/2$), and the space BV of
functions of bounded variation, which satisfies
$\mathcal{B}_{1,1}^{1}\subseteq\text{BV}\subseteq\mathcal{B}_{1,\infty}^{1}$
111As will become clear from our main results, the parameter $q$ of the Besov
space does not affect minimax convergence rates, and so an inclusion of the
form
$\mathcal{B}_{p,1}^{\sigma}\subseteq\mathcal{S}\subseteq\mathcal{B}_{p,\infty}^{\sigma}$
is effectively an equivalence for our purposes.. For standard references on
Besov spaces, see Sawano (2018) or Chapter 14 of Leoni (2017).

Figure 1: Three example functions, $f,g,h:\mathcal{X}=[-2.5,2.5]\to[0,1]$ with
maxima at $x=0$. $f$ and $g$ have the same Lipschitz norm
($\|f\|_{\mathcal{C}^{1}(\mathcal{X})}=\|g\|_{\mathcal{C}^{1}(\mathcal{X})}$),
but $g$ has a larger Besov $\mathcal{B}_{p,q}^{1}$ norm
($\|f\|_{\mathcal{B}_{p,p}^{1}(\mathcal{X})}<\|g\|_{\mathcal{B}_{p,p}^{1}(\mathcal{X})}$)
for $p<\infty$. The function $h$ has infinite Lipschitz norm
($\|h\|_{\mathcal{C}^{1}(\mathcal{X})}=\infty$), but has finite
$\tfrac{1}{2}$-Hölder and $1$-Sobolev norms
($\|h\|_{\mathcal{C}^{1/2}(\mathcal{X})},\|h\|_{\mathcal{B}^{1}_{2,2}(\mathcal{X})}<\infty$).

To illustrate the relative behaviors of Hölder and other Besov smoothness
conditions, Figure 1 shows an example of two functions with identical
Lipschitz norm and different Besov norms. Compared to Hölder smoothness, which
is determined only by the least smooth part of a function, Besov smoothness is
determined by a function’s integrated smoothness over the entire domain. For
example, for a differentiable function, the total variation norm is equivalent
to the average magnitude of the gradient of the function (Leoni, 2017).
However, since both norms ultimately measure smoothness, Besov norms also
exhibit certain relationships with Hölder norms, exemplified by the following
continuous embedding of Besov spaces into Hölder spaces, which will be
utilized later in this paper:

######  Theorem 3 (Special Case of Theorem 1.15, Haroske and Schneider
(2009)).

Let $\sigma\in(d/p,\infty)$, $p,q\in[1,\infty]$. Then,
$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})\subseteq\mathcal{C}^{\sigma-d/p}(\mathcal{X})$
and

$\sup_{f\in\mathcal{B}_{p,q}^{\sigma}(\mathcal{X}):f\neq
0}\frac{\|f\|_{\mathcal{C}^{\sigma-d/p}(\mathcal{X})}}{\|f\|_{\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})}}<\infty.$

Moreover, for any $s>\sigma-d/p$,
$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})\not\subseteq\mathcal{C}^{s}(\mathcal{X})$.

###  2.2 Problem Setup

Fix real values $p,q\in[1,\infty]$ and $s\in(0,\infty)$, and fix a Besov class
$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)$ over $X:=[0,1]^{d}$. Fix an
unknown function $f\in\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)$. To allow
randomized optimization algorithms, also fix a random variable $Z$ taking
values in a set $\mathcal{Z}$. Finally, fix an evaluation time horizon
$T\in\mathbb{N}$ (which may be known or unknown).

Formally, for $t\in[T]$, a sequential $0^{th}$-order optimization procedure
$\widehat{X}=(\widehat{X}_{1},...,\widehat{X}_{T})$ can be characterized as a
sequence of queries
$\widehat{X}_{n+1}:\mathcal{Z}\times\mathcal{X}^{t}\times\mathbb{R}^{d}\to\mathcal{X}$
mapping the previous queries $(X_{1},...,X_{t})$ and their observed costs
$Y_{1}=f(X_{1}),...,Y_{t}=f(X_{t})$ to a new query point. For notational
compactness, we will simply write
$\widehat{X}_{t+1}=\widehat{X}_{t+1}(Z,X_{1},...,X_{t},Y_{1},...,Y_{t})$, with
the understanding that $\widehat{X}_{t+1}$ may depend on algorithmic
randomness $Z$, previous queries $X_{1},...,X_{t}$, and observations
$Y_{1},...,Y_{t}$.

In this paper, we will separately consider two problem settings, namely the
noiseless case, in which we query the objective $f$ directly, and the noisy
case, in which observations are corrupted by IID additive noise.

In the noiseless case, the goal is to minimize either the expected simple
regret

$R_{S}(\widehat{X},f):=\mathop{\mathbb{E}}_{Z}\left[\sup_{x\in\mathcal{X}}f(x)-f(X_{T})\right]$

or the expected cumulative regret

$R_{C}(\widehat{X},f):=\mathop{\mathbb{E}}_{Z}\left[\sum_{n=1}^{T}\sup_{x\in\mathcal{X}}f(x)-f(X_{n})\right].$

Specifically, we are interested in bounding (up to constant factors) the
minimax risk

$\displaystyle M_{S,T}(\mathcal{F})$
$\displaystyle:=\inf_{\widehat{x},Z}\sup_{f\in\mathcal{F}}R_{S}(\widehat{x}_{T},f)$
$\displaystyle\text{ or }\quad M_{C,T}(\mathcal{F})$
$\displaystyle:=\inf_{\widehat{x},Z}\sup_{f\in\mathcal{F}}R_{C}(\widehat{X}_{T},f),$

where the infima are taken over all optimization procedures $\widehat{X}$ and
randomness sources $Z$.

In the noisy case, we additionally fix a noise probability distribution $P$,
and, rather than observing the true value of $f$ at the query point, our
observations are corrupted by additive noise sampled IID from $P$; i.e.,
$Y_{t}=f(X_{t})+\epsilon_{t}$, where
$\epsilon_{1},...,\epsilon_{T}\stackrel{{\scriptstyle IID}}{{\sim}}P$. Our
goal is again to minimize either the expected simple regret

$R_{S}(\widehat{X},f,P):=\mathop{\mathbb{E}}_{Z,\epsilon}\left[\sup_{x\in\mathcal{X}}f(x)-f(X_{T})\right]$

or the expected cumulative regret

$R_{C}(\widehat{X},f,P):=\mathop{\mathbb{E}}_{Z,\epsilon}\left[\sum_{n=1}^{T}\sup_{x\in\mathcal{X}}f(x)-f(X_{n})\right],$

where the expectations are now additionally taken with respect to the noise
$\epsilon$. Assuming further that the unknown distribution $P$ lies in some
regularity class $\mathcal{P}$, we are interested in characterizing the
minimax risk

$\displaystyle M_{S,T}(\mathcal{F},\mathcal{P})$
$\displaystyle:=\inf_{\widehat{x},Z}\sup_{f\in\mathcal{F},P\in\mathcal{P}}R_{S}(\widehat{X}_{T},f,P)$
$\displaystyle\text{ or }\quad M_{C,T}(\mathcal{F},\mathcal{P})$
$\displaystyle:=\inf_{\widehat{x},Z}\sup_{f\in\mathcal{F},P\in\mathcal{P}}R_{C}(\widehat{X}_{T},f,P).$

In particular, we shall consider the case when
$\mathcal{P}=\mathsf{subG}(\eta^{2})$ is the class of sub-Gaussian random
variables with variance proxy $\eta^{2}$.

##  3 Related Work

Bandit optimization of smooth functions is a well-studied problem with a
number of theoretical and practical results, with perhaps the earliest results
due to Kiefer and Wolfowitz (1952), who proposed a predecessor of stochastic
gradient descent based on finite-differences. Here, we review the most
relevant modern work, consisting primarily of work on the cases of objective
functions lying in a Hölder space or RKHS.

First, we note that, under _simple regret_ , upper bounds can be obtained
using “pure exploration” algorithms that randomly sample $\mathcal{X}$, apply
standard (e.g., local polynomial (Fan and Gijbels, 1996) or wavelet (Donoho
and Johnstone, 1998)) nonparametric regression methods, and then compute the
optimum of the estimated regression function. In the Hölder case, it has long
been known that pure exploration algorithms can achieve the minimax optimal
rate, and recent work has therefore focused on showing that more sophisticated
adaptive algorithms strictly outperform pure exploration under stronger
assumptions (e.g., local convexity or other shape constraints near the
optimum) on the objective function (Minsker, 2012, 2013; Grill et al., 2015;
Wang et al., 2018a). Our lower bounds confirm that, in both noiseless and
noisy conditions, in terms of simple regret, pure exploration continues to be
minimax rate-optimal under more general Besov smoothness assumptions, at least
in the absence of further shape constraints.

Under _cumulative regret_ , a number of papers have studied the case when the
objective function lies in the Hölder class $\mathcal{C}^{s}$ (Kleinberg,
2005; Kleinberg et al., 2008; Bubeck et al., 2011a, b) (see Chapter 4 of
Slivkins (2019) for further review). For example, in the $s$-Hölder case with
$s\in(0,1]$, if $\mathcal{P}$ is any class of uniformly bounded noise
distributions, Bubeck et al. (2011a) showed upper bounds on cumulative regret
of order

$M_{C,T}\left(\mathcal{C}^{s}(\mathcal{X},L),\mathcal{P}\right)\lesssim
T^{\frac{s+d}{2s+d}}\left(\log T\right)^{\frac{1}{2s+d}}.$ (2)

This rate matches lower bounds such as ours in $T$, which, together with the
embedding Theorem 3, implies that the rates are minimax-optimal. Our lower
bounds for the noisy case match (2) when $s\in(0,1]$, showing that this rate
is optimal not only for Hölder classes but also for more general Besov
classes. We note that Bubeck et al. (2011a) additionally showed that, for the
upper bound (2), Lipschitz smoothness (i.e., with $s=1$) is needed only in
weak local form, near the global optimum.

To the best of our knowledge, under cumulative regret, the $s$-Hölder case
with $s>1$ is open, in that a gap remains existing lower and upper
bounds222For the significantly different setting where the objective function
is strongly convex, Akhavan et al. (2020) recently derived a dimension-
independent minimax cumulative regret rate of $\asymp T^{1/s}$ that applies
for $s>1$, in the presence of noise.. Our lower bounds extend smoothly from
$s\in(0,1]$ to $s>1$, and we conjecture that our lower bounds are tight for
all $s>0$. However, existing upper bounds for $s>1$ are polynomially larger in
$T$. In particular, Grant and Leslie (2020) consider the Hölder case when $s$
is an integer (i.e., when the $(s-1)^{st}$ derivative of the objective
function is Lipschitz). Under a particular prior on the objective function,
they upper bounded the Bayesian cumulative risk of a Thompson Sampling
(Thompson, 1933) algorithm by order $T^{\frac{2s^{2}+7s+1}{4s^{2}+6s+2}}$. To
the best of our knowledge, these are the only known upper bounds on cumulative
regret that improve with $s>1$, and they approach the parametric rate of
$T^{1/2}$ as $s\to\infty$. However these rates are polynomially worse than the
minimax rate for all finite $s$, and, moreover, they are proven for Bayesian
risk, a weaker notion of risk than worst-case/minimax risk (Wasserman, 2013).
333Concurrent with the present article, Liu et al. (2020) provide an
algorithm, based on running a UCB meta-algorithm over a collection of local
polynomial bandit algorithms, and bound the worst-case risk of this algorithm
in the presence of noise for all Hölder exponents $s>0$. Their upper bound
exceeds our lower bound by a multiplicative factor of $\left(\log
T\right)^{\frac{5s+2d}{2s+d}}$, providing the minimax rate, up to
polylogarithmic factors in $T$, for all $s>0$.

Beyond the Hölder setting, work has largely focused on the setting in which
$\mathcal{F}$ is an RKHS. Due to the equivalence between RKHSs of Matérn
kernels and Hilbert-Sobolev spaces (see, e.g., Example 2.6 of Kanagawa et al.
(2018)) this intersects with our results in the special case $p=q=2$. In
particular, for $\sigma>d/2$, in the absence of noise, Bull (2011) obtained a
simple regret rate

$M_{S,T}\left(\mathcal{B}_{2,2}^{\sigma}(\mathcal{X},L)\right)\asymp
T^{1/2-\sigma/d},$ (3)

while, in the presence of noise, Scarlett et al. (2017) proved lower bounds of
orders

$M_{S,T}\left(\mathcal{B}_{2,2}^{\sigma}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)\gtrsim\left(\frac{\eta^{2}}{T}\right)^{\frac{\sigma-d/2}{2\sigma}}$

and

$M_{C,T}\left(\mathcal{B}_{2,2}^{\sigma}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)\gtrsim\eta^{\frac{\sigma-d/2}{\sigma}}T^{\frac{\sigma+d/2}{2\sigma}},$

under simple and cumulative regrets, respectively. Our results generalize
these lower bounds to other Besov classes with $p,q\neq 2$, and, in the noisy
case, tighten these bounds (by polylogarithmic factors in $T$ and $\eta$) to
the correct minimax rate. Moreover, although this was less apparent in the
context of RKHSs (and appears to have been missed by Scarlett et al. (2017)),
expressing this problem in the language of Besov spaces makes it clear that
these rates can be directly compared to those for the Hölder case. Prior to
this, the best known upper bound appears to have been a rate of
$T^{\frac{\sigma+d^{2}+d/2}{2\sigma+d^{2}}}$ due to Srinivas et al. (2010)
based on a Gaussian Process UCB algorithm; this rate is sub-optimal by
polynomial factors in $T$.

On a more historical note, one of the primary motivators for this work is a
series of classical results due to Nemirovski (1985), Donoho and Johnstone
(1998), and others, who investigated nonparametric _estimation_ problems, such
as regression and density estimation, over Besov spaces. They showed that
minimax rates for nonparametric estimation over
$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})$ are significantly faster than over
$\mathcal{C}^{\sigma-d/p}(\mathcal{X})$. However, in contrast to minimax
optimality over $\mathcal{C}^{\sigma-d/p}(\mathcal{X})$, for which simpler
methods (e.g., kernel, basis, or spline regression) are optimal, leveraging
accelerated convergence rates over $\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})$
requires relatively complex methods (such as thresholding wavelet methods or
spatially adaptive smoothing splines) that are “spatially adaptive”, in the
sense of devoting selectively greater representational power to less smooth
parts of the estimand. The equivalence we show for _optimization_ over Besov
and Hölder spaces shows that this phenomenon fails to carry over from
estimation; in particular, spatial adaptivity confers no asymptotic benefit
for bandit optimization over Besov spaces. This may be related to the finding
of Bubeck et al. (2011a) that Lipschitz smoothness may be needed only in a
weak sense, locally near the global optimum; i.e., allocation of
representational power to the function near a point $x$ should correspond to
the plausibility of that point being a global optimum, rather than to relative
smoothness of $f$ near that point.

Finally, we note that many papers have also considered stronger shape
assumptions such as convexity of the objective function, either globally or
near an optimum Auer et al. (2007); Cope (2009); Bubeck et al. (2018); Wang et
al. (2018b); Golovin et al. (2019). In contrast to the rates discussed above
and elsewhere in this paper, which scale exponentially with the number $d$ of
optimization variables, rates in these cases scale much more favorably
(polynomially) with $d$. However, convexity and related shape assumptions may
be too strong for many black-box optimization settings.

To summarize, our main contributions are:

  1. 1.

We identify upper bounds on the regret of existing algorithms under Besov
smoothness conditions by embedding the Besov space in a Hölder space of lower
smoothness. For many Besov spaces, these bounds are the first available,
while, for others, such as Hilbert-Sobolev spaces, they improve on state-of-
the-art results.

  2. 2.

We derive novel information-theoretic lower bounds on minimax regret under
Besov smoothness conditions, generalizing and tightening a number of existing
lower bounds. These lower bounds match the above upper bounds, thereby
identifying the minimax rate for these problems.

These results are provided in both noiseless and noisy conditions, under both
simple and cumulative losses, providing a comprehensive picture for this class
of problems. Moreover, as discussed in Section 6, a number of big-picture
conclusions follow from our results.

##  4 Noiseless Case

Here, we state and prove our main results, identifying minimax rates of simple
and cumulative regret over Besov spaces, for the noiseless case.

######  Theorem 4 (Minimax Rates, Noiseless Case).

For $\sigma>d/p$,

$\displaystyle M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)$
$\displaystyle\asymp
M_{S,T}\left(\mathcal{C}^{\sigma-d/p}(\mathcal{X},L)\right)$
$\displaystyle\asymp T^{1/p-\sigma/d},$

and, for $\sigma\in(d/p,d/p+1]$,

$\displaystyle M_{C,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)$
$\displaystyle\asymp
M_{C,T}\left(\mathcal{C}^{\sigma-d/p}(\mathcal{X},L)\right)$
$\displaystyle\asymp T^{1+1/p-\sigma/d}.$

Before proving this result, we make a few remarks:

######  Remark 5.

When $\sigma\leq d/p$, the picture is quite simple, as
$\mathcal{B}_{p,q}^{\sigma}\left(\mathcal{X},L\right)$ contains functions with
unbounded singularities, and it is easy to show that
$M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}\left(\mathcal{X},L\right)\right)=M_{C,T}\left(\mathcal{B}_{p,q}^{\sigma}\left(\mathcal{X},L\right)\right)=\infty$.
When $\sigma>d/p+1$, under cumulative regret, our lower bounds continue to
hold, and we conjecture that they are tight, while upper bounds are either
non-existent or loose. We leave it to future work to tighten these upper
bounds.

######  Remark 6.

Since the instantaneous regret at time $t$ of any algorithm for minimizing
cumulative regret will always be at least the instantaneous regret of the the
minimax-optimal procedure for minimizing simple regret at time $t$, one can
check that, for any $T$, $M_{C,T}\geq\sum_{t=1}^{T}M_{S,t}$. Also, given an
algorithm $\widehat{X}$ with expected cumulative regret $R$, one can construct
an algorithm with expected simple regret $\frac{R}{T-1}$ by sampling $X_{T}$
uniformly at random from the first $T-1$ queries of $\widehat{X}$. Thus, in
what follows, we prove lower bounds on simple regret and upper bounds on
cumulative regret, while lower bounds on cumulative regret and upper bounds on
simple regret follow via these inequalities.

We now turn to proving Theorem 4. We begin by constructing a “worst-case”
subset of Besov functions:

######  Lemma 7 (Construction of $\Theta_{j}$).

For any positive integer $j\in\mathbb{N}$, there exists a subset
$\Theta_{j}\subseteq\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)$ of Besov
functions with the following properties:

  1. (a)

$|\Theta_{j}|=2^{d\cdot j}$.

  2. (b)

For $f\neq g\in\Theta_{j}$, $f$ and $g$ have disjoint supports.

  3. (c)

For all $f\in\Theta_{j}$,
$\sup_{x\in\mathcal{X}}f(x)=L\phi^{*}2^{j(d/p-\sigma)}$, where
$\phi^{*}=\sup_{x\in\mathcal{X}}\phi(x)$ is a constant depending only on the
choice of father wavelet.

  4. (d)

For all $f\in\Theta_{j}$,

$\max_{f\in\Theta_{j}}\|f\|_{\infty}=L\|\phi\|_{\infty}2^{j(d/p-\sigma)}.$

###### Proof.

Let
$\Theta_{j}:=\\{L2^{-j(\sigma+d(1/2-1/p))}\phi_{j,\lambda}:\lambda\in[2^{j}]^{d}\\}$
be a rescaling of the wavelet basis functions at resolution $j$ by the
constant $L2^{-j(\sigma+d(1/2-1/p))}$. Properties (a) and (b) are immediate
from the construction of the wavelet basis. For property (c), note that, for
$f\in\Theta_{j}$,

$\displaystyle\max_{x\in\mathcal{X}}f(x)$
$\displaystyle=L\phi_{j,\lambda}^{*}2^{-j(\sigma+d(1/2-1/p))}$
$\displaystyle=L(\phi^{*}2^{dj/2})2^{-j(\sigma+d(1/2-1/p))}$
$\displaystyle=L\phi^{*}2^{j(d/p-\sigma)}.$

Checking property (d) is nearly identical, replacing $\phi^{*}$ with
$\|\phi\|_{\infty}$. Finally, to verify
$\Theta_{j}\subseteq\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)$, Eq. (1) gives,
for $\lambda\in[2^{j}]^{d}$:

$\displaystyle\|L2^{-j(\sigma+d(1/2-1/p))}\phi_{j,\lambda}\|_{\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})}$
$\displaystyle=2^{j(\sigma+d(1/2-1/p))}\langle
L2^{-j(\sigma+d(1/2-1/p))}\phi_{j,\lambda},\phi_{j,\lambda}\rangle=L.$

∎

Lemma 7 constructed a large family $\Theta_{j}$ of objective functions with
disjoint supports. We now use this lemma to prove the main Theorem 4. The
strategy is to show, based on the size of $\Theta_{j}$, that, for any
optimization strategy $\widehat{x}$, there is some objective function
$f\in\Theta_{j}$ such that none of the first $T$ points sampled by
$\widehat{x}$ lies in the support of $f$. Thus, when the true objective
function is $f$, the simple regret incurred by $\widehat{x}$ is
$\sup_{x\in\mathcal{X}}f(x)$.

###### Proof of Theorem 4.

For $s\in(0,1]$, existing results (e.g., Munos (2011); Malherbe and Vayatis
(2017)) imply the upper bound
$M_{C,T}\left(\mathcal{C}^{s}(\mathcal{X},L)\right)\lesssim T^{1-s/d}$ for the
Hölder case. Our upper bound on
$M_{C,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)$ follows from
the continuous embedding
$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})\hookrightarrow\mathcal{C}^{\sigma-d/p}(\mathcal{X})$
(Theorem 3). We now turn to lower bounds on expected simple regret.

Fix any optimization procedure $\widehat{x}$ and any random variable $Z$, and
let $\Theta=\Theta_{j}$ denote the set of functions constructed in Lemma 7
with $j=\left\lceil\frac{1}{d}\log_{2}(2T)\right\rceil$, so that
$2T\leq|\Theta|\leq 2^{d+1}T$ and

$\min_{f\in\Theta_{j}}\max_{x\in\mathcal{X}}f(x)\geq
L\phi^{*}2^{(d+1)(1/p-\sigma/d)}T^{1/p-\sigma/d}.$

For each $f\in\Theta$, let $E_{f}$ denote the event that
$f(X_{1})=\cdots=f(X_{T})=0$. Since the elements of $\Theta$ have disjoint
support, $\sum_{f\in\Theta}1_{E_{f}}\geq|\Theta|-T\geq T$, and so, since
$|\Theta|\leq 2^{d+1}T$,

$2^{d+1}T\cdot\max_{f\in\Theta}\Pr_{Z}[E_{f}]\geq\sum_{f\in\Theta}\Pr_{Z}[E_{f}]=\mathop{\mathbb{E}}_{Z}\left[\sum_{f\in\Theta}1_{E_{f}}\right]\geq
T.$

In particular, there exists a function444Note that $f_{\widehat{x},\mu_{Z}}$
depends only on the distribution $\mu_{Z}$ of $Z$, not on any particular
instance of $Z$. $f_{\widehat{x},\mu_{Z}}\in\Theta$ such that
$\Pr_{Z}[E_{f_{\widehat{x},\mu_{Z}}}]\geq 2^{-(d+1)}$. It follows that

$\displaystyle M_{T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)$
$\displaystyle=\inf_{\widehat{x},Z}\sup_{f\in\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)}\mathop{\mathbb{E}}_{Z}\left[\max_{x\in\mathcal{X}}f(x)-f(x_{T})\right]$
$\displaystyle\geq\inf_{\widehat{x},Z}\mathop{\mathbb{E}}_{Z}\left[\max_{x\in\mathcal{X}}f_{\widehat{x},\mu_{Z}}(x)-f_{\widehat{x},\mu_{Z}}(x_{T})\right]$
$\displaystyle\geq\inf_{\widehat{x},Z}\Pr_{Z}[E_{f_{\widehat{x},\mu_{Z}}}]\cdot\max_{x\in\mathcal{X}}f_{\widehat{x},\mu_{Z}}(x)$
$\displaystyle\geq 2^{-(d+1)}\cdot
L\phi^{*}2^{(d+1)(1/p-\sigma/d)}T^{1/p-\sigma/d}.$

∎

##  5 Noisy Case

In this section, we state and prove lower bounds on the minimax simple and
cumulative regrets in the case where the observed rewards are noisy. Our
approach will utilize the following version of Fano’s lemma:

######  Lemma 8.

(Fano’s Lemma; Simplified Form of Theorem 2.5 of Tsybakov (2008)) Fix a family
$\mathcal{P}$ of distributions over a sample space $\mathcal{X}$ and fix a
pseudo-metric $\rho:\mathcal{P}\times\mathcal{P}\to[0,\infty]$ over
$\mathcal{P}$. Suppose there exist $P_{0}\in\mathcal{P}$ and a set
$\Theta\subseteq\mathcal{P}$ such that

$\frac{1}{|\Theta|}\sum_{P\in\Theta}D_{KL}(P_{0},P)\leq\frac{\log|\Theta|}{16},$

where $D_{KL}:\mathcal{P}\times\mathcal{P}\to[0,\infty]$ denotes Kullback-
Leibler divergence. Then,

$\inf_{\widehat{P}}\sup_{P\in\mathcal{P}}\Pr\left[\rho(P,\widehat{P})\geq\frac{1}{2}\inf_{Q,R\in\Theta}\rho(Q,R)\right]\geq
1/8,$

where the first $\inf$ is taken over all estimators $\widehat{P}$.

######  Theorem 9 (Minimax Rates, Noisy Case).

For $\sigma>d/p$,

$\displaystyle
M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)$
$\displaystyle\asymp
M_{S,T}\left(\mathcal{C}^{\sigma-d/p}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)$
$\displaystyle\asymp\max\left\\{\left(\frac{\eta^{2}\log\frac{T}{\eta^{2}}}{T}\right)^{\frac{\sigma-d/p}{2(\sigma-d/p)+d}},M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)\right\\},$

and, for $\sigma\in(d/p,d/p+1]$,

$\displaystyle
M_{C,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)$
$\displaystyle\asymp
M_{C,T}\left(\mathcal{C}^{\sigma-d/p}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)$
$\displaystyle\asymp\max\Bigg{\\{}T^{\frac{d+\sigma-d/p}{2(\sigma-d/p)+d}}\left(\eta^{2}\log\frac{T}{\eta^{2}}\right)^{\frac{\sigma-d/p}{2(\sigma-d/p)+d}},$
$\displaystyle\hskip
42.67912ptM_{C,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)\Bigg{\\}}.$

As in the noiseless case, consistent optimization is possible only for
$\sigma>d/p$, and, while our lower bounds hold for $\sigma>d/p+1$, existing
upper bounds under cumulative regret appear loose (see Remark 5 and Related
Work in Section 3). Also, as noted in Remark 6, lower bounds on minimax
cumulative regret and upper bounds on minimax simple regret follow from lower
bounds on minimax simple regret and upper bounds on minimax cumulative regret,
respectively.

###### Proof.

For $s\in(0,1]$, Corollary 1 of Auer et al. (2007) implies the upper bound

$M_{C,T}\left(\mathcal{C}^{s}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)\lesssim
T^{\frac{s+d}{2s+d}}\left(\eta^{2}\log\frac{T}{\eta^{2}}\right)^{-\frac{s}{2s+d}}$

for the Hölder case. Our upper bound on
$M_{C,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L),\mathsf{subG}(\eta^{2})\right)$
follows from the continuous embedding
$\mathcal{B}_{p,q}^{\sigma}(\mathcal{X})\hookrightarrow\mathcal{C}^{\sigma-d/p}(\mathcal{X})$
(Theorem 3).

For proving a minimax lower bound over sub-Gaussian noise distributions with
variance parameter at most $\eta^{2}$, we may assume homoskedastic Gaussian
noise with variance $\eta^{2}$. Let $\Theta_{j}$ be as given in Lemma 7 with
$j$ specified in Eq. (5). Let $P_{0}$ denote the distribution of the
observations $(Y_{1},...,Y_{T})$ when the true function is constant $0$ on
$\mathcal{X}$, and, for each $f\in\Theta_{j}$, let $P_{f}$ denote the
distribution of $(Y_{1},...,Y_{T})$ when the true function is $f$. For
$f\in\Theta_{j}\cup\\{0\\}$ and each $t\in[T]$, let $P_{0,t}$ denote the
conditional distribution of $(Y_{i},...,Y_{T})|(Y_{1},...,Y_{t-1})$. By a
standard formula for KL divergence between two Gaussians and the construction
of $\Theta_{j}$ (specifically part (d) of Lemma 7), the maximum information
gain $D_{j}^{*}$ of a single query is:

$\displaystyle D_{j}^{*}$
$\displaystyle:=\sup_{x\in\mathcal{X},t\in[T],f\in\Theta_{j}}D_{KL}\left(P_{0,t}\|P_{f,t}\right)$
$\displaystyle=\max_{f\in\Theta_{j}}\frac{\|f\|_{\infty}^{2}}{2\eta^{2}}=\frac{L^{2}\|\phi\|_{\infty}^{2}2^{2j(d/p-\sigma)}}{2\eta^{2}}$
(4)

Then, by the chain rule for KL divergence, for any
$f\in\Theta_{j}\cup\\{0\\}$,

$\displaystyle D_{\text{KL}}(P_{0}\|P_{f})$
$\displaystyle=\sum_{t=1}^{T}\mathop{\mathbb{E}}_{P_{0}}\left[D(P_{0,t}\|P_{f,t})\middle|Y_{1},...,Y_{t-1}\right]$
$\displaystyle=\sum_{t=1}^{T}\mathop{\mathbb{E}}_{P_{0}}\left[D(P_{0,t}\|P_{f,t})\middle|Y_{1},...,Y_{t-1},f(X_{t})>0\right]$
$\displaystyle\hskip 30.0pt\cdot\Pr_{P_{0}}[f(X_{t})>0|Y_{1},...,Y_{t-1}]$
$\displaystyle\leq
D_{j}^{*}\sum_{t=1}^{T}\Pr_{P_{0}}[f(X_{i})>0|Y_{1},...,Y_{t-1}].$

Since the elements of $\Theta_{j}$ have disjoint support,

$\displaystyle\sum_{f\in\Theta_{j}}\Pr_{P_{0}}[f(X_{T})>0|Y_{1},...,Y_{T}]$
$\displaystyle=\mathop{\mathbb{E}}_{P_{0}}\left[\sum_{f\in\Theta_{j}}1_{\\{f(X_{T})>0\\}}|Y_{1},...,Y_{T}\right]\leq
1,$

and so, since $|\Theta_{j}|=2^{dj}$ (by part (a) of Lemma 7),

$\displaystyle\frac{1}{|\Theta_{j}|}\sum_{f\in\Theta_{j}}D_{\text{KL}}(P_{0}\|P_{f})$
$\displaystyle\leq\frac{1}{|\Theta_{j}|}\sum_{f\in\Theta_{j}}D_{j}^{*}\sum_{i=1}^{T}\Pr_{0}[f(X_{i})>0|Y_{1},...,Y_{i-1}]$
$\displaystyle\leq\frac{D_{j}^{*}T}{|\Theta_{j}|}=\frac{L^{2}\|\phi\|_{\infty}^{2}2^{2j(d/p-\sigma)}T}{2\eta^{2}2^{dj}}.$

Defining $r:=d+2(\sigma-d/p)>0$ and
$z:=\frac{L^{2}\|\phi\|_{\infty}^{2}}{2}\frac{T}{\eta^{2}}$,

$\frac{1}{|\Theta_{j}|}\sum_{f\in\Theta_{j}}D_{\text{KL}}(P_{0}\|P_{f})\leq
z2^{-rj}.$

Let

$j=\frac{1}{r}\log_{2}\frac{8rz}{d\log_{2}8rz}.$ (5)

Since $z\to\infty$ as $\frac{T}{\eta^{2}}\to\infty$, for sufficiently large
$\frac{T}{\eta^{2}}$, we have $\log_{2}(8rz)\geq 2d\log_{2}\log_{2}(8rz)$.
Therefore,

$\displaystyle\frac{1}{|\Theta_{j}|}\sum_{f\in\Theta_{j}}D_{\text{KL}}(P_{0}\|P_{f})\leq
z2^{-rj}$ $\displaystyle=\frac{d}{8r}\log_{2}(8rz)$
$\displaystyle\leq\frac{d}{16r}\left(\log_{2}(8rz)-d\log_{2}\log_{2}(8rz)\right)$
$\displaystyle=\frac{d}{16r}\log_{2}\frac{8rz}{d\log_{2}(8rz)}=\frac{dj}{16}=\frac{\log_{2}|\Theta_{j}|}{16}.$

Therefore, we can apply Fano’s Lemma (Lemma 8) to the class
$\\{P_{f}:f\in\Theta_{j}\\}$ equipped with the discrete metric
$\rho(P,Q)=1_{\\{P\neq Q\\}}$, giving

$\inf_{\widehat{x},Z}\sup_{f\in\Theta_{j}}\Pr\left[\rho(P_{f},P_{\widehat{x},Z})\geq\frac{1}{2}\inf_{f,g\in\Theta_{j}}\rho(P_{f},P_{g})\right]\geq\frac{1}{8},$

where $P_{\widehat{x},Z}$ is any distribution “estimate” (i.e., any
distribution computed as a function of the data
$\\{(X_{t},Y_{t})\\}_{t\in[T]}$ and $Z$). In particular, since the elements of
$\Theta_{j}$ have disjoint support, we can construct a distribution estimate
$P_{\widehat{x},Z}:=P_{f}$, where $f$ is the unique element of $\Theta_{j}$
with $f(X_{T})>0$ (if $f(X_{T})=0$ for all $f\in\Theta_{j}$, then $f$ can be
selected at random). For this particular estimate, we have

$\displaystyle\inf_{\widehat{x},Z}\sup_{f\in\Theta_{j}}\Pr\left[P_{f}\neq
P_{\widehat{x},Z}\right]$
$\displaystyle=\inf_{\widehat{x},Z}\sup_{f\in\Theta_{j}}\Pr\left[\rho(P_{f},P_{\widehat{x},Z})\geq\frac{1}{2}\inf_{f,g\in\Theta_{j}}\rho(P_{f},P_{g})\right]\geq\frac{1}{8}.$

Therefore, since
$\Theta_{j}\subseteq\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)$,

$\displaystyle
M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)\geq
M_{S,T}\left(\Theta_{j}\right)$
$\displaystyle=\inf_{\widehat{x},Z}\sup_{f\in\Theta_{j}}\mathop{\mathbb{E}}_{f}\left[\sup_{x\in\mathcal{X}}f(x)-f(X_{n})\right]$
$\displaystyle\geq\inf_{\widehat{x},Z}\sup_{f\in\Theta_{j}}\mathop{\mathbb{E}}_{f}\left[L\phi^{*}2^{j(d/p-\sigma)}1_{\\{f(X_{n})=0\\}}\right]$
$\displaystyle=L\phi^{*}2^{j(d/p-\sigma)}\inf_{\widehat{x},Z}\sup_{f\in\Theta_{j}}\Pr\left[P_{f}\neq
P_{\widehat{x},Z}\right]$
$\displaystyle\geq\frac{L\phi^{*}2^{j(d/p-\sigma)}}{8}$

Plugging in $j$ from Equation 5 gives the desired result:

$\displaystyle M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)$
$\displaystyle\qquad\geq\frac{L\phi^{*}}{8}\left(\frac{4rL^{2}\|\phi\|_{\infty}^{2}\frac{T}{\eta^{2}}}{d\log_{2}\left(4rL^{2}\|\phi\|_{\infty}^{2}\frac{T}{\eta^{2}}\right)}\right)^{\frac{d/p-\sigma}{2(\sigma-d/p)+d}},$

where $r=d+2\sigma-2d/p>0$ is a constant. ∎

To illustrate results of Theorem 9, Figure 2 shows a phase diagram of minimax
simple regret rates, as a function of parameters $\sigma$ and $p$ of the Besov
space in which the objective function is assumed to lie. The figure highlights
the main conclusion of our paper: minimax rates of bandit optimization over a
Besov space correspond precisely to rates over the smallest Hölder space
within which that Besov space embeds.

Comparing with the noiseless simple regret rate $T^{1/p-\sigma/d}$ from
Theorem 4, up to logarithmic factors, the noiseless rate dominates for
$\eta\lesssim T^{1/p-\sigma/d}$, while the noisy rate dominates for
$\eta\gtrsim T^{1/p-\sigma/d}$. This threshold matches the convergence rate in
the noiseless case; a natural interpretation is that the optimal error
decreases at the fast noiseless rate $T^{1/p-\sigma/d}$ until it is of the
same order as the noise $\eta$, at which point convergence slows to the noisy
rate
$T^{\frac{d+\sigma-d/p}{2(\sigma-d/p)+d}}\eta^{\frac{2\sigma-2d/p}{2(\sigma-d/p)+d}}$.

Figure 2: Phase diagram showing exponent of minimax simple regret rate in the
noisy case (i.e., $\alpha$ such that
$M_{S,T}\left(\mathcal{B}_{p,q}^{\sigma}(\mathcal{X},L)\right)\asymp
T^{-\alpha}$), plotted as a function of the Besov space parameters $\sigma$
and $1/p$, in dimension $d=1$. Dashed lines indicate Hölder ($p=\infty$) and
(Hilbert-)Sobolev ($p=2$) special cases studied in prior work. The contours,
corresponding to level sets of $\sigma-d/p$, demonstrate the relationship
between Besov spaces and the embedding Hölder space.

##  6 Conclusion

In this paper, we derived novel lower bounds on minimax simple and cumulative
regrets over general Besov spaces, in both the presence and absence of noise.
Via continuous embedding of Besov spaces in Hölder spaces, these rates match
existing upper bounds, implying optimality of existing proposed algorithms and
identifying optimal rates over all Besov spaces.

Our results suggest that, up to constant factors, $0^{th}$-order methods
cannot benefit from the additional $d/p$ orders of smoothness in
$\mathcal{B}_{p,q}^{s+d/p}(\mathcal{X})$ over $\mathcal{C}^{s}(\mathcal{X})$.
From a theoretical perspective, this suggests that Hölder spaces are “natual”
spaces over which to study continuum-armed bandit optimization, since, optimal
algorithms over Hölder spaces are automatically optimal over other Besov
spaces, while the converse may not hold. From a practical perspective, our
results suggest there may be little (at most constant-factor) benefit to
employing spatially adaptive methods (i.e., those that more finely represent
areas of the search space where the objective is less smooth), contrasting
from results in regression and other nonparametric estimation problems, where
spatial adaptivity is needed to obtain minimax rates (Donoho and Johnstone,
1998).

Finally, to the best of our knowledge, the results of this paper are the first
to consider modeling continuum-armed bandits using wavelets. Future work
should explore efficient bandit algorithms based on fast wavelet methods,
which can provide significant computational advantages over kernel methods
(Mallat, 2008).

#### Acknowledgements

The author would like to thank Ananya Uppal and Arkady Epshteyn for helpful
insights and discussions, and Sivaraman Balakrishnan for pointing out an error
in an earlier version of the paper.

## References

  * Agrawal (1995) Rajeev Agrawal.  The continuum-armed bandit problem.  _SIAM Journal on Control and Optimization_ , 33(6):1926–1951, 1995.
  * Akhavan et al. (2020) Arya Akhavan, Massimiliano Pontil, and Alexandre B Tsybakov.  Exploiting higher order smoothness in derivative-free optimization and continuous bandits.  _arXiv preprint arXiv:2006.07862_ , 2020.
  * Auer et al. (2007) Peter Auer, Ronald Ortner, and Csaba Szepesvári.  Improved rates for the stochastic continuum-armed bandit problem.  In _International Conference on Computational Learning Theory_. Springer, 2007.
  * Bubeck et al. (2011a) Sébastien Bubeck, Rémi Munos, Gilles Stoltz, and Csaba Szepesvári.  X-armed bandits.  _Journal of Machine Learning Research_ , 12(May), 2011a.
  * Bubeck et al. (2011b) Sébastien Bubeck, Gilles Stoltz, and Jia Yuan Yu.  Lipschitz bandits without the Lipschitz constant.  In _International Conference on Algorithmic Learning Theory_. Springer, 2011b.
  * Bubeck et al. (2018) Sébastien Bubeck, Michael Cohen, and Yuanzhi Li.  Sparsity, variance and curvature in multi-armed bandits.  In _Algorithmic Learning Theory_ , 2018.
  * Bull (2011) Adam D Bull.  Convergence rates of efficient global optimization algorithms.  _Journal of Machine Learning Research_ , 12(Oct):2879–2904, 2011.
  * Cope (2009) Eric Cope.  Regret and convergence bounds for immediate-reward reinforcement learning with continuous action spaces.  _IEEE Transactions on Automatic Control_ , 54(6):1243–1253, 2009.
  * Daubechies (1992) Ingrid Daubechies.  _Ten lectures on wavelets_ , volume 61.  SIAM, 1992.
  * Donoho and Johnstone (1998) David L Donoho and Iain M Johnstone.  Minimax estimation via wavelet shrinkage.  _The Annals of Statistics_ , 26(3):879–921, 1998\.
  * Donoho et al. (1995) David L Donoho, Iain M Johnstone, Gérard Kerkyacharian, and Dominique Picard.  Wavelet shrinkage: asymptopia?  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 1995.
  * Fan and Gijbels (1996) Jianqing Fan and Irene Gijbels.  _Local polynomial modelling and its applications: monographs on statistics and applied probability_.  CRC Press, 1996.
  * Feurer and Hutter (2019) Matthias Feurer and Frank Hutter.  Hyperparameter optimization.  In _Automated Machine Learning_ , pages 3–33. Springer, 2019.
  * Golovin et al. (2017) Daniel Golovin, Benjamin Solnik, Subhodeep Moitra, Greg Kochanski, John Karro, and D Sculley.  Google Vizier: A service for black-box optimization.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1487–1495, 2017.
  * Golovin et al. (2019) Daniel Golovin, John Karro, Greg Kochanski, Chansoo Lee, Xingyou Song, and Qiuyi Zhang.  Gradientless descent: High-dimensional zeroth-order optimization.  _arXiv preprint arXiv:1911.06317_ , 2019.
  * Grant and Leslie (2020) James A Grant and David S Leslie.  On Thompson sampling for smoother-than-Lipschitz bandits.  _arXiv preprint arXiv:2001.02323_ , 2020.
  * Grill et al. (2015) Jean-Bastien Grill, Michal Valko, and Rémi Munos.  Black-box optimization of noisy functions with unknown smoothness.  In _Advances in Neural Information Processing Systems_ , pages 667–675, 2015.
  * Haroske and Schneider (2009) Dorothee D Haroske and Cornelia Schneider.  Besov spaces with positive smoothness on $\mathbb{R}^{n}$, embeddings and growth envelopes.  _Journal of Approximation Theory_ , 161(2):723–747, 2009.
  * Kanagawa et al. (2018) Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur.  Gaussian processes and kernel methods: A review on connections and equivalences.  _arXiv preprint arXiv:1807.02582_ , 2018.
  * Kandasamy (2018) Kirthevasan Kandasamy.  _Tuning Hyperparameters without Grad Students: Scaling up Bandit Optimisation_.  PhD thesis, Carnegie Mellon University, 2018.
  * Kiefer and Wolfowitz (1952) Jack Kiefer and Jacob Wolfowitz.  Stochastic estimation of the maximum of a regression function.  _The Annals of Mathematical Statistics_ , 1952.
  * Kleinberg et al. (2008) Robert Kleinberg, Aleksandrs Slivkins, and Eli Upfal.  Multi-armed bandits in metric spaces.  In _Proceedings of the fortieth annual ACM symposium on Theory of computing_ , pages 681–690, 2008.
  * Kleinberg (2005) Robert D Kleinberg.  Nearly tight bounds for the continuum-armed bandit problem.  In _Advances in Neural Information Processing Systems_ , 2005.
  * Kochanski et al. (2017) Greg Kochanski, Daniel Golovin, John Karro, Benjamin Solnik, Subhodeep Moitra, and D Sculley.  Bayesian optimization for a better dessert.  In _Proceedings of NIPS Workshop on Bayesian Optimization_ , 2017\.
  * Larson et al. (2019) Jeffrey Larson, Matt Menickelly, and Stefan M Wild.  Derivative-free optimization methods.  _Acta Numerica_ , 28:287–404, 2019.
  * Leoni (2017) Giovanni Leoni.  _A first course in Sobolev spaces_.  American Mathematical Soc., 2017.
  * Li et al. (2017) Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  _The Journal of Machine Learning Research_ , 18(1):6765–6816, 2017.
  * Liu et al. (2020) Yusha Liu, Yining Wang, and Aarti Singh.  Smooth bandit optimization: Generalization to hölder space.  _arXiv preprint arXiv:2012.06076_ , 2020.
  * Malherbe and Vayatis (2017) Cedric Malherbe and Nicolas Vayatis.  Global optimization of lipschitz functions.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , pages 2314–2323, 2017.
  * Mallat (2008) Stéphane Mallat.  _A wavelet tour of signal processing_.  Elsevier, 2008.
  * Meyer (1992) Yves Meyer.  _Wavelets and operators_ , volume 1.  Cambridge university press, 1992.
  * Minsker (2012) Stanislav Minsker.  _Non-asymptotic bounds for prediction problems and density estimation._ PhD thesis, Georgia Institute of Technology, 2012.
  * Minsker (2013) Stanislav Minsker.  Estimation of extreme values and associated level sets of a regression function via selective sampling.  In _Conference on Learning Theory_ , pages 105–121, 2013.
  * Munos (2011) Rémi Munos.  Optimistic optimization of a deterministic function without the knowledge of its smoothness.  In _Advances in neural information processing systems_ , pages 783–791, 2011.
  * Nakamura et al. (2017) Nathan Nakamura, Jason Seepaul, Joseph B Kadane, and B Reeja-Jayan.  Design for low-temperature microwave-assisted crystallization of ceramic thin films.  _Applied Stochastic Models in Business and Industry_ , 33(3):314–321, 2017.
  * Nemirovski (1985) Arkadi S Nemirovski.  Nonparametric estimation of smooth regression functions.  _Izv. Akad. Nauk. SSR Teckhn. Kibernet_ , 3:50–60, 1985\.
  * Pandey et al. (2007) Sandeep Pandey, Deepak Agarwal, Deepayan Chakrabarti, and Vanja Josifovski.  Bandits for taxonomies: A model-based approach.  In _Proceedings of the 2007 SIAM International Conference on Data Mining_ , pages 216–227. SIAM, 2007.
  * Sawano (2018) Yoshihiro Sawano.  _Theory of Besov spaces_ , volume 56.  Springer, 2018.
  * Scarlett et al. (2017) Jonathan Scarlett, Ilija Bogunovic, and Volkan Cevher.  Lower bounds on regret for noisy gaussian process bandit optimization.  In _Conference on Learning Theory_ , pages 1723–1742, 2017.
  * Slivkins (2019) Aleksandrs Slivkins.  Introduction to multi-armed bandits.  _Foundations and Trends in Machine Learning_ , 12(1-2):1–286, 2019.
  * Srinivas et al. (2010) Niranjan Srinivas, Andreas Krause, Sham Kakade, and Matthias Seeger.  Gaussian process optimization in the bandit setting: no regret and experimental design.  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , pages 1015–1022, 2010.
  * Thompson (1933) William R Thompson.  On the likelihood that one unknown probability exceeds another in view of the evidence of two samples.  _Biometrika_ , 1933.
  * Tsybakov (2008) Alexandre B Tsybakov.  _Introduction to nonparametric estimation_.  Springer Science & Business Media, 2008.
  * Wang et al. (2018a) Yining Wang, Sivaraman Balakrishnan, and Aarti Singh.  Optimization of smooth functions with noisy observations: Local minimax rates.  In _Advances in Neural Information Processing Systems_ , pages 4338–4349, 2018a.
  * Wang et al. (2018b) Yining Wang, Simon Du, Sivaraman Balakrishnan, and Aarti Singh.  Stochastic zeroth-order optimization in high dimensions.  In _International Conference on Artificial Intelligence and Statistics_ , 2018b.
  * Wasserman (2013) Larry Wasserman.  _All of statistics: a concise course in statistical inference_.  Springer Science & Business Media, 2013.