[1]Nicoletta D’Angelo

[1]Department of Economics, Business and Statistics, University of Palermo,
Palermo, Italy

2]Sustainable Mobility Center (Centro Nazionale per la Mobilità
Sostenibile—CNMS)

3]MOX - Department of Mathematics, Politecnico di Milano, Milano, Italy

# Spatio-temporal point process modelling of fires in Sicily exploring human
and environmental factors

<EMAIL_ADDRESS>Alessandro Albano<EMAIL_ADDRESS>Andrea Gilardi<EMAIL_ADDRESS>Giada Adelfio<EMAIL_ADDRESS>*  [  [

###### Abstract

In 2023, Sicily faced an escalating issue of uncontrolled fires, necessitating
a thorough investigation into their spatio-temporal dynamics. Our study
addresses this concern through point process theory. Each wildfire is treated
as a unique point in both space and time, allowing us to assess the influence
of environmental and anthropogenic factors by fitting a spatio-temporal
separable Poisson point process model, with a particular focus on the role of
land usage. First, a spatial log-linear Poisson model is applied to
investigate the influence of land use types on wildfire distribution,
controlling for other environmental covariates. The results highlight the
significant effect of human activities, altitude, and slope on spatial fire
occurrence. Then, a Generalized Additive Model with Poisson-distributed
response further explores the temporal dynamics of wildfire occurrences,
confirming their dependence on various environmental variables, including the
maximum daily temperature, wind speed, surface pressure, and total
precipitation.

###### keywords:

Fires, Land usage, Point processes, Spatial analysis, Intensity estimation.

##  1 Introduction

In recent years, Sicily has experienced an alarming increase in uncontrolled
wildfires, establishing it as the Italian region with the highest frequency of
fire events and, consequently, the largest burned area. The occurrence of
these fires is primarily associated with ignition sources, forest fuels, and
environmental conditions [23, 25, 34]. Ignition sources are typically
categorized into natural causes, such as lightning and geological factors, and
human causes, both accidental and intentional [3, 35]. Human-driven factors,
particularly arson (but also, unintentional wildfires), often ignited for
purposes like creating new pasture resources or burning stubble, stand out as
the primary contributors to wildfires in Sicily, especially in areas where
vegetation interfaces with urban structures [20].

Wildfires pose a significant threat to ecosystems, human settlements, and
economic activities, demanding a comprehensive understanding of their spatio-
temporal dynamics for effective mitigation and management. The island of
Sicily, nestled in the heart of the Mediterranean Sea, has a rich history and
diverse landscapes, making it vulnerable to the escalating impacts of climate
change, including an increased frequency and intensity of wildfires. In the
year 2023, Sicily faced a pronounced wildfire season, underscoring the urgent
need for advanced analytical methodologies to unravel the underlying patterns
and drivers of these destructive events.

Traditional approaches to wildfire analysis often rely on descriptive
statistics and basic spatial visualization, providing limited insights into
the effect of environmental, climatic, and anthropogenic factors that
contribute on wildfire occurrences. This paper advocates for the application
of point process methodology to explore the spatial and temporal dynamics of
Sicilian wildfires in 2023.

Point process methodology, rooted in statistical theory, enables the modeling
of events occurring in space and time, making it an ideal tool for in-depth
analysis of wildfire patterns. By treating each wildfire occurrence as a point
in a spatial and temporal domain, we aim to discern potential predictors that
may contribute to the ignition of wildfires across the Sicilian landscape.

Moreover, traditional analyses often fall short of capturing the relationships
between environmental, climatic, and anthropogenic factors that contribute to
wildfires. By exploiting spatio-temporal models, we aim to assess which
variables influence the probability of fire occurrence. Including covariates
such as land usage, altitude, slope, temperature, precipitation, surface
pressure, and wind speed provides a comprehensive framework to understand the
underlying mechanisms influencing wildfire patterns.

As noted by [11], fires and land usage are intrinsically connected, but
research investigating this dynamic is limited. Specifically, the influence of
land use on fire ignition is a well-documented phenomenon. Previous papers
have shown that land use, particularly the presence of roads, can
significantly impact the ignition of fires [33], with the influence of roads
being much stronger in less dense land cover areas. As human activities
increasingly encroach upon natural landscapes, understanding the impact of
land usage becomes crucial for effective wildfire management. Through the
proposed modelling, we seek to quantify the relationship between different
land use types and the probability of fire occurrence, finding areas where
human activities may be making the Sicilian landscape worse.

The proposed spatio-temporal separable Poisson point process model seeks to
determine the specific contributions of land usage categories on fire
occurrences while controlling for other environmental covariates. As a matter
of fact, understanding the influence of land usage on wildfire probability is
crucial for designing targeted interventions, land-use planning, and
developing resilient ecosystems in the face of evolving environmental
challenges. We first fit a spatial log-linear Poisson model in order to
investigate the influence of land use types on wildfire occurrence,
controlling for other environmental covariates. The results highlight the
significant effect of human activities, altitude, and slope on spatial fire
occurrence. Secondly, a Generalized Additive Model with Poisson-distributed
response is fitted to model the temporal fire occurrences, confirming their
dependence on various environmental variables, including the maximum daily
temperature, wind speed, surface pressure, and total precipitation. This study
seeks to contribute not only to the scientific understanding of wildfire
dynamics but also to the practical tools necessary for mitigating their
impacts on communities and ecosystems.

All the analyses are carried out through the R statistical software [32],
making particular use of the packages stopp [14], sf [30], and stars [31].

The structure of the paper is as follows. Section 2 presents the data. Section
3 provides an overview of spatio-temporal point processes and the model
employed in the paper, while Section 4 illustrates the actual model fitting
through the spatial and temporal intensity estimation. Section 5 is devoted to
the discussion and conclusions.

##  2 Data

This section focuses on the description of the data employed in the paper,
namely the fire data, which will represent our spatio-temporal point pattern
under analysis, and some spatial and/or temporal variables which will serve as
covariates, assumed to influence the overall occurrence of points, as typical
in point process analysis.

###  2.1 Fires point pattern

To conduct our analysis, we rely on data obtained from the Fire Information
for Resource Management System (FIRMS) platform, accessible for download at
the following URL: https://firms.modaps.eosdis.nasa.gov/download/. This source
provides near real-time active fire locations to natural resource managers.

The variables considered in our study are:

  * •

Latitude: Center of 1 km fire pixel, but not necessarily the actual location
of the fire as one or more fires can be detected within the 1 km pixel.

  * •

Longitude: Center of 1 km fire pixel, but not necessarily the actual location
of the fire as one or more fires can be detected within the 1 km pixel.

  * •

Acq_Date (Acquisition Date): Date of MODIS acquisition.

  * •

Acq_Time (Acquisition Time): Time of acquisition/overpass of the satellite (in
UTC).

In particular, Longitude and Latitude will be used as spatial coordinates of
our point pattern, while Acq_Date and Acq_Time will be combined in a unique
variable serving as the time occurrence of the fire, with the smallest detail
as the hour of occurrence within a day.

Focusing exclusively on the spatial distribution, Figure 1 presents fire
counts in Italy throughout the year 2023.

Figure 1: Choropleth map displaying the counts of the wildfires that occurred
in each Italian region during 2023.

In particular, it shows the counts by region of the 26,724 fires that occurred
in Italy in 2023. Notably, the graphical representation emphasizes a
pronounced concentration of fires in the Sicilian territory, with the second-
highest incidence observed in the Southern region of Puglia (Apulia).
Following this, Figure 2 compares fire counts that occurred in Italy and
Sicily in 2023.

Figure 2: Barplot comparing the number of fires that occurred in the whole
country (yellow bars) and just in the Sicily region (red bars).

It is evident that the summer months, from July to October, stand out as
critical periods with the highest number of fires. Specifically, July emerges
as the most challenging month, recording a total of 8,842 wildfires in Italy,
of which 3,814 occurred exclusively in Sicily. Finally, Figure 3 shows the
spatial distribution of Sicilian fires in 2023. The inset map (which is
reported in the top-left) displays the whole observation window, which
includes the mainland of the region and three smaller islands (namely
Pantelleria, Lampedusa, and Linosa, respectively), located in the far South
and reported in the bottom-left as inset maps. The figure notably indicates
clustering behaviour, particularly in the western parts of the mainland.

Figure 3: Spatial distribution of fires recorded in the Sicily region during
2023. The inset map reported in the top-left displays the whole observation
window which includes the mainland and three smaller islands (namely
Pantelleria, Lampedusa, and Linosa, respectively), located in the far South
and reported in the bottom-left as inset maps.

###  2.2 Land use

The main focus of our work is the analysis of the fire occurrences depending
on the land usage since, as noted later, this type of information can be used
as a spatial covariate for our purposes. The land use data come from the
CORINE Land Cover (CLC) dataset111URL: https://land.copernicus.eu/en/cart-
downloads; Data downloaded during 2023-11, which is a comprehensive land cover
and land-use database created by the European Environment Agency (EEA) to
facilitate environmental monitoring and assessment across Europe. It is
usually employed to provide insights into the changing landscapes of Europe,
aiding in several projects such as environmental management and spatial
planning, in addition to wildfire analysis as exemplified in this paper.
Combining satellite imagery (obtained from Sentinel-2 and Landsat-8 projects)
and ground-based information, CLC classifies land cover into a hierarchical
raster system of over 40 classes, including urban areas, forests, water
bodies, and agricultural lands. The detailed classification according to CLC
comes in Table 1, whereas the product documentation and the nomenclature
guidelines can be browsed on the project’s website222URL:
https://land.copernicus.eu/content/corine-land-cover-nomenclature-
guidelines/html/. Last access: Feb. 2024..

Table 1: Land Use Classification according to Corine Land Cover Legend. Level 1 | Level 2 | Level 3
---|---|---
Artificial surfaces | Urban fabric | Continuous urban fabric
|  | Discontinuous urban fabric
| Industrial, commercial | Industrial or commercial units
| and transport units | Road and rail networks and associated land
|  | Port areas
|  | Airports
| Mine, dump | Mineral extraction sites
| and construction sites | Dump sites
|  | Construction sites
| Artificial, non-agricultural | Green urban areas
| vegetated areas | Sport and leisure facilities
Agricultural areas | Arable Land | Non-irrigated arable land
|  | Permanently irrigated land
| Permanent crops | Vineyards
|  | Fruit trees and berry plantations
|  | Olive groves
| Heterogeneous | Annual crops associated with permanent crops
| agricultural areas | Complex cultivation patterns
|  | Land principally occupied by agriculture
|  | with significant areas of natural vegetation
Forest and semi-natural areas | Forests | Broad-leaved forest
|  | Coniferous forest
|  | Mixed forest
| Scrub and/or herbaceous | Natural grasslands
| vegetation associations | Moors and heathland
|  | Sclerophyllous vegetation
|  | Transitional woodland-shrub
| Open spaces with little | Beaches dunes sands
| or no vegetation | Bare rocks
|  | Sparsely vegetated areas
|  | Burnt areas
Wetlands | Inland wetlands | Inland marshes
| Marine wetlands | Salt marshes
|  | Salines
Water bodies | Inland waters | Water courses
|  | Water bodies
| Marine waters | Coastal lagoons
|  | Sea and ocean

Figure 4 depicts the downloaded data on land usage categorized by level 1 of
the CLC Legend after merging the 4th and 5th categories, i.e. wetlands and
water bodies (“waterbodies” from now on). In particular, we represent the
macro land usage classification, whose categories are:

  * •

Artificial surfaces

  * •

Agricultural areas

  * •

Forest and semi-natural areas

  * •

Water bodies

Figure 4: Land usage in Sicily as reported in the CORINE Land Cover dataset.
The inset map included in the top-left zooms around the area of Palermo (which
in the capital of the region), whereas the three inset maps displayed in the
bottom-left showcase the land usage for Pantelleria, Lampedusa, and Linosa.

As shown in Table 2, the majority of the Sicilian territory is constituted by
agricultural areas ($68\%$, mainly in the souther parts of the region),
followed by forests and semi-natural areas ($26\%$, typically located in
proximity of some mountain chains in the North and North-East areas). Then,
artificial surfaces and water bodies represent the smallest portion ($5\%$ and
$1\%$, respectively).

Table 2: Area percentage, number of fires, and number of fires over the area percentage of the macro land usage. | Area perc. | Points | Points perc.
---|---|---|---
Artificial surfaces | 0.050 | 858 | 0.081
Agricultural areas | 0.683 | 6529 | 0.613
Forest and semi-natural areas | 0.259 | 3224 | 0.303
Water bodies | 0.007 | 45 | 0.004

###  2.3 Environmental covariates

The other covariates considered in this work are now introduced, detailing the
sources from where we downloaded the raw data and explaining the pre-
processing operations adopted to convert them into an usable format.

First, we downloaded the Digital Elevation Model (DEM) data for the Sicily
region from the webpage333URL:
https://tinitaly.pi.ingv.it/Download_Area1_1.html. Last access: Nov. 2023. of
the National Institute of Geophysics and Vulcanology [37]. Starting from the
DEM of the whole country, we selected the tiles belonging to the Sicily region
and combined them to obtain a unique raster object with 10m resolution that
represents the altitude in the area of interest. Then, we used the Horn’s
formula [26, 38] through the GDAL DEM utility command-line tool [24] to derive
the slope from the DEM data. More precisely, the Horn’s formula says that the
slope $\phi$ can be derived as

$\phi=\text{arctan}\left(\sqrt{\frac{\partial\text{Alt}}{\partial\text{lon}}+\frac{\partial\text{Alt}}{\partial\text{lat}}}\right)$
(1)

where Alt denotes the altitude dimension and
$\partial\text{Alt}/\partial\text{lon}$ and
$\partial\text{Alt}/\partial\text{lat}$ represent variations in the East-West
and North-South axis, respectively. The two partial derivatives in Equation
(1) can be approximated by using finite-order methods such as

$\frac{\partial\text{Alt}}{\partial\text{lon}}=\frac{(\text{Alt}_{1}+2\text{Alt}_{2}+\text{Alt}_{3})-(\text{Alt}_{7}-2\text{Alt}_{6}+\text{Alt}_{5})}{8\Delta\text{lon}}$
(2)

where $\text{Alt}_{1},\dots,\text{Alt}_{7}$ denote the altitude values in a
$3\times 3$ grid surrounding the cell under analysis (i.e. the cell in the
middle of the grid), following the scheme detailed in Figure 5 and
$\Delta\text{lon}$ represents the spacing between points in the horizontal
direction. Similar considerations hold for the estimation of
$\partial\text{Alt}/\partial\text{lat}$ and we refer to Weih Jr and Mattson,
[38] and the GDAL DEM webpage444URL: https://gdal.org/programs/gdaldem.html
for more details.

Figure 5: Grid map displaying the relationship between the raster cells used
to approximate the variations in the axis directions as used in Equation (2).

The altitude and slope covariates are displayed in the top panel of Figure 6.
As we can see, the first map clearly highlights the Etna volcano (located in
the eastern part of the region), the mountain chains in the North and North-
East areas (Monti Nebrodi), and the Catania valley surrounding the homonymous
city. The slope surface is generally flat or almost flat, apart from some
regions in the North-East (near Messina).

Figure 6: Top panels: Spatial covariates representing the elevation and slope
in Sicily: $Z_{elev}(u)$ and $Z_{slope}(u)$; Bottom panels: Daily means of the
surface pressure and skin temperature in Sicily in 2023.

Finally, Copernicus ERA5 data contains many environmental covariates that were
downloaded from https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-
era5-land?tab=overview. All of them are obtained in a $38\times 34$ spatial
grid with $10\text{km}\times 10\text{km}$ raster cells and hourly resolution
for each day of the year. The values are averaged daily over the whole year to
match the temporal resolution of fires data. The bottom panel of Figure 6
contains an example of the spatio-temporal resolution of such variables, and
the complete description of the available covariates comes as follows.

The 10u and 10v variables are the eastward and northward components,
respectively, of the 10m wind. They represent the horizontal speed of air
moving towards the East and North, at a height of ten metres above the surface
of the Earth, in metres per second. In particular, the U wind component is
parallel to the x-axis (i.e. longitude) while the V wind component is parallel
to the y-axis (i.e. latitude). In particular, a positive U wind comes from the
West, and a positive V wind comes from the South.

The surface pressure sp is the pressure (force per unit area) of the
atmosphere on the surface of land, sea and in-land water. It is a measure of
the weight of all the air in a column vertically above the area of the Earth’s
surface represented at a fixed point. The units of this variable are Pascals
(Pa).

The total precipitation tp is the accumulated liquid and frozen water,
comprising rain and snow, that falls to the Earth’s surface, not including
fog, dew or the precipitation that evaporates in the atmosphere before it
lands at the surface of the Earth. This variable is the total amount of water
accumulated over a particular time period, which depends on the data
extracted. The units of this variable are depth in metres of water equivalent.

The dew point temperature d2m is the temperature to which the air, at 2 metres
above the surface of the Earth, would have to be cooled for saturation to
occur. The air temperature at 2m above the surface of land, sea or in-land
waters is indicated by t2m. The skin temperature skt is the temperature of the
surface of the Earth, that is, the theoretical temperature that is required to
satisfy the surface energy balance. It represents the temperature of the
uppermost surface layer, which has no heat capacity and so can respond
instantaneously to changes in surface fluxes. All the soil temperatures stl1,
stl2, stl3, and stl4 are the temperature of the soil in the middle of that
given layer. The Advancing global NWP through international collaboration
(ECMWF) Integrated Forecasting System (IFS) has a four-layer representation of
soil, where the surface is at 0cm. Layer 1: 0 - 7cm; Layer 2: 7 - 28cm; Layer
3: 28 - 100cm; Layer 4: 100 - 289cm. Soil temperature is set at the middle of
each layer, and heat transfer is calculated at the interfaces between them. It
is assumed that there is no heat transfer out of the bottom of the lowest
layer. All the variables representing temperatures have units of kelvin (K),
but have been converted to degrees Celsius (C) by subtracting 273.15 in this
paper. They are calculated by interpolating between the lowest model level and
the Earth’s surface, taking account of the atmospheric conditions.

##  3 Spatio-temporal point processes

We consider a spatio-temporal point process with no multiple points as a
random countable subset $X$ of $\mathbb{R}^{2}\times\mathbb{R}$, where a point
$(u,t)\in X$ corresponds to an event at $u\in\mathbb{R}^{2}$ occurring at time
$t\in\mathbb{R}$. A typical realisation of a spatio-temporal point process $X$
on $\mathbb{R}^{2}\times\mathbb{R}$ is a finite set
$\\{(u_{i},t_{i})\\}^{n}_{i=1}$ of distinct points within a bounded spatio-
temporal region $W\times T\subset\mathbb{R}^{2}\times\mathbb{R}$, with area
$|W|>0$ and length $|T|>0$, where $n\geq 0$ is not fixed in advance. In this
context, $N(A\times B)$ denotes the number of points of a set $(A\times B)\cap
X$, where $A\subseteq W$ and $B\subseteq T$. As usual [13], when $N(W\times
T)<\infty$ with probability 1, which holds e.g. if $X$ is defined on a bounded
set, we call $X$ a finite spatio-temporal point process.

For a given event $(u,t)$, the events that are close to $(u,t)$ in both space
and time, for each spatial distance $r$ and time lag $h$, are given by the
corresponding spatio-temporal cylindrical neighbourhood of the event $(u,t)$,
which can be expressed by the Cartesian product as

$b((u,t),r,h)=\\{(v,s):||u-v||\leq r,|t-s|\leq h\\},\quad\quad(u,t),(v,s)\in
W\times T,$

where $||\cdot||$ denotes the Euclidean distance in $\mathbb{R}^{2}$ and
$|\cdot|$ is the absolute value. Note that $b((u,t),r,h)$ is a cylinder with
centre $(u,t)$, radius $r$, and height $2h$.

Product densities $\lambda^{(k)},\,k\in\mathbb{N}\text{ and }k\geq 1$,
arguably the main tools in the statistical analysis of point processes, may be
defined through the so-called Campbell Theorem (see [13]), that constitutes an
essential result in spatio-temporal point process theory. It states that,
given a spatio-temporal point process $X$, for any non-negative function $f$
on $(\mathbb{R}^{2}\times\mathbb{R})^{k}$

$\mathbb{E}\Bigg{[}\sum_{\zeta_{1},\dots,\zeta_{k}\in
X}^{\neq}f(\zeta_{1},\dots,\zeta_{k})\Bigg{]}=\int_{\mathbb{R}^{2}\times\mathbb{R}}\dots\int_{\mathbb{R}^{2}\times\mathbb{R}}f(\zeta_{1},\dots,\zeta_{k})\lambda^{(k)}(\zeta_{1},\dots,\zeta_{k})\prod_{i=1}^{k}\text{d}\zeta_{i},$

where $\neq$ indicates that the sum is over distinct values. In particular,
for $k=1$ these functions are called the intensity function $\lambda$. Broadly
speaking, the intensity function describes the rate at which the events occur
in the given spatio-temporal region, representing the point process analogues
of the mean function of a real-valued process. Then, the first-order intensity
function is formally defined as

$\lambda(u,t)=\lim_{|\text{d}u\times\text{d}t|\rightarrow
0}\frac{\mathbb{E}[N(\text{d}u\times\text{d}t)]}{|\text{d}u\times\text{d}t|},$

where $\text{d}u\times\text{d}t$ defines a small region around the point
$(u,t)$ and $|\text{d}u\times\text{d}t|$ is its volume.

###  3.1 Spatio-temporal Poisson point processes with separable intensity

The description of the observed point pattern intensity is a crucial issue
when dealing with spatio-temporal point pattern data, and specifying a
statistical model is a very effective way compared to analyzing data by
calculating summary statistics. Formulating and adapting a statistical model
to the data allows taking into account effects that otherwise could introduce
distortion in the analysis [5].

When dealing with intensity estimation for spatio-temporal point processes, it
is quite common to assume that the intensity function $\lambda(u,t)$ is
separable [19, 22]. Under this assumption, the intensity function is given by
the product

$\lambda(u,t)={\lambda}(u){\lambda}(t)$ (3)

where ${\lambda}(u)$ and ${\lambda}(t)$ are non-negative functions on $W$ and
$T$, respectively. Suitable estimates of $\lambda(u)$ and $\lambda(t)$ in
Equation (3) depend on the characteristics of each application. This
formulation can include a combination of a parametric spatial point pattern
model, potentially depending on the spatial coordinates and/or spatial
covariates, and a parametric log-linear model for the temporal component.
Also, non-parametric kernel estimate forms are legit. The spatio-temporal
intensity is therefore obtained by multiplying the purely spatial and purely
temporal intensities, previously fitted separately. The resulting intensity is
normalised, to make the estimator unbiased, making the expected number of
points

$\mathbb{E}\bigg{[}\int_{W\times
T}\hat{\lambda}(u,t)d_{2}(u,t)\bigg{]}=\int_{W\times
T}\lambda(u,t)d_{2}(u,t)=n,$

and the final intensity function is obtained as

$\hat{\lambda}(u,t)=\frac{\hat{\lambda}(u)\hat{\lambda}(t)}{\int_{W\times
T}\hat{\lambda}(u,t)d_{2}(u,t)}.$

###  3.2 Spatial component

A general spatial log-linear Poisson model [12], generalizes both homogeneous
and inhomogeneous models, such that:

$\lambda(u)=\exp(\boldsymbol{\theta}^{\top}\textbf{Z}(u)),$ (4)

where $u\in W$, and $\textbf{Z}(u)=\\{Z_{1}(u),\ldots,Z_{p}(u)\\}$ are known
covariate functions,
$\boldsymbol{\theta}^{\top}=(\theta_{1},\ldots,\theta_{p})$ are unknown
parameters, and therefore
$\boldsymbol{\theta}^{\top}\textbf{Z}(u)=(\theta_{1}Z_{1}(u)+\ldots+\theta_{p}Z_{p}(u))$.
These models have an especially convenient structure, since the log intensity
is a linear function of the parameters and covariates can be quite general
functions, making them a very wide class of models. The estimation of the
point process parameters is carried out through the maximization of the log-
likelihood, defined by:

$\text{log}\textrm{L}(\boldsymbol{\theta})=\sum_{i}\text{log}\lambda(x_{i};\boldsymbol{\theta})-\int_{W}\lambda(u;\boldsymbol{\theta})du$
(5)

where the sum is over all points $x_{i}$ in the point process x [13]. For
estimation purposes, we use a finite quadrature approximation of the log-
likelihood, following [7], that is, starting from the observed points
$u_{1},\dots,u_{n}$, generate $m$ additional dummy points
${u}_{n+1}\dots,{u}_{m+n}$ to form a set of $n+m$ quadrature points (where
$m>n$). Then, integral in (5) is approximated by the Riemann sum

$\int_{W}\lambda({u};\boldsymbol{\theta})\text{d}u\approx\sum_{k=1}^{n+m}a_{k}\lambda({u}_{k};\boldsymbol{\theta})$

where $a_{1},\dots,a_{m}$ are quadrature weights and
$\sum_{k=1}^{n+m}a_{k}=l(W)$, with $l(\cdot)$ the Lebesgue measure. In order
to distinguish between dummy points and observed points, we introduce the
indicator $e_{k}$, which assumes values equal to 1 if the related point is a
data point, zero otherwise. Then, writing $y_{k}=e_{k}/a_{k}$, the log-
likelihood (5) of the template model can be approximated by $\log
L(\theta)\approx\sum_{j}a_{k}(y_{k}\log\lambda({u}_{k};\boldsymbol{\theta})-\lambda({u}_{k};\boldsymbol{\theta}))+\sum_{k}a_{k}.$
Apart from the constant $\sum_{k}a_{k}$, this expression is formally
equivalent to the weighted log-likelihood of a Poisson regression model with
responses $y_{k}$ and means $\lambda({u}_{k},t_{k};\theta)=\exp(\theta
Z(u_{k},t_{k}))$. This means that the model can be maximised using standard
GLM software, but also that covariate values must be known in every data and
dummy point location. Indeed, the spatial covariates are referred to as those
variables with observable values, at least in principle, at each spatial
location in the spatial window, and for the inferential purposes given above,
their values must be known at each point of the data point pattern and at
least at some other locations.

####  3.2.1 Diagnostics

For an inhomogeneous Poisson process model, with fitted intensity
$\hat{\lambda}(u)$, the predicted number of points falling in any region $W$
is $\int_{W}\hat{\lambda}(u)$d$u$. Hence, the residual in each region
$W\subset\mathbb{R}^{2}$ is the ‘observed minus predicted’ number of points
falling in $W$ [4], that is $R(W)=n(\textbf{x}\cap
W)-\int_{W}\hat{\lambda}(u)du$, where x is the observed point pattern,
$n(\textbf{x}\cap W)$ the number of points of x in the region $W$, and
$\hat{\lambda}(u)$ is the intensity of the fitted model. A simple residual
visualization can be obtained by smoothing them. The ‘smoothed residual
fields’ are defined as

$s(u)=\tilde{\lambda}(u)-\lambda^{\dagger}(u)$ (6)

where
$\tilde{\lambda}(u)=e(u)\sum_{i=1}^{n(\textbf{x})}\kappa(u-\textbf{x}_{i})$ is
the non-parametric, kernel estimate of the fitted intensity
$\hat{\lambda}(u)$, while $\lambda^{\dagger}(u)$ is a correspondingly-smoothed
version of the (typically parametric) estimate of the intensity of the fitted
model, $\lambda^{\dagger}(u)=e(u)\int_{W}\kappa(u-v)\hat{\lambda}(v)dv$. Here,
$\kappa$ is the smoothing kernel and $e(u)$ is the edge correction. The
smoothing bandwidth for the kernel estimation of the raw residuals is selected
by cross-validation, as the value that minimises the Mean Squared Error
criterion defined by [18] by the method of [6].

The differences in Equation (6) should be approximately zero when the fitted
model is close to the real one. Therefore, the best model is the one with the
lowest values of the smoothed raw residuals. Their graphical representation
gives further insights on which spatial regions the model could not be a good
fit.

###  3.3 Temporal component

Temporal count data analysis often involves modelling the occurrence of events
over time, and a suitable statistical approach is essential to capture the
temporal dynamics adequately. In this context, the Poisson regression model,
which belongs to the class of Generalized Linear Models (GLMs), provides a
robust framework for analyzing count data exhibiting temporal variation.

Consider a set of covariates
$\boldsymbol{Z}(t)=\\{Z_{1}(t),\ldots,Z_{q}(t)\\}$ influencing the event rate.
The temporal count Poisson model is then formulated as follows:

$\lambda(t)=\exp(\boldsymbol{\beta}^{T}\boldsymbol{Z}(t)).$

Here, $t\in T$ and
$\boldsymbol{\beta}=\\{\beta_{0},\beta_{1},\ldots\beta_{q}\\}$ are the model
parameters associated with the temporal covariates.

The model parameters are estimated through the maximum likelihood estimation
(MLE) method within the GLM framework [27].

After fitting the model, assessing its performance involves standard GLM
techniques such as residual analysis, goodness-of-fit tests, and validation on
independent datasets. Interpretation of the model parameters provides insights
into the impact of covariates and temporal components on the event rate over
time.

A generalized additive model (GAM) is a GLM in which the linear predictor is
given by a sum of smooth functions of the covariates plus a conventional
parametric component. In such a case, we may write $\lambda(t)=\exp(f(Z(t)))$
where $f(\cdot)$ is a smooth function of the covariate $Z$. The smooth terms
can be functions of any number of covariates, and the control over the
smoothness of the functions is solved by using the Generalized Cross
Validation (GCV) criterion [40].

##  4 Results

In this section, we present the results obtained fitting the separable spatio-
temporal intensity function presented in Section 3 to the data introduced in
Section 2. Specifically, in Section 4.1, a spatial log-linear Poisson point
process model is applied to investigate the influence of land use types on the
purely spatial fire distribution, accounting for other environmental
covariates. In Section 4.2, a Poisson Generalized Additive Model is fitted to
the temporal fire occurrences, depending on both non-parametric components of
the temporal coordinates and parametric formulation of environmental
covariates.

###  4.1 Spatial intensity

Starting from the maximal model and proceeding with a backwards procedure
based on the comparison of the AIC values and the sequential performance of
tests for the significance of the model parameters, the chosen spatial model
has a linear predictor that includes a non-parametric term for the spatial
coordinates and parametric expression for the spatial covariates, as follows

$\lambda(u)=\exp(f(u)+\boldsymbol{\theta}_{1}Z_{land\\_use}(u)+\theta_{2}Z_{elev}(u)+\theta_{3}Z_{slope}(u)).$
(7)

Here, $f(\cdot)$ is a nonparametric function for $u\in W$, estimated through
thin plate regression splines [39] with 30 knots.

Table 3 presents the results from the fitted spatial Poisson model that
investigates the influence of various factors on fire occurrence. In addition
to the intercept, the first three rows correspond to the categories of the
land use variable, with Artificial surfaces serving as the baseline. Also note
that the elevation has been converted from meters to kilometres, to obtain
estimated coefficients on the same magnitude as the other ones.

| Estimate | Std. Error | z value | p value
---|---|---|---|---
Intercept | -15.129 | 0.040 | -376.580 | $<1\times 10^{-3}$
Agricultural areas | -0.286 | 0.040 | -7.179 | $<1\times 10^{-3}$
Forest and semi-natural areas | -0.540 | 0.046 | -11.727 | $<1\times 10^{-3}$
Wetlands and water bodies | -0.595 | 0.154 | -3.870 | $<1\times 10^{-3}$
Elevation | 1.675 | 0.040 | 42.320 | $<1\times 10^{-3}$
Slope | 0.006 | 0.001 | 4.768 | $<1\times 10^{-3}$
Table 3: Estimated coefficients of the fitted spatial model in Equation (7).

Considering different land usages, agricultural areas show a significant
negative estimated parameter, signifying a considerable decrease in the fire
occurrence intensity moving from artificial surfaces to agricultural ones.
Forest and semi-natural areas also show a negative effect, although with a
higher magnitude compared to agricultural areas. Then, wetlands and water
bodies also display a significant negative parameter, albeit with a wider
confidence interval, indicating greater uncertainty. This is likely due to the
sparsity of the area covered in Sicily by this type of land. Lastly, both
elevation and slope exhibit a significant effect, with their increase being
associated with higher intensity. Note that an attempt has been made to assess
the significance of the ERA5 spatio-temporal covariates, averaged over time.
However, these purely spatial covariates did not influence the overall spatial
occurrence of fires and therefore were not included in the chosen model. This
is likely due to the very coarse grid of the obtained spatial covariates if
compared to those representing elevation and slope (top panels of Figure 6).

Among the perks of fitting a semi-parametric model as the one in Equation (7)
there is the smoothed prediction map that can be obtained, as shown in the
left panel of Figure 7. The smoothed raw residuals in the right panel of
Figure 7, computed as detailed in Section 3.2.1, indicate a good model
fitting, exhibiting almost all residuals approaching zero.

Figure 7: Left panel: Intensity predicted according to the fitted spatial
model in Equation (7); Right panel: Smoothed raw residuals of the fitted
model.

###  4.2 Temporal intensity

Following the same backward procedure adopted for the spatial component, the
chosen temporal model comes with a linear predictor that includes a non-
parametric term for temporal coordinates and parametric expression for the
temporal covariates, as follows:

$\lambda(t)=\exp(f(t)+\beta_{1}Z_{v10}(t)+\beta_{2}Z_{stl2}(t)+\beta_{3}Z_{sp}(t)+\beta_{4}Z_{tp}(t)),$
(8)

with $f(\cdot)$ a nonparametric function for $t\in T$, estimated through
penalized regression basis splines [40] with 50 knots.

| Estimate | Std. Error | z value | p value
---|---|---|---|---
Intercept | 77.314 | 5.471 | 14.133 | $<1\times 10^{-3}$
Wind speed from South | 0.232 | 0.008 | 30.426 | $<1\times 10^{-3}$
Temperature | 0.440 | 0.024 | 18.672 | $<1\times 10^{-3}$
Surface pressure | -0.001 | 0.000 | -16.251 | $<1\times 10^{-3}$
Total precipitation | -14.848 | 2.266 | -6.553 | $<1\times 10^{-3}$
Table 4: Estimated coefficients of the fitted temporal model in Equation (8).

Figure 8: Daily maxima of the wind components, temperatures, precipitation and
surface pressure

In order to fit this model, the ERA5 environmental spatio-temporal covariates
have been averaged with respect to the spatial components, making them purely
temporal covariates. Moreover, the results of the models fitted with the daily
means of such covariates have been compared to those containing the daily
maxima, and the latter yielded better fitting. The resulting temporal
covariates are shown in Figure 8. As evident, the temperatures appear highly
correlated among themselves and with the time of the year. In detail, the
deeper the surface layers of detection, the smoother the temperature
variability. Figure 9 further illustrates such high correlations, especially
among temperature variables. Table 4 outlines the outcomes of the fitted
model, whose detailed breakdown follows.

Figure 9: Correlations among the temporal covariates and the daily fire counts
over 2023. On the main diagonal, the univariate distributions of the
variables. On the lower panels, the scatterplots of each pair of variables
with a smoothing function overlapped in red. On the upper panels, the
correlations among the variables.

Wind speed demonstrates a notable positive influence, implying that higher
wind speeds from the South correlate with increased fire counts, highlighting
the role of the Scirocco wind in wildfire dynamics. This Mediterranean wind
that comes from the Sahara brings heat and dust from African coastal regions.
Temperature shows a significant positive relationship with fire counts.
Pressure exhibits a negative estimated parameter, suggesting that higher
atmospheric pressure is associated with lower fire counts, hinting at its
potential suppressive effect. Precipitation reveals a substantial negative
effect, indicating fewer fire counts with higher daily precipitation,
highlighting its mitigating effect.

The significant smooth terms suggest that time plays a crucial role in the
observed residual variability of fire incidents, not taken into account by the
environmental covariates.

An adjusted R-squared of 0.738 and an 83.9% of deviance explained indicates
that the model explains substantial variability in fire counts. Figure 10
shows the intensity predicted according to the fitted model in Equation (8),
demonstrating a good ability to capture the temporal variability of fire
counts.

Figure 10: Temporal intensity predicted according to the fitted temporal model
in Equation (8) (red dashed line) and the observed fire counts (solid line).

##  5 Discussion and conclusions

The proposed analysis of the spatial and temporal distribution of fire
occurrences in Sicily during 2023 has provided insights into the factors
influencing these events, with a particular emphasis on the role of land
usage. Our findings highlight the significance of human activities in
contributing to the prevalence of fires. The identification of artificial
surfaces as a key contributor to the increased probability of fire occurrence
highlights the urgent need for targeted intervention and policy measures. The
implications of this research are particularly pertinent, considering the
well-established fact that the majority of fires in Sicily are human-induced
[20]. These results have broader implications for regional planning, resource
allocation, and the development of proactive measures to manage fire risks
effectively. Recognizing the human-centric nature of fire occurrences prompts
us to advocate for sustainable land-use practices and the implementation of
policies aimed at minimizing the risk of fires in areas characterized by
artificial surfaces. Our analysis extends beyond land usage, incorporating
environmental covariates that significantly influence the occurrence of fires
in Sicily during 2023. Our findings reveal a significant role of these
environmental factors over the temporal event occurrences. Particularly, wind
speed from the South and temperature emerge as crucial variables, emphasising
the importance of climatic conditions in amplifying fire risks. Conversely,
our study identifies surface pressure and total precipitation as environmental
covariates that negatively influence fire occurrences as mitigating factors.
Additionally, our research highlights the significance of terrain
characteristics in explaining the spatial distribution of fire points. In
particular, both elevation and slope exhibit a positive effect on the
occurrence of fires, indicating that areas with higher elevations and steeper
slopes are more susceptible.

Note that other FIRMS variables could have been used in this paper. These are
the brightness temperature in Kelvin, the along scan and track pixel sizes,
the type of satellite, a confidence value attached to each individual
hotspot/fire pixel, the collection/version source, the Fire Radiative Power,
in megawatts, and whether the fire occurred during daytime or nighttime. The
main reason for not using this type of covariates is that, in point process
theory, they are referred to as “marks” of the observed point pattern [13].
While the spatial covariates are referred to as those variables with
observable values, at least in principle, at each spatial location in the
spatio-temporal window, marks are characteristics of the events, which are
attached to the points of the observed patterns and can be studied to explore
their eventual aid in describing the phenomenon under study. Formally,
however, their inclusion in a regression-type model is quite different if
compared to the will to include the so-called spatial or spatio-temporal
covariates. Indeed, the latter covariates are characteristics of the spatial
or spatio-temporal region, and their inclusion in a regression-type point
process model is more straightforward, not needing any particular assumption
on, for instance, their distribution. Moreover, the literature on the so-
called marked point process methodologies is rather specific to particular
physical phenomena, like the seismic ones with the ETAS models [29], and quite
limited for the cases in the presence of multiple marks. This is why we
decided not to include marks in our analysis, even though we believe it
represents a very interesting topic, both from the methodological and the
applied points of view.

Moving to other possible future research paths, a major topic that could be
addressed is the definition and application of tailored diagnostic procedures
to assess the need for more complex models in case of residual clustering
behaviour of points not taken into account by the already considered factors.
Indeed, the employment of most known diagnostic tools based on second-order
summary statistics [1, 2], both global and local, becomes computationally hard
as the number of points increases. In a preliminary way, it is worth noticing
that one could first run a test for first-order separability prior to fitting
the model. In particular, [36], [16], and [21] show that the intensity of
forest fire occurrences varies in space and time in a nonseparable way. We,
however, believe that our proposed model overcomes computational burden
issues, typical of more complex non-separable spatio-temporal models.

Moreover, the inclusion and test of spatially (or spatio-temporally) varying
covariates in intensity function has been of particular interest in [17], [8,
10, 9] and [28]. Following these, future work includes the assessment of
separability assumption, as well as exploring more complex models, like the
log-Gaussian Cox processes, multitype Poisson models, and local ones [15].

## Fundings

The research work of Nicoletta D’Angelo and Giada Adelfio was supported by:

  * •

Targeted Research Funds 2024 (FFR 2024) of the University of Palermo (Italy);

  * •

Mobilità e Formazione Internazionali - Miur INT project “Sviluppo di
metodologie per processi di punto spazio-temporali marcati funzionali per la
previsione probabilistica dei terremoti”;

  * •

European Union - NextGenerationEU, in the framework of the GRINS -Growing
Resilient, INclusive and Sustainable project (GRINS PE00000018 – CUP
C93C22005270001).

Andrea Gilardi acknowledges the support by MUR, grant Dipartimento di
Eccellenza 2023-2027. His research work is funded by the European Union -
NextGenerationEU, in the framework of the GRINS - Growing Resilient, INclusive
and Sustainable project (GRINS PE00000018 – CUP D43C22003110001).

The research work of Alessandro Albano has been supported by the European
Union - NextGenerationEU - National Sustainable Mobility Center CN00000023,
Italian Ministry of University and Research Decree n. 1033— 17/06/2022, Spoke
2, CUP B73C2200076000 and Targeted Research Funds 2024 (FFR 2024) of the
University of Palermo (Italy).

The views and opinions expressed are solely those of the authors and do not
necessarily reflect those of the European Union, nor can the European Union be
held responsible for them.

## References

  * Adelfio and Schoenberg, [2009] Adelfio, G. and Schoenberg, F. P. (2009).  Point process diagnostics based on weighted second-order statistics and their asymptotic properties.  Annals of the Institute of Statistical Mathematics, 61(4):929–948.
  * Adelfio et al., [2020] Adelfio, G., Siino, M., Mateu, J., and Rodríguez-Cortés, F. J. (2020).  Some properties of local weighted second-order statistics for spatio-temporal point processes.  Stochastic Environmental Research and Risk Assessment, 34(1):149–168.
  * Aldersley et al., [2011] Aldersley, A., Murray, S. J., and Cornell, S. E. (2011).  Global and regional analysis of climate and human drivers of wildfire.  Science of the Total Environment, 409(18):3472–3481.
  * Alm, [1998] Alm, S. E. (1998).  Approximation and simulation of the distributions of scan statistics for poisson processes in higher dimensions.  Extremes, 1(1):111–126.
  * Baddeley et al., [2015] Baddeley, A., Rubak, E., and Turner, R. (2015).  Spatial point patterns: methodology and applications with R.  Chapman and Hall/CRC.
  * Berman and Diggle, [1989] Berman, M. and Diggle, P. (1989).  Estimating weighted integrals of the second-order intensity of a spatial point process.  Journal of the Royal Statistical Society: Series B (Methodological), 51(1):81–92.
  * Berman and Turner, [1992] Berman, M. and Turner, T. R. (1992).  Approximating point process likelihoods with glim.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 41(1):31–38.
  * Borrajo et al., [2017] Borrajo, M., González-Manteiga, W., and Martínez-Miranda, M. (2017).  Testing first-order intensity model in non-homogeneous poisson point processes with covariates.  arXiv preprint arXiv:1709.07716.
  * [9] Borrajo, M., González-Manteiga, W., and Martínez-Miranda, M. (2020a).  Testing for significant differences between two spatial patterns using covariates.  Spatial Statistics, 40:100379.
  * [10] Borrajo, M. I., González-Manteiga, W., and Martínez-Miranda, M. (2020b).  Bootstrapping kernel intensity estimation for inhomogeneous point processes with spatial covariates.  Computational Statistics & Data Analysis, 144:106875.
  * Butsic et al., [2015] Butsic, V., Kelly, M., and Moritz, M. A. (2015).  Land use and wildfire: A review of local interactions and teleconnections.  Land, 4(1):140–156.
  * Cox, [1972] Cox, D. (1972).  The statistical analysis of dependencies in point processes.  Stochastic Point Processes. Wiley: New York, pages 55–66.
  * Daley and Vere-Jones, [2007] Daley, D. J. and Vere-Jones, D. (2007).  An Introduction to the Theory of Point Processes. Volume II: General Theory and Structure.  Springer-Verlag, New York, second edition.
  * D’Angelo and Adelfio, [2023] D’Angelo, N. and Adelfio, G. (2023).  stopp: Spatio-Temporal Point Pattern Methods, Model Fitting, Diagnostics, Simulation, Local Tests.  R package version 0.1.0.
  * D’Angelo et al., [2023] D’Angelo, N., Adelfio, G., and Mateu, J. (2023).  Locally weighted minimum contrast estimation for spatio-temporal log-gaussian cox processes.  Computational Statistics & Data Analysis, 180:107679.
  * Díaz-Avalos et al., [2013] Díaz-Avalos, C., Juan, P., and Mateu, J. (2013).  Similarity measures of conditional intensity functions to test separability in multidimensional point processes.  Stochastic Environmental Research and Risk Assessment, 27:1193–1205.
  * Díaz-Avalos et al., [2014] Díaz-Avalos, C., Juan, P., and Mateu, J. (2014).  Significance tests for covariate-dependent trends in inhomogeneous spatio-temporal point processes.  Stochastic environmental research and risk assessment, 28:593–609.
  * Diggle, [1985] Diggle, P. (1985).  A kernel method for smoothing point process data.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 34(2):138–147.
  * Diggle, [2013] Diggle, P. J. (2013).  Statistical analysis of spatial and spatio-temporal point patterns.  Chapman and Hall/CRC.
  * Ferrara et al., [2019] Ferrara, C., Salvati, L., Corona, P., Romano, R., and Marchi, M. (2019).  The background context matters: Local-scale socioeconomic conditions and the spatial distribution of wildfires in italy.  Science of the Total Environment, 654:43–52.
  * Fuentes-Santos et al., [2018] Fuentes-Santos, I., González-Manteiga, W., and Mateu, J. (2018).  A first-order, ratio-based nonparametric separability test for spatiotemporal point processes.  Environmetrics, 29(1):e2482.
  * Gabriel and Diggle, [2009] Gabriel, E. and Diggle, P. J. (2009).  Second-order analysis of inhomogeneous spatio-temporal point process data.  Statistica Neerlandica, 63(1):43–51.
  * Ganteaume et al., [2013] Ganteaume, A., Camia, A., Jappiot, M., San-Miguel-Ayanz, J., Long-Fournel, M., and Lampin, C. (2013).  A review of the main driving factors of forest fire ignition over europe.  Environmental management, 51:651–662.
  * GDAL/OGR contributors, [2023] GDAL/OGR contributors (2023).  GDAL/OGR Geospatial Data Abstraction software Library.  Open Source Geospatial Foundation.
  * Hantson et al., [2015] Hantson, S., Pueyo, S., and Chuvieco, E. (2015).  Global fire size distribution is driven by human impact and climate.  Global Ecology and Biogeography, 24(1):77–86.
  * Horn, [1981] Horn, B. K. (1981).  Hill shading and the reflectance map.  Proceedings of the IEEE, 69(1):14–47.
  * McCullagh and Nelder, [1989] McCullagh, P. and Nelder, J. A. (1989).  Generalised Linear Models.  Chapman and Hall, 2nd edition.
  * Myllymäki et al., [2021] Myllymäki, M., Kuronen, M., and Mrkvička, T. (2021).  Testing global and local dependence of point patterns on covariates in parametric models.  Spatial Statistics, 42:100436.
  * Ogata, [1988] Ogata, Y. (1988).  Statistical models for earthquake occurrences and residual analysis for point processes.  Journal of the American Statistical Association, 83(401):9–27.
  * Pebesma, [2018] Pebesma, E. (2018).  Simple Features for R: Standardized Support for Spatial Vector Data.  The R Journal, 10(1):439–446.
  * Pebesma and Bivand, [2023] Pebesma, E. and Bivand, R. (2023).  Spatial Data Science: With applications in R.  Chapman and Hall/CRC, London.
  * R Core Team, [2024] R Core Team (2024).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Ricotta et al., [2018] Ricotta, C., Bajocco, S., Guglietta, D., and Conedera, M. (2018).  Assessing the influence of roads on fire ignition: does land cover matter?  Fire, 1(2):24.
  * Ricotta and Di Vito, [2014] Ricotta, C. and Di Vito, S. (2014).  Modeling the landscape drivers of fire recurrence in sardinia (italy).  Environmental management, 53(6):1077–1084.
  * Rodrigues and De la Riva, [2014] Rodrigues, M. and De la Riva, J. (2014).  An insight into machine-learning algorithms to model human-caused wildfire occurrence.  Environmental Modelling & Software, 57:192–201.
  * Schoenberg, [2004] Schoenberg, F. P. (2004).  Testing separability in spatial-temporal marked point processes.  Biometrics, pages 471–481.
  * Tarquini et al., [2023] Tarquini, S., Isola, I., Favalli, M., Battistini, A., and Dotta, G. (2023).  TINITALY, a digital elevation model of Italy with a 10 meters cell size, 1.1 edition.  Istituto Nazionale di Geofisica e Vulcanologia (INGV). https://doi.org/10.13127/tinitaly/1.1.
  * Weih Jr and Mattson, [2004] Weih Jr, R. C. and Mattson, T. L. (2004).  Modeling slope in a geographic information system.  Journal of the Arkansas Academy of Science, 58(1):100–108.
  * Wood, [2003] Wood, S. N. (2003).  Thin plate regression splines.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 65(1):95–114.
  * Wood, [2017] Wood, S. N. (2017).  Generalized additive models: an introduction with R.  CRC press.


# Exponential tractability of $L_{2}$-approximation
with function values

David Krieg1,2 and Paweł Siedlecki3 and Mario Ullrich1 and Henryk
Woźniakowski3,4 1 Institut für Analysis, Johannes Kepler Universität, Linz,
Austria  2 Johann Radon Institute for Computational and Applied Mathematics
(RICAM), Austrian Academy of Sciences, Linz, Austria.  3 Faculty of
Mathematics, Informatics and Mechanics, University of Warsaw, Warsaw, Poland
4 Department of Computer Science, Columbia University, New York, United States
of America<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We study the complexity of high-dimensional approximation in the $L_{2}$-norm
when different classes of information are available; we compare the power of
function evaluations with the power of arbtirary continuous linear
measurements. Here, we discuss the situation when the number of linear
measurements required to achieve an error $\varepsilon\in(0,1)$ in dimension
$d\in\mathbb{N}$ depends only poly-logarithmically on $\varepsilon^{-1}$. This
corresponds to an exponential order of convergence of the approximation error,
which often happens in applications. However, it does not mean that the high-
dimensional approximation problem is easy, the main difficulty usually lies
within the dependence on the dimension $d$. We determine to which extent the
required amount of information changes if we allow only function evaluation
instead of arbitrary linear information. It turns out that in this case we
only lose very little, and we can even restrict to linear algorithms. In
particular, several notions of tractability hold simultaneously for both types
of available information.

Keywords. Approximation, Multivariate problems, Tractability, Complexity

##  1\. Exposition

We want to approximate real- or complex-valued functions defined on some
(nonempty) set $\mathcal{D}$, and belonging to a space $F$. We assume that $F$
is a separable Banach space of functions defined on $\mathcal{D}$, such that
function evaluation $f\mapsto f(x)$ is continuous on $F$ for each
$x\in\mathcal{D}$ and $F$ is continuously embedded in
$L_{2}=L_{2}(\mathcal{D},\mu)$ for some measure $\mu$. Formally, the
approximation problem is given as

$\operatorname{APP}:F\to L_{2},\quad\operatorname{APP}(f):=f,$

which might be understood as a continuous embedding into $L_{2}$. The class of
all spaces $F$ satisfying the assumptions above will be denoted by
$\mathcal{A}$. In particular, for each $F\in\mathcal{A}$ we have some
associated nonempty set $\mathcal{D}$, measure $\mu$ on $\mathcal{D}$ and
continuous embedding $\operatorname{APP}$.

We approximate $\operatorname{APP}$ by using functionals from the class
$\Lambda^{\mathrm{std}}$ consisting of all function evaluations, or from the
class $\Lambda^{\mathrm{all}}=F^{\ast}$ of all continuous linear functionals.

Below $B_{F}$ denotes the closed unit ball in $F$. Let us define, for
$n\in\mathbb{N}$, the

  * •

_$n$ -th linear sampling width_ as

$\displaystyle
e_{n}(F,L_{2})\,:=\,\inf_{\begin{subarray}{c}x_{1},\dots,x_{n}\in\mathcal{D}\\\
\varphi_{1},\dots,\varphi_{n}\in L_{2}\end{subarray}}\,\sup_{f\in
B_{F}}\,\Big{\|}f-\sum_{i=1}^{n}f(x_{i})\,\varphi_{i}\Big{\|}_{L_{2}},$

  * •

_$n$ -th sampling width_ as

$\displaystyle
g_{n}(F,L_{2})\,:=\,\inf_{\begin{subarray}{c}x_{1},\dots,x_{n}\in\mathcal{D}\\\
\phi\colon\mathbb{R}^{n}\to L_{2}\\\ \end{subarray}}\,\sup_{f\in
B_{F}}\,\Big{\|}f-\phi(f(x_{1}),\dotsc,f(x_{n}))\Big{\|}_{L_{2}},$

  * •

_$n$ -th linear width_ as

$\displaystyle a_{n}(F,L_{2})\,:=\,\inf_{\begin{subarray}{c}T\colon L_{2}\to
L_{2}\\\ \operatorname{rank}(T)\,\leq\,n\end{subarray}}\,\sup_{f\in
B_{F}}\,\big{\|}f-Tf\big{\|}_{L_{2}},$

  * •

_$n$ -th Gelfand width_ as

$\displaystyle
c_{n}(F,L_{2})\,:=\,\inf_{\begin{subarray}{c}\phi\colon\mathbb{R}^{n}\to
L_{2}\\\ N\in(F^{*})^{n}\end{subarray}}\,\sup_{f\in
B_{F}}\,\big{\|}f-\phi\circ N(f)\big{\|}_{L_{2}}.$

These quantities represent the _minimal worst case errors_ that can be
achieved with linear or nonlinear algorithms using at most $n$ function values
or linear measurements, respectively.

We also define the information-based complexity of the problem
$\operatorname{APP}$ for the classes $\Lambda^{\mathrm{std}}$ and
$\Lambda^{\mathrm{all}}$, respectively, as the minimal number of evaluations
from $\Lambda^{\mathrm{std}}$ or $\Lambda^{\mathrm{all}}$ necessary to obtain
the absolute precision of approximation at most $\varepsilon$, i.e., as

$n^{\rm std}(\varepsilon,F)\,:=\,\min\big{\\{}n\colon
g_{n}(F,L_{2})\leq\varepsilon\big{\\}}$

and

$n^{\rm all}(\varepsilon,F)\,:=\,\min\big{\\{}n\colon
c_{n}(F,L_{2})\leq\varepsilon\big{\\}}.$

Note that, since $g_{n}(F,L_{2})\leq e_{n}(F,L_{2})$, we have

$n^{\rm std}(\varepsilon,F)\,\leq\,\min\big{\\{}n\colon
e_{n}(F,L_{2})\leq\varepsilon\big{\\}}\,\,=:\,n^{\text{\rm std-
lin}}(\varepsilon,F),$

and all our upper bounds are proven for $n^{\text{\rm std-
lin}}(\varepsilon,F)$. There is a lot of literature on the size of these
quantities for specific classes $F$. We refer to the monographs [10, 33, 35,
36, 42, 43, 45] for more details and literature on the subject.

Here, we are specifically interested in the comparison of these quantities for
general classes $F$. That is, since $n^{\rm all}(\varepsilon,F)\leq n^{\rm
std}(\varepsilon,F)$ is obvious for all $F\in\mathcal{A}$, we ask for an upper
bound on $n^{\rm std}(\varepsilon,F)$ based on knowledge of the function
$n^{\rm all}(\varepsilon,F)$. However, it is known that such a bound cannot
hold without certain assumptions on $F$, see [36, Chapter 26] and references
therein, and even then, the involved “constants” depend in a non-trivial way
on $F$. One approach to obtain qualitative statements on the relation of the
complexities is to consider a whole sequence of spaces
$(F_{d})_{d\in\mathbb{N}}$, where $d$ can be interpreted as the dimension of
the underlying domain. We then assume a certain bound on $n^{\rm
all}(\varepsilon,F_{d})$, depending only on $\varepsilon\in(0,1)$ and
$d\in\mathbb{N}$, and ask for an upper bound on $n^{\rm
std}(\varepsilon,F_{d})$, hopefully not much worse than the bound on $n^{\rm
all}(\varepsilon,F_{d})$.

In the present paper, we allow arbitrary Banach spaces of functions $F_{d}$,
but we assume that $n^{\rm all}(\varepsilon,F_{d})$ depends only poly-
logarithmically on $\varepsilon^{-1}$. That is, we assume that there exist
$A_{d},B_{d}>0$ such that

(1) $n^{\rm all}(\varepsilon,F_{d})\leq
A_{d}\left(1+\ln\varepsilon^{-1}\right)^{B_{d}}\qquad\text{for all
}0<\varepsilon\leq 1,$

and study how this translates into bounds on $n^{\rm std}(\varepsilon,F_{d})$.
Note that the above bound (1) on the complexity implies that

(2) $c_{n}(F_{d},L_{2})\leq e\,\exp(-(n/A_{d})^{1/B_{d}})\qquad\text{for all
}n\geq A_{d},$

whereas (2) implies that (1) holds with $+1$ added on the right hand side. The
assumption (1) is therefore equivalent to the existence of a (possibly non-
linear) algorithm based on arbitrary linear information that converges
exponentially fast. We will show that in this case, we do not lose much when
we only allow linear algorithms and function evaluations as information. One
of our main results may be stated as follows.

######  Theorem (see Corollary 5).

Assume that $F_{d}\in\mathcal{A}$ for every $d\in\mathbb{N}$ and

$\displaystyle n^{\rm
all}(\varepsilon,F_{d})\,\leq\,c\,d^{q}\,(1+\ln\varepsilon^{-1})^{p}$

for some $p,c>0$, $q\geq 0$, and all $\varepsilon\in(0,1)$. Then

$\displaystyle n^{\text{\rm std-
lin}}(\varepsilon,F_{d})\,\leq\,C\,d^{q}\,(1+\ln
d)^{p}\,(1+\ln\varepsilon^{-1})^{p}$

for all $\varepsilon\in(0,1)$ and $d\in\mathbb{N}$, and some $C>0$ that
depends only on $c$, $p$ and $q$.

This shows that every Banach space that is assumed to be _approximable in high
dimensions_ (in the above sense) with an exponential order by some algorithm
and information can practically already be treated with linear algorithms
based on function values. In particular, this improves upon Theorem 26.21 from
[36] and solves Open Problem 128 therein. Let us add that we do not know if
the additional $(1+\ln d)^{p}$ is necessary.

There are many appearances of the assumption (1) in the literature. Besides
the detailed study of certain weighted Hilbert spaces of analytic functions
[8, 24, 25, 30, 47], it appears naturally in the context of approximation with
(increasingly flat) Gaussian kernels [12, 17, 26, 41], or in tensor product
approximations [14, 16], or for certain smoothness spaces on complex spheres
[7]. Moreover, it is a typical assumption for the construction of greedy bases
[4, 5, 15]. Let us also add that there is quite some study on the _stability_
of algorithms that can achieve an exponential convergence, see [1, 2, 3, 39]
for details.

When it comes to the study of the _tractability_ of the problem, i.e., the
precise dependence of the error on the dimension, especially when we only
allow function evaluations, there is much less to cite and we are only aware
of the Hilbert space references from above. As an explicit example, let us
mention the Gaussian space on $\mathbb{R}^{d}$ with reproducing kernel
$K(x,y)=\exp(-\|x-y\|_{2}^{d})$, which satisfies a relation of the form (1)
for $L_{2}$-approximation with respect to the Gaussian probability measure
$\mu$, see [41]. In the Hilbert case, there are some general results which
make the situation somewhat simpler. For example, it is known that linear
algorithms are always optimal and one may work with the singular value
decomposition of the embedding $\operatorname{APP}$. We refer to [33, Chapter
4] and [36, Chapter 26].

A bit more is known in the case of _algebraic tractability_ , i.e., when the
complexity depends polynomially on $\varepsilon^{-1}$ instead of
$\ln\varepsilon^{-1}$. In addition to general Hilbert space results, see [36,
Chapter 26], and characterizations for weighted Korobov spaces, see [11],
there are also quite sharp results for the classical smoothness spaces
$C^{k}(\Omega_{d})$ of $k$-times differentiable functions on certain
$d$-dimensional domains, possibly for $k=\infty$. See [27, 34, 46] for details
on approximation, and [20, 21, 22, 23] for numerical integration in the same
classes. However, a comparison as proven here in the case of exponential
convergence is not possible in this case, see the end of Section 2. In any
case, it is open to determine the precise behavior of $n^{\rm
std}(\varepsilon,F_{d})$ for most classical spaces, while $n^{\rm
all}(\varepsilon,F_{d})$ is more often known.

Our results are based on the following (special case of) Theorem 3 from [9],
see also [18, 28, 29, 32, 44], which allows us to treat more general classes
of functions.

######  Theorem 1.

For each $0<r<2$, there is a universal constant $b\in\mathbb{N}$, depending
only on $r$, such that the following holds. For all $F\in\mathcal{A}$ and
$n\geq 2$, we have

$e_{bn}(F,L_{2})\,\leq\,\left(\frac{1}{n}\sum_{k\geq
n}a_{k}(F,L_{2})^{r}\right)^{1/r}.$

Additionally, we use the following fundamental result from [38], see also [6,
31].

######  Theorem 2.

For all $F\in\mathcal{A}$ and $n\geq 1$, we have

$a_{n}(F,L_{2})\,\leq\,\left(1+\sqrt{n}\,\right)c_{n}(F,L_{2}).$

######  Remark 3.

We would like to stress that the proof of Theorem 1 in [9] is non-
constructive, and we do not know how to explicitly construct evaluation points
$x_{1},\dotsc,x_{bn}\in\mathcal{D}$ together with some
$\varphi_{1},\dotsc,\varphi_{bn}\in L_{2}$ satisfying

$\sup_{f\in
B_{F}}\,\Big{\|}f-\sum_{i=1}^{bn}f(x_{i})\,\varphi_{i}\Big{\|}_{L_{2}}\leq\left(\frac{1}{n}\sum_{k\geq
n}a_{k}(F,L_{2})^{r}\right)^{1/r}.$

However, for problems with known operators $T$ achieving the infimum as in the
definition of linear widths $a_{n}(F,L_{2})$, we are able to specify
algorithms utilizing i.i.d. sampling of evaluation points from some known
distribution, and satisfying inequalities similar to the one above with high
probability, see Theorem 8 of [29].

##  2\. Exponential tractability of approximation

The notions of tractability are defined as follows. Let us fix, for every
$d\in\mathbb{N}$, some space $F_{d}\in\mathcal{A}$. For each
$F_{d}\in\mathcal{A}$ we have some associated set $\mathcal{D}_{d}$ equipped
with a measure $\mu_{d}$, and a continuous embedding
$\operatorname{APP}_{d}:F_{d}\to L_{2}(\mathcal{D}_{d},\mu_{d})$. The index
$d\in\mathbb{N}$ is an arbitrary parameter, but it usually stands for the
dimension of the domain $\mathcal{D}_{d}$. _A multivariate approximation
problem_ is simply a sequence of embeddings

$\widetilde{\mathrm{APP}}=\big{(}\mathrm{APP}_{d}:F_{d}\to
L_{2}(\mathcal{D}_{d},\mu_{d})\big{)}_{d\in\mathbb{N}}.$

Moreover, tractability notions are defined relative to the considered class of
information operations, i.e., we can consider tractability for
$\Lambda^{\mathrm{std}}$ or $\Lambda^{\mathrm{all}}$. Therefore, for
$\mathrm{x}\in\left\\{\mathrm{std},\mathrm{all}\right\\}$, we say that
$\widetilde{\mathrm{APP}}$ is

  * •

_exponentially strongly polynomially tractable_ (EXP-SPT) for the class
$\Lambda^{\mathrm{x}}$ if and only if

$n^{\rm x}(\varepsilon,F_{d})\,\leq\,C\,(1+\ln\varepsilon^{-1})^{p}$

for some $C,p>0$ and for all $d\in\mathbb{N}$ and $\varepsilon\in(0,1)$,

  * •

_exponentially polynomially tractable_ (EXP-PT) for the class
$\Lambda^{\mathrm{x}}$ if and only if

$n^{\rm x}(\varepsilon,F_{d})\,\leq\,C\,d^{q}\,(1+\ln\varepsilon^{-1})^{p}$

for some $C,p,q>0$ and for all $d\in\mathbb{N}$ and $\varepsilon\in(0,1)$,

  * •

_exponentially quasi-polynomially tractable_ (EXP-QPT) for the class
$\Lambda^{\mathrm{x}}$ if and only if

$n^{\rm x}(\varepsilon,F_{d})\,\leq\,C\,\exp(t(1+\ln
d)(1+\ln(1+\ln\varepsilon^{-1})))$

for some $C,t>0$ and for all $d\in\mathbb{N}$ and $\varepsilon\in(0,1)$,

  * •

_exponentially uniformly weakly tractable_ (EXP-UWT) for the class
$\Lambda^{\mathrm{x}}$ if and only if for all $\alpha,\beta>0$ we have

$\lim_{d+\varepsilon^{-1}\to\infty}\,\frac{\ln n^{\rm
x}(\varepsilon,F_{d})}{d^{\,\alpha}+(1+\ln\varepsilon^{-1})^{\beta}}=0,$

  * •

_exponentially weakly tractable (EXP-WT)_ for the class $\Lambda^{\mathrm{x}}$
if and only if

$\lim_{d+\varepsilon^{-1}\to\infty}\,\frac{\ln n^{\rm
x}(\varepsilon,F_{d})}{d+(1+\ln\varepsilon^{-1})}=0.$

It is easy to see that we have the following logical relation between the
tractability notions defined above

$\text{EXP-SPT}\implies\text{EXP-PT}\implies\text{EXP-QPT}\implies\text{EXP-
UWT}\implies\text{EXP-WT}.$

For a multivariate approximation problem we prove that _exponential strong
polynomial tractability_ (EXP-SPT), _exponential polynomial tractability_
(EXP-PT), _exponential uniform weak tractability_ (EXP-UWT) and _exponential
weak tractability_ (EXP-WT) for the class $\Lambda^{\mathrm{all}}$ are each
equivalent to the corresponding tractability property for the class
$\Lambda^{\mathrm{std}}$. Moreover, _exponential quasi-polynomial
tractability_ (EXP-QPT) for $\Lambda^{\mathrm{all}}$ implies _exponential
uniform weak tractability_ (EXP-UWT) for $\Lambda^{\mathrm{std}}$, i.e, the
next tractability notion in the tractability hierarchy considered here.
Whether the equivalence of _exponential quasi-polynomial tractability_ (EXP-
QPT) for the classes $\Lambda^{\mathrm{all}}$ and $\Lambda^{\mathrm{std}}$
holds remains an open problem.

These equivalences are in sharp contrast to the results for algebraic
tractability. See, e.g., [19, 37, 40] for examples where the problem is
algebraically tractable for $\Lambda^{\mathrm{all}}$ but the curse of
dimensionality holds for $\Lambda^{\mathrm{std}}$. In particular, [37, Example
5] shows that for the tensor product $W_{2,d}^{s}$ of certain univariate
periodic Sobolev spaces, $s>1/2$, we have QPT for $\Lambda^{\mathrm{all}}$,
but the curse of dimensionality for $\Lambda^{\mathrm{std}}$.

##  3\. Results

We now present our results. The first results are concerned with EXP-(S)PT and
EXP-QPT. Both are direct corollaries of the following theorem.

######  Theorem 4.

Assume that $F\in\mathcal{A}$ satisfies

$\displaystyle n^{\rm
all}(\varepsilon,F)\,\leq\,A\,\bigl{(}1+\ln\varepsilon^{-1}\bigr{)}^{B}$

for some $B>0$ and $A\geq 1$ and all $\varepsilon\in(0,1)$. Then

$\displaystyle n^{\rm std}(\varepsilon,F)\,\leq\,n^{\text{\rm std-
lin}}(\varepsilon,F_{d})\,\leq\,C\,\bigl{(}1+\ln\varepsilon^{-1}\bigr{)}^{B}$

for all $\varepsilon\in(0,1)$, where

$C=3b\,A\,\left(\ln(36A)\,(1+B^{3})\right)^{B}$

and $b$ is the absolute constant from Theorem 1 in the case $r=1$.

###### Proof.

Observe that the inequality

$n^{\rm all}(\varepsilon,F)\leq A\bigl{(}1+\ln\varepsilon^{-1}\bigr{)}^{B}$

implies that

$c_{n}(F,L_{2})\leq e\,\exp(-(n/A)^{1/B}).$

We obtain from Theorem 2, and $1+n^{1/2}\leq 2n^{1/2}$, that

$a_{n}(F,L_{2})\leq 2\,e\,n^{1/2}\,\exp(-(n/A)^{1/B}).$

Applying first Lemma 8 and then Lemma 9 from the Appendix, we deduce that

$\displaystyle\sum_{k\geq n}a_{k}(F,L_{2})$ $\displaystyle\leq 2\,e\sum_{k\geq
n}k^{1/2}\exp(-(k/A)^{1/B})$ $\displaystyle\leq
6\,A^{1/B}\,B\,\max(3B/2,1)\,(n-1)^{3/2-1/B}\,\exp(-((n-1)/A)^{1/B})$

for all $n\geq n_{0}(A,B):=A\,\max(3B/2,1)^{B}+1$. In particular,
$(a_{n}(F,L_{2}))\in\ell_{1}$. It follows from Theorem 1 that there exists an
absolute constant $b\in\mathbb{N}$ such that

$\displaystyle e_{bn}(F,L_{2})$ $\displaystyle\leq\,n^{-1}\,\sum_{k\geq
n}a_{k}(F,L_{2})$
$\displaystyle\leq\,6\,A^{1/B}\,B\,\max(3B/2,1)\,(n-1)^{1/2-1/B}\,\exp(-((n-1)/A)^{1/B})$

for all $n\geq n_{0}(A,B)$.

In the case $B\leq 2$ we have $(n-1)^{1/2-1/B}\leq A^{1/2-1/B}$, and thus

$e_{bn}(F,L_{2})\,\leq\,36\,A^{1/2}\,\exp(-((n-1)/A)^{1/B}).$

If $B>2$, then Lemma 10 with $u=1/2-1/B$ yields for any $\delta>0$ that

$\displaystyle e_{bn}(F,L_{2})\leq
9\,A^{1/2}\,B^{2}\,\delta^{1-B/2}\,\exp\left(((B/2-1)\delta-1)((n-1)/A)^{1/B}\right)$

and taking $\delta=2/B$ yields

$\displaystyle
e_{bn}(F_{d},L_{2})\,\leq\,36\,A^{1/2}(B/2)^{B/2+1}\exp\left(-\frac{2}{B}\left(\frac{n-1}{A}\right)^{1/B}\right).$

If we put $B_{0}:=\max\\{B/2,1\\}$, we have for all $B>0$ and $n\geq
n_{0}(A,B)$ the bound

$e_{bn}(F_{d},L_{2})\,\leq\,36\,A^{1/2}B_{0}^{B_{0}+1}\exp\left(-\frac{1}{B_{0}}\left(\frac{n-1}{A}\right)^{1/B}\right)$

which is smaller than $\varepsilon$ if

$n\,\geq\,A\,B_{0}^{B}\left(\ln\left(36\,A^{1/2}B_{0}^{B_{0}+1}\varepsilon^{-1}\right)\right)^{B}+1.$

Thus

$\displaystyle n^{\rm std}(\varepsilon,F_{d})$
$\displaystyle\leq\,b\,\max\left\\{A\,B_{0}^{B}\left(\ln\left(36\,A^{1/2}B_{0}^{B_{0}+1}\varepsilon^{-1}\right)\right)^{B}+2,\
n_{0}(A,B)\right\\}$
$\displaystyle\leq\,3b\,\max\left\\{A\,B_{0}^{B}R^{B}\left(1+\ln(\varepsilon^{-1})\right)^{B},\
A\,(3B_{0})^{B}\right\\}$
$\displaystyle\leq\,3b\,A\,B_{0}^{B}R^{B}\left(1+\ln(\varepsilon^{-1})\right)^{B}$

with

$R\,:=\,\ln(36)+\frac{\ln A}{2}+(B_{0}+1)\ln B_{0}\,\leq\,\ln(36A)\,B_{0}^{2}$

which gives the desired estimate.
∎

######  Corollary 5.

Assume that $F_{d}\in\mathcal{A}$ for every $d\in\mathbb{N}$ and

$\displaystyle n^{\rm
all}(\varepsilon,F_{d})\,\leq\,c\,d^{q}\,(1+\ln\varepsilon^{-1})^{p}$

for some $p,c>0$, $q\geq 0$, and all $\varepsilon\in(0,1)$. Then

$\displaystyle n^{\rm std}(\varepsilon,F_{d})\,\leq\,n^{\text{\rm std-
lin}}(\varepsilon,F_{d})\,\leq\,C\,d^{q}\,(1+\ln
d)^{p}\,(1+\ln\varepsilon^{-1})^{p}$

for all $\varepsilon\in(0,1)$ and $d\in\mathbb{N}$, and some $C>0$ that
depends only on $c$, $p$ and $q$.

In particular, if $\widetilde{\operatorname{APP}}$ is exponentially (strongly)
polynomially tractable for the class $\Lambda^{\mathrm{all}}$ then it is
exponentially (strongly) polynomially tractable for $\Lambda^{\mathrm{std}}$.

###### Proof.

We use Theorem 4 with $A=c\,d^{q}+1$ and $B=p$.
∎

We now turn to the assumption that $\widetilde{\mathrm{APP}}$ is exponentially
quasi-polynomially tractable for the class $\Lambda^{\mathrm{all}}$. This is
the only case where we do not know if it implies the same property for
$\Lambda^{\mathrm{std}}$.

For convenience, let us write $\ln_{+}(x):=1+\ln(x)$.

######  Corollary 6.

Assume that $F_{d}\in\mathcal{A}$ for every $d\in\mathbb{N}$ and

$\displaystyle n^{\rm
all}(\varepsilon,F_{d})\,\leq\,c\,\exp\left(t\cdot\ln_{+}d\cdot\ln_{+}\ln_{+}\varepsilon^{-1}\right)$

for some $c,t>0$ and all $\varepsilon\in(0,1)$. Then

$\displaystyle n^{\rm
std}(\varepsilon,F_{d})\,\leq\,c\,\exp\left(t\cdot\ln_{+}d\cdot\Bigl{(}\ln_{+}\ln_{+}\varepsilon^{-1}\,+\,4\ln\bigl{(}t\,\ln_{+}d\bigr{)}\,+\,C\Bigr{)}\right)$

for all $\varepsilon\in(0,1)$ and $d>(e+\frac{1}{c})^{1/t}e^{-1}$, and some
$C>0$ that depends only on $c$.

In particular, if $\widetilde{\mathrm{APP}}$ is exponentially quasi-
polynomially tractable for the class $\Lambda^{\mathrm{all}}$, then it is
exponentially uniformly weakly tractable for the class $\Lambda^{std}$.

###### Proof.

Note that

$\displaystyle c\,\exp\left(t\,(1+\ln
d)\,(1+\ln(1+\ln\varepsilon^{-1}))\right)\,=\,c\,e^{t}\,d^{\,t}\,\left(1+\ln\varepsilon^{-1}\right)^{t(1+\ln
d)}.$

Hence, we can apply Theorem 4 with $A=c\,e^{t}\,d^{\,t}$ and $B=t\,\ln_{+}d$,
i.e., $A=c\,e^{B}$. Note that $A,B\geq 1$ for $d>(e+\frac{1}{c})^{1/t}e^{-1}$.
We obtain that there exists an absolute constant $b>0$ such that

$\displaystyle n^{\rm
std}(\varepsilon,F)\,\leq\,C\,\bigl{(}1+\ln\varepsilon^{-1}\bigr{)}^{B}$

for all $\varepsilon\in(0,1)$, with

$\begin{split}C\,&=\,3b\,A\,\left(\ln(36A)\,(1+B^{3})\right)^{B}\leq\,3b\,A\,(2B)^{3B}\left(\ln(36c)+B\right)^{B}\\\
&\leq\,c\,\exp\left(B\Bigl{(}c^{\prime}+4\ln(B)\Bigr{)}\right)\end{split}$

where $c^{\prime}>0$ only depends on $c$. This proves the bound.

Now, since $\ln n^{\rm std}(\varepsilon,F_{d})$ depends only logarithmically
on $d$ and double-logarithmically on $\varepsilon^{-1}$, we obtain

$\lim_{d+\varepsilon^{-1}\to\infty}\,\frac{\ln n^{\rm
std}(\varepsilon,F_{d})}{d^{\,\alpha}+(1+\ln\varepsilon^{-1})^{\beta}}=0$

for all $\alpha,\beta>0$, i.e., $\widetilde{\mathrm{APP}}$ is exponentially
uniformly weakly tractable for the class $\Lambda^{\mathrm{std}}$.
∎

We finally discuss EXP-UWT and EXP-WT.

######  Theorem 7.

Assume that $F_{d}\in\mathcal{A}$ for every $d\in\mathbb{N}$. If the problem
$\widetilde{\mathrm{APP}}$ is exponentially (uniformly) weakly tractable for
the class $\Lambda^{\mathrm{all}}$, then it is exponentially (uniformly)
weakly tractable for the class $\Lambda^{\mathrm{std}}$.

###### Proof.

Assume that there are $0<\alpha,\beta\leq 1$ such that

$\lim_{d+\varepsilon^{-1}\to\infty}\,\frac{\ln n^{\rm
all}(\varepsilon,F_{d})}{d^{\alpha}+(1+\ln\varepsilon^{-1})^{\beta}}=0.$

It is enough to show that

(3) $\lim_{d+\varepsilon^{-1}\to\infty}\,\frac{\ln n^{\rm
std}(\varepsilon,F_{d})}{d^{\alpha}+(1+\ln\varepsilon^{-1})^{\beta}}=0.$

By assumption, for every $0<h\leq 1/16$, there is some $v_{0}\in\mathbb{N}$
such that

$0\leq\frac{\ln n^{\rm
all}(\varepsilon,F_{d})}{d^{\alpha}+(1+\ln\varepsilon^{-1})^{\beta}}\leq h$

for all $\varepsilon\in(0,1)$ and $d\in\mathbb{N}$ with
$d^{\alpha}+(1+\ln\varepsilon^{-1})^{\beta}\geq v_{0}$. It follows that

$c_{n}(F_{d},L_{2})\leq e\exp\left(-\left(\frac{\ln
n}{h}-d^{\alpha}\right)^{1/\beta}\right)$

for all $n\geq\exp(hv_{0})$. From Theorem 2, we get

$a_{n}(F_{d},L_{2})\leq 2e\exp\left(-\left(\left(\frac{\ln
n}{h}-d^{\alpha}\right)^{1/\beta}-\frac{1}{2}\ln n\right)\right).$

For all $n\geq\exp(2hd^{\alpha})$ and $h\leq 1/16$, we have

$\left(\frac{\ln n}{h}-d^{\alpha}\right)^{1/\beta}-\frac{1}{2}\ln
n\geq\frac{1}{2}\left(\frac{\ln
n}{h}-d^{\alpha}\right)^{1/\beta}+\frac{1}{8}\frac{\ln n}{h}$

and hence we have for all $n\geq\max\\{\exp(hv_{0}),\exp(2hd^{\alpha})\\}$
that

$a_{n}(F_{d},L_{2})\leq 2e\exp\left(-\frac{1}{2}\left(\frac{\ln
n}{h}-d^{\alpha}\right)^{1/\beta}\right)\cdot n^{-1/(8h)}.$

It follows from Theorem 1 that for some absolute constant $b\in\mathbb{N}$ and
all $n\geq\max\\{\exp(hv_{0}),\exp(2hd^{\alpha})\\}$, we have

$\displaystyle e_{bn}(F_{d},L_{2})$ $\displaystyle\leq\frac{1}{n}\sum_{k\geq
n}a_{k}(F_{d},L_{2})\leq\sum_{k\geq n}a_{k}(F_{d},L_{2})$ $\displaystyle\leq
2e\exp\left(-\frac{1}{2}\left(\frac{\ln
n}{h}-d^{\alpha}\right)^{1/\beta}\right)\sum_{k\geq n}k^{-1/(8h)}$
$\displaystyle\leq 2e\exp\left(-\frac{1}{2}\left(\frac{\ln
n}{h}-d^{\alpha}\right)^{1/\beta}\right),$

where we again used that $h\leq 1/16$. It follows that

$\displaystyle n^{\rm std}(\varepsilon,F_{d})\leq
D\exp\left(4h\left((1+\ln\varepsilon^{-1})^{\beta}+d^{\alpha}\right)\right)$

for some absolute constant $D>0$ and all $\varepsilon\in(0,1)$ and
$d\in\mathbb{N}$ such that $d+\varepsilon^{-1}$ is sufficiently large. This
implies

$0\leq\lim_{d+\varepsilon^{-1}\to\infty}\,\frac{\ln n^{\rm
std}(\varepsilon,F_{d})}{d^{\alpha}+(1+\ln\varepsilon^{-1})^{\beta}}\leq 4h.$

Since $h\in(0,1/16)$ can be chosen arbitrarily close to $0$, we obtain (3).
This allows us to conclude our statement. Indeed, for uniform weak
tractability we take arbitrary $\alpha$ and $\beta$ from $(0,1)$, and for weak
tractability we take $\alpha=\beta=1$. ∎

##  4\. Appendix: technical lemmas

The following lemmas are used in the proofs of our results.

######  Lemma 8.

Let $A$ and $B$ be arbitrary positive real numbers. For $n\geq A(B/2)^{B}$ we
have the following inequality

$\displaystyle\sum_{k\geq
n+1}k^{1/2}\exp(-(k/A)^{1/B})\leq\int_{n}^{\infty}t^{1/2}\exp(-(t/A)^{1/B})dt.$

###### Proof.

It is enough to show that the function $f:(0,\infty)\to\mathbb{R}$ given by
$f(t)=t^{1/2}\exp(-(t/A)^{1/B})$ is decreasing on $(A(B/2)^{B},\infty)$.
Indeed, for $t>A(B/2)^{B}$ we have

$\displaystyle f^{\prime}(t)$
$\displaystyle=\left(\exp\left(\frac{1}{2}\ln(t)-(t/A)^{1/B}\right)\right)^{\prime}=$
$\displaystyle=\exp\left(\frac{1}{2}\ln(t)-(t/A)^{1/B}\right)\left(\frac{1}{2t}-(1/AB)(t/A)^{1/B-1}\right)<0.$

∎

######  Lemma 9.

Let $A$ and $B$ be arbitrary positive real numbers. For every $n\geq
A\,\max(3B/2,1)^{B}$ we have the following inequality

$\displaystyle\int_{n}^{\infty}t^{1/2}\exp(-(t/A)^{1/B})dt\,\leq\,A^{1/B}\,B\,\max(3B/2,1)\,n^{3/2-1/B}\,\exp(-(n/A)^{1/B}).$

###### Proof.

Using integration by substitution, with $u=(t/A)^{1/B}$, we obtain that

$\int_{n}^{\infty}t^{1/2}\exp(-(t/A)^{1/B})dt\,=\,A^{3/2}\,B\,\Gamma(3B/2,(n/A)^{1/B})$

where, for $a\in\mathbb{R}$ and $x>0$,
$\Gamma(a,x)=\int_{x}^{\infty}v^{a-1}\exp(-v)dv$ is the incomplete gamma
function.

It is known (see, e.g., Satz 4.4.3 in [13]) that for $a\geq 1$ and $x>a$ we
have

$\Gamma(a,x)\leq a\,x^{a-1}\,\exp(-x).$

If, on the other hand, $0<a<1$ and $x>1$ then since $v^{a-1}\leq x^{a-1}$ for
$v\geq x$ we have

$\displaystyle\Gamma(a,x)$
$\displaystyle=\int_{x}^{\infty}v^{a-1}\exp(-v)dv\leq
x^{a-1}\int_{x}^{\infty}\exp(-v)dv=x^{a-1}\exp(-x).$

Therefore, for every $a>0$ and $x>\max(a,1)$, the following bound holds

$\Gamma(a,x)\,\leq\,\max(a,1)\,x^{a-1}\,\exp(-x).$

Thus for $n>A\,\max(3B/2,1)^{B}$, and taking $a=3B/2$ and $x=(n/A)^{1/B}$, we
have

$\displaystyle\int_{n}^{\infty}t^{1/2}\exp(-(t/A)^{1/B})dt\,\leq\,A^{1/B}\,B\,\max(3B/2,1)\,n^{3/2-1/B}\,\exp(-(n/A)^{1/B}).$

∎

######  Lemma 10.

For every $A,B,n,\delta,u>0$ we have the following inequality

$\displaystyle
n^{u}\,\exp\left(-(n/A)^{1/B}\right)\,\leq\,A^{u}\,\delta^{-uB}\,\exp\left((uB\delta-1)(n/A)^{1/B}\right).$

###### Proof.

Let $x=\delta(n/A)^{1/B}$. Then $n^{u}=A^{u}\,\delta^{-uB}\,x^{uB}$. Using the
fact that $\ln(x)\leq x$ for all $x>0$ we obtain that

$\displaystyle\ln(n^{u})$ $\displaystyle=\ln(A^{u}\delta^{-uB})+uB\ln(x)\leq$
$\displaystyle\leq\,\ln(A^{u}\delta^{-uB})+uBx\,=\,\ln(A^{u}\delta^{-uB})+uB\delta(n/A)^{1/B}$

Hence, taking exponentials of both sides we derive that

$\displaystyle n^{u}\leq A^{u}\delta^{-uB}\exp(uB\delta(n/A)^{1/B})$

and thus

$\displaystyle n^{u}\exp(-(n/A)^{1/B})\leq
A^{u}\delta^{-uB}\exp((uB\delta-1)(n/A)^{1/B})$

as claimed. ∎

Funding David Krieg is supported by the Austrian Science Fund (FWF) Project
F5506, which is part of the Special Research Program “Quasi-Monte Carlo
Methods: Theory and Applications”.

Declarations

Conflict of interest The authors declare no competing interests.

## References

  * [1] B. Adcock, A. C. Hansen, and A. Shadrin.  A stability barrier for reconstructions from Fourier samples.  SIAM Journal on Numerical Analysis, 52(1):125–139, 2014.
  * [2] B. Adcock, D. Huybrechs, and J. Martin-Vaquero.  On the numerical stability of Fourier extensions.  Foundations of Computational Mathematics, 14:635–687, 2014.
  * [3] B. Adcock, R. B. Platte, and A. Shadrin.  Optimal sampling rates for approximating analytic functions from pointwise samples.  IMA Journal of Numerical Analysis, 39(3):1360–1390, 2018.
  * [4] P. Binev, A. Cohen, W. Dahmen, R. DeVore, G. Petrova, and P. Wojtaszczyk.  Convergence rates for greedy algorithms in reduced basis methods.  SIAM Journal on Mathematical Analysis, 43(3):1457–1472, 2011.
  * [5] A. Buffa, Y. Maday, A. T. Patera, C. Prud’homme, and G. Turinici.  A priori convergence of the greedy algorithm for the parametrized reduced basis method.  ESAIM: Mathematical Modelling and Numerical Analysis – Modélisation Mathématique et Analyse Numérique, 46(3):595–603, 2012.
  * [6] J. Creutzig and P. Wojtaszczyk.  Linear vs. nonlinear algorithms for linear problems.  Journal of Complexity, 20:807–820, 2004.
  * [7] J.J. A. Deimer and A. T. Sergio.  Estimates for n-widths of sets of smooth functions on complex spheres.  Journal of Complexity, 64:101537, 2021.
  * [8] J. Dick, P. Kritzer, F. Pillichshammer, and H. Woźniakowski.  Approximation of analytic functions in Korobov spaces.  Journal of Complexity, 30(2):2–28, 2014.
  * [9] M. Dolbeault, D. Krieg, and M. Ullrich.  A sharp upper bound for sampling numbers in $L_{2}$.  Applied and Computational Harmonic Analysis, 63:113–134, 2023.
  * [10] D. Dũng, N. Temlyakov V. and T. Ullrich.  Hyperbolic Cross Approximation.  Advanced Courses in Mathematics - CRM Barcelona, Springer International Publishing, 2018.
  * [11] A. Ebert and F. Pillichshammer.  Tractability of approximation in the weighted Korobov space in the worst-case setting – a complete picture.  Journal of Complexity, 67:101571, 2021.
  * [12] G. E. Fasshauer, F. J. Hickernell, and H. Woźniakowski.  On dimension-independent rates of convergence for function approximation with Gaussian kernels.  SIAM Journal on Numerical Analysis, 50(1):247–271, 2012.
  * [13] W. Gabcke.  Neue Herleitung und explizite Restabschätzung der Riemann-Siegel-Formel.  PhD thesis, University of Göttingen, 1979.
  * [14] M. Griebel and J. Oettershagen.  On tensor product approximation of analytic functions.  Journal of Approximation Theory, 207:348–379, 2016.
  * [15] B. Haasdonk.  Convergence rates of the pod–greedy method.  ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 47(3):859–873, 2013.
  * [16] W. Hackbusch and B. N. Khoromskij.  Tensor-product approximation to operators and functions in high dimensions.  Journal of Complexity, 23(4):697–714, 2007.  Festschrift for the 60th Birthday of Henryk Woźniakowski.
  * [17] T. Hangelbroek and A. Ron.  Nonlinear approximation using Gaussian kernels.  Journal of Functional Analysis, 259(1):203–219, 2010.
  * [18] A. Hinrichs, D. Krieg, E. Novak, J. Prochno, and M. Ullrich.  Random sections of ellipsoids and the power of random information.  Transactions of the American Mathematical Society, 374(12):8691 – 8713, 2021.
  * [19] A. Hinrichs, D. Krieg, E. Novak, and J. Vybíral.  Lower bounds for the error of quadrature formulas for Hilbert spaces.  Journal of Complexity, 65:101544, 2021.
  * [20] A. Hinrichs, E. Novak, and M. Ullrich.  On weak tractability of the Clenshaw–Curtis Smolyak algorithm.  Journal of Approximation Theory, 183:31–44, 2014.
  * [21] A. Hinrichs, E. Novak, M. Ullrich, and H. Woźniakowski.  The curse of dimensionality for numerical integration of smooth functions.  Mathematics of Computation, 83(290):2853–2863, 2014.
  * [22] A. Hinrichs, E. Novak, M. Ullrich, and H. Woźniakowski.  Product rules are optimal for numerical integration in classical smoothness spaces.  Journal of Complexity, 38:39–49, 2017.
  * [23] A. Hinrichs, J. Prochno, and M. Ullrich.  The curse of dimensionality for numerical integration on general domains.  Journal of Complexity, 50:25–42, 2019.
  * [24] C. Irrgeher, P. Kritzer, F. Pillichshammer, and H. Woźniakowski.  Approximation in Hermite spaces of smooth functions.  Journal of Approximation Theory, 207:98 – 126, 2016.
  * [25] C. Irrgeher, P. Kritzer, F. Pillichshammer, and H. Woźniakowski.  Tractability of multivariate approximation defined over Hilbert spaces with exponential weights.  Journal of Approximation Theory, 207:301–338, 2016.
  * [26] T. Karvonen and S. Särkkä.  Worst-case optimal approximation with increasingly flat Gaussian kernels.  Advances in Computational Mathematics, 46:1–17, 2020.
  * [27] D. Krieg.  Uniform recovery of high-dimensional $C^{r}$-functions.  Journal of Complexity, 50:116–126, 2019.
  * [28] D. Krieg and M. Ullrich.  Function values are enough for $L_{2}$-approximation.  Foundations of Computational Mathematics, 21(4):1141–1151, 2021\.
  * [29] D. Krieg and M. Ullrich.  Function values are enough for $L_{2}$-approximation: Part II.  Journal of Complexity, 66:101569, 2021.
  * [30] P. Kritzer and H. Woźniakowski.  Simple characterizations of exponential tractability for linear multivariate problems.  Journal of Complexity, 51:110–128, 2019.
  * [31] P. Mathé.  _s_ -numbers in information-based complexity.  Journal of Complexity, 6(1):41–66, 1990.
  * [32] N. Nagel, M. Schäfer, and T. Ullrich.  A new upper bound for sampling numbers.  CoRR, abs/2010.00327, 2020.
  * [33] E. Novak and H. Woźniakowski.  Tractability of multivariate problems. Vol. 1: Linear information, volume 6 of EMS Tracts in Mathematics.  European Mathematical Society (EMS), Zürich, 2008.
  * [34] E. Novak and H. Woźniakowski.  Approximation of infinitely differentiable multivariate functions is intractable.  Journal of Complexity, 25:398–404, 2009.
  * [35] E. Novak and H. Woźniakowski.  Tractability of multivariate problems. Volume II: Standard information for functionals, volume 12 of EMS Tracts in Mathematics.  European Mathematical Society (EMS), Zürich, 2010.
  * [36] E. Novak and H. Woźniakowski.  Tractability of multivariate problems. Volume III: Standard information for operators, volume 18 of EMS Tracts in Mathematics.  European Mathematical Society (EMS), Zürich, 2012.
  * [37] E. Novak and H. Woźniakowski.  Tractability of multivariate problems for standard and linear information in the worst case setting: part i.  Journal of Approximation Theory, 207:177–192, 2016.
  * [38] A. Pietsch.  Operator ideals.  Deutscher Verlag, Wiss., Berlin, 1978.
  * [39] Rodrigo B. Platte, Lloyd N. Trefethen, and Arno B. J. Kuijlaars.  Impossibility of fast stable approximation of analytic functions from equispaced samples.  SIAM Review, 53(2):308–318, 2011.
  * [40] J. Vybíral.  A variant of Schur’s product theorem and its applications.  Advances in Mathematics, 368:107140, 2020.
  * [41] I. Sloan and H. Woźniakowski.  Multivariate approximation for analytic functions with Gaussian kernels.  Journal of Complexity, 45:1–21, 2018.
  * [42] V. N. Temlyakov.  Approximation of periodic functions.  Computational Mathematics and Analysis Series, Nova Science Publishers, Inc., Commack, NY, 1993.
  * [43] V. N. Temlyakov.  Multivariate Approximation, volume 32 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, 2018.
  * [44] M. Ullrich.  On the worst-case error of least squares algorithms for $L_{2}$-approximation with high probability.  Journal of Complexity, 60, 2020.
  * [45] H. Wendland.  Scattered data approximation, volume 17 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2005.
  * [46] G. Xu.  On weak tractability of the Smolyak algorithm for approximation problems.  Journal of Approximation Theory, 192:347–361, 2015.
  * [47] J. Zhang.  A note on EC-tractability of multivariate approximation in weighted Korobov spaces for the standard information class.  Journal of Complexity, 67, 2021.


# Discussion of Martingale Posterior Distributions

David Rossell

Here is a brief discussion of Martingale posterior distributions, by Fong et
al. (2021). Congratulations on a thought-provoking piece. Building Bayesian
inference from a (likelihood, predictive) pair, rather than a (likelihood,
prior), enriches the paradigm and provides new ways to think about, formulate
and solve problems. A few respectful remarks.

First, although the authors never claim this, it is worth emphasizing that the
framework is not prior-free. There is a posterior and a likelihood, hence the
prior is proportional to their ratio. The key is that said prior is data-
dependent, providing an interesting avenue to develop objective Bayes methods,
at the cost of loosing the coherence property in belief updating. Inspecting
the prior can be informative. Figure 1 shows a Bernoulli example where truly
$\theta=0.5$ but the implied prior places little mass around that value, and a
Gaussian example where the prior is centered around the sample mean111code at
https://github.com/davidrusi/paper_examples/tree/main/2022_Rossell_martingale_posteriors).
This apparently erratic prior behavior might be problematic for model choice
via Bayes factors, e.g. returning a very small integrated likelihood in the
Bernoulli example.

Second, while sometimes it is easier to elicit a predictive than a prior, in
my experience the reverse is often true. For example in regression a prior on
parameters defines a prior on the $R^{2}$ coefficient, an easy-to-interpret
quantity, whereas eliciting predictives may be less intuitive for non-
statisticians. Further, note that computational considerations elegantly
discussed by the authors severely restrict the range of predictives one may
consider in practice, limiting the flexibility of the framework.

Third, I am afraid I disagree on the frameworks’ computational convenience.
Doing a single optimization may be faster than sampling, but the framework
requires solving many optimizations. This is not cheaper than posterior
sampling in a standard (likelihood, prior) construction, also the latter
offers fast non-sampling based tools, e.g. Laplace approximations and
extensions. It would be interesting to consider analogues for the predictive
framework.

Finally, a remark on assuming that at $n=\infty$ there is no uncertainty left.
In some settings this is not true, e.g. in high-dimensional regression with
$p\gg n$ (one adds higher-order polynomial terms as $n$ grows, say) and a
normal prior on the parameters there remains posterior uncertainty even as
$n\rightarrow\infty$. The proposed framework does not account for such
uncertainty, unless suitable adjustments are made.

Bernoulli | Gaussian
---|---
|
Figure 1. Likelihood, posterior and prior densities for Bernoulli(0.5) and
Normal(5,1) examples with $n=100$

## References

  * Fong et al. [2021] Edwin Fong, Chris Holmes, and Stephen G Walker.  Martingale posterior distributions.  _arXiv_ , 2103.15671:1–62, 2021.