[en-GB]ord=omit 2.5cm2.5cm2.5cm2.5cm

# Statistical visualisation of tidy and geospatial data in R via kernel
smoothing methods in the eks package

Tarn Duong111Paris, France F-75000. Email<EMAIL_ADDRESS>
###### Abstract

Kernel smoothers are essential tools for data analysis due to their ability to
convey complex statistical information with concise graphical visualisations.
Their inclusion in the base distribution and in the many user-contributed add-
on packages of the R statistical analysis environment caters well to many
practitioners. Though there remain some important gaps for specialised data
types, most notably for tibbles (tidy data) within the tidyverse, and for
simple features (geospatial data) within geospatial analysis. The proposed eks
package fills in these gaps. In addition to kernel density estimation, this
package also caters for more complex data analysis situations, such as density
derivative estimation, density-based classification (supervised learning) and
mean shift clustering (unsupervised learning). We illustrate with experimental
data how to obtain and to interpret the statistical visualisations for these
kernel smoothing methods.

Keywords: classification, clustering, ggplot2, GIS, kernel density estimation,
sf

##  1 Introduction

Kernel smoothers form an essential suite of statistical techniques for data
analysis in the 21st century due to their ability to convey complex
statistical information in a concise and intuitive visual format. This ability
arises from their shared characteristic of transforming data samples into
smoothed estimates. Kernel smoothers have provided insight in data analysis
problems in many situations. A small recent selection of these includes: the
identification of important biomedical functions, such as characterising
different sub-cellular structures in single cells (Schauer et al.,, 2010) or
characterising a single cell population in mixed cell samples (Chacón et al.,,
2011); the evaluation of predicted extreme temperatures to calibrate climate
models (Béranger et al.,, 2019); the estimation of the home range of animal
movements (Baíllo and Chacón,, 2021); or the detection of traffic anomalies
from traffic flows (Kalair and Connaughton,, 2021).

A major access point to kernel smoothers in the R statistical programming
environment is the ks (‘kernel smoothing’) add-on package (Duong,, 2007),
which implements density estimation, density derivative estimation,
classification (unsupervised learning), clustering (unsupervised learning),
and inferential methods. This package utilises the base R graphics engine to
generate its statistical graphics. Whilst it remains the most comprehensive
graphics engine in R, the ggplot2 graphics engine (Wickham,, 2016) has gained
popularity, as part of the ‘tidyverse’, especially with data analysis
practitioners. Despite the dramatic rise in the number of analysis methods
available in the tidyverse, nonetheless it comprises a limited range of
natively implemented kernel smoothers. The first goal of the eks (‘extended
kernel smoothing’) package (Duong,, 2023) is to provide access to a
comprehensive suite of kernel smoothers in the tidyverse.

There is an analogous lack of kernel smoothers for geospatial data analysis.
Since the term ‘geospatial’ data analysis refers to many different yet
overlapping concepts, we employ it in this paper to refer to data analysis
which is compatible with ‘Geographical Information Systems’ (GIS). Within R,
the sf package (Pebesma,, 2018) provides geospatial/GIS functionality via its
robust implementation of the ‘simple features’ GIS standard data format (OGC,,
2021). The eks package relies on this simple features implementation, in order
to take advantage of the ability to facilitate visualisations in both ggplot2
and base R graphical engines, and input/output to external GIS software (such
as ArcGIS and QGIS). The second goal of the eks package is to provide access
to a comprehensive suite of kernel smoothers for geospatial analysis.

Thus a wide range of kernel smoothers is now available for tidy and geospatial
data, and for ggplot2 and base R graphical visualisations. The user is able to
select and combine these components, with their differing strengths and
applicabilities, in order to construct suitable data analysis workflows. To
illustrate kernel smoothers, we employ the tidy data set air from the ks
package, and the geospatial data set grevilleasf from the eks package, as
shown in Figure 1.

|
---|---
Figure 1: Scatter plots of data sets. (Left) Tidy air quality measurements air
($n=30239$), in the Châtelet underground train station in Paris, France.
(Right) Geospatial Grevillea locations grevilleasf ($n=22203$), in the
biodiversity hotspot of south-western Western Australia.

A tidy data set is a data matrix where (i) each variable forms a column and
(ii) each observation forms a row, and it is also known as a ‘long’ data set
(Wickham,, 2014). The first three records of the air data set are

    R> data(air, package="ks")
    R> air
            date  time no no2 pm10 co2 temp humi
    1 2013-01-01 01:00  6  31  182 776 21.5 46.2
    2 2013-01-01 02:00  6  30  166 800 21.6 47.3
    3 2013-01-01 03:00  4  27  124 799 21.8 47.0
    ...

These are the hourly mean air quality measurements from 01 January 2013 to 31
December 2016 in the Châtelet underground train station, which is a major hub
in the Paris transport network (RATP,, 2016). We focus on the concentrations
of carbon dioxide CO2 (g/m3) (co2) and of particulate matter less than 10
$\mu$m in diameter PM10 (parts per million) (pm10), and the hourly interval
(time). The concentrations of CO2 indicate the renewal rate of fresh air, and
of PM10 the potential to affect adversely respiratory health. There are
$n=30239$ complete (co2, pm10) measurements.

The geospatial data grevilleasf data set consists of 22203 plants from 238
different Grevillea species collected in Western Australia. The south-west
corner of Western Australia is one of the 25 ‘biodiversity hotspots’ which are
‘areas featuring exceptional concentrations of endemic species and
experiencing exceptional loss of habitat’ identified in Myers et al., (2000)
to assist in formulating priorities in biodiversity conservation policies. The
geodetic coordinates (degrees) of the Grevillea locations into planar
coordinates (metres) using the GDA2020/MGA zone 50 (EPSG:7850) projection.
They are encoded as a simple feature in the column geometry, which is a
special data structure that cannot be treated like the usual floating point
variables in data frames or tibbles, and require specialised methods
implemented by the sf package. So we refer to grevilleasf solely as a
geospatial data set, and omit any mention of its tidy status, to emphasise its
distinct geospatial characteristics.

    R> data(grevilleasf, package="eks")
    R> grevilleasf
    Simple feature collection with 22303 features and 2 fields
    Geometry type: POINT
    Dimension:     XY
    Bounding box:  xmin: 73519.97 ymin: 6120859 xmax: 1795868 ymax: 8451928
    Projected CRS: GDA2020 / MGA zone 50
                      name   species                 geometry
    1    Grevillea robusta   robusta POINT (390106.5 6462671)
    2   Grevillea speciosa  speciosa POINT (382689.2 6457387)
    3    Grevillea robusta   robusta POINT (390089.8 6462603)
    ...

This paper focuses on the software implementation of the kernel smoothers, and
is complementary to Chacón and Duong, (2018) which focuses on the underlying
statistical framework. In Section 2, we explore kernel density estimation, in
Section 3 classification (supervised learning), in Section 4 density gradient
estimation, and in Section 5 clustering (unsupervised learning). We illustrate
each case first for tidy data with ggplot2 graphics, and then for geospatial
data with ggplot2 and base R graphics. In Section 6, we briefly mention kernel
smoothers in other data analysis settings, which are implemented in the eks
package but have been omitted for brevity, and we end with some concluding
remarks.

##  2 Density estimation

Density estimation is a fundamental statistical analysis tool, since it
supplies much information about the data set at hand. Our data ${\bf
X}_{1},\dots,{\bf X}_{n}$ is a random sample drawn from the common density
function $f$. The goal of density estimation, as its name suggests, is to
estimate this unknown density. Kernel density estimates are a popular choice
among the many available smoothed density estimation methods, since they
posses an intuitive construction. It is most widely used kernel smoother, and
can be considered to be a smoothed version of the histogram. For an arbitrary
estimation point ${\bm{x}}$, the kernel density estimate is

$\hat{f}_{\bf H}({\bm{x}})=n^{-1}\sum_{i=1}^{n}K_{\bf H}({\bm{x}}-{\bf
X}_{i}).$ (1)

Throughout the eks package, the kernel function is the $2$-dimensional
Gaussian density function $K_{\bf H}({\bm{x}})=(2\pi)^{-1}|{\bf
H}|^{-1/2}\exp(-\tfrac{1}{2}{\bm{x}}^{\top}{\bf H}^{-1}{\bm{x}})$. Equation
(1) tells us that to compute a kernel density estimate, we place a Gaussian
function, with variance ${\bf H}$, at each data point ${\bf X}_{i}$, and then
we sum these kernel functions. This way, the data sample ${\bf
X}_{1},\dots,{\bf X}_{n}$ are transformed into a smooth surface $\hat{f}_{\bf
H}$. Chacón and Duong, (2018, Chapter 2) contains a more detailed overview of
kernel density estimates.

The bandwidth matrix ${\bf H}$ in Equation (1) is the crucial tuning
parameter. A bandwidth matrix which is too small leads to an undersmoothed
density estimate since it does not offer sufficient reduction in the
complexity of the observed data. On the other hand, a bandwidth matrix which
is too large leads to an oversmoothed density estimate that obscures important
details in the observed data. Thus it is critical to find an optimal trade-off
between this under- and oversmoothing. Many possible solutions for optimal
smoothing are implemented in the ks package, and are thus available in the eks
package, including the plug-in, unbiased cross validation and smoothed cross
validation bandwidths.

###  2.1 Tidy density estimation

To illustrate density estimation for tidy data, we focus on a single hourly
interval of the air quality measurements. Figure 2 compares the density
estimates, for the $n=1285$ measurements from 11:00 to 12:00, with an optimal
bandwidth and a sub-optimal one. The optimal bandwidth is computed from the
eks package and the sub-optimal one from the ggalt package (Rudis et al.,,
2017). The former, known as the bivariate plug-in bandwidth matrix (Duong and
Hazelton,, 2003), is the default optimal bandwidth in the eks package, and it
is obtained from a call to ks::Hpi. For the air quality measurements, the
optimal ks::Hpi matrix is [342.1, 97.2; 97.2 365.2]. The presence of non-zero
off-diagonal entries in the optimal matrix appropriately orients the kernel
functions, and the resulting density estimate is unimodal, as shown in the
centre panel of Figure 2. The default bandwidth in ggalt, which is widely used
in the tidyverse, is obtained from the element-wise application of the
univariate plug-in bandwidth KernSmooth::dpik. For the air quality
measurements, this bandwidth is [218.6, 0; 0, 124.9]. Since this sub-optimal
matrix only applies smoothing in the coordinate axis directions, it yields an
undersmoothed density estimate with spurious bimodal structure on the right
panel.

Figure 2: Filled contour plots of density estimates for air quality
measurements 11:00–12:00 ($n=1285$) with quartile probability contour levels.
(Left) Scatter plot. (Centre) Optimally smoothed. (Right) Undersmoothed.

In Figure 2, the heights of the contour regions are calculated according to
the probability contours method (Bowman and Foster,, 1993, Hyndman,, 1996).
The pink region is the smallest region that contains 25% of the probability
mass, the orange region plus the enclosed pink region is the smallest region
that contains 50% of the probability mass, and the yellow region plus the
enclosed orange and pink regions is the smallest region that contains 75% of
the probability mass. Since these are relative heights, they facilitate the
choice of the contour levels, since it involves selecting values from 0% to
100%, rather than from the range of the density values. These probability
contours can also be considered as a multivariate extension of the univariate
percentiles, e.g., the 50% contour region is a bivariate equivalent to the
median. Due to their intuitive properties, these probability contours are
employed throughout in eks, with the quartile contour levels (25%, 50%, 75%)
being the default values. In addition to their intuitive interpretation, these
probability contours are straightforward to compute: the kernel density
estimate is evaluated at the $n$ observed data values $\hat{f}_{\bf H}({\bf
X}_{1}),\dots,\hat{f}_{\bf H}({\bf X}_{n})$, then we compute $\tau_{\alpha}$
as the $\alpha$-quantile of these evaluated values, and the $\alpha$
probability contour region is the level set of the density estimate at
$\tau_{\alpha}$, i.e. $\\{{\bm{x}}:\hat{f}_{\bf H}({\bm{x}})>\tau_{\alpha}\\}$
(Hyndman,, 1996). These probability contours are also implemented in the
ggdensity package (Otto and Kahle,, 2023), though with a similar sub-optimal
bandwidths as in ggalt.

The R code snippets included here are intended to give an overall idea of the
syntax of the eks package, rather than a complete code to reproduce the
figures. The latter is provided in the companion R script. The code snippet to
compute the density estimate with the optimal bandwidth, in the centre panel
in Figure 2, is

    R> ## tidy density estimate
    R> air2 <- ungroup(filter(air, time=="11:00"))
    R> air2 <- select(air2, co2, pm10)
    R> t1 <- tidy_kde(air2, H=H1)
    R> ggplot(t1, aes(x=co2, y=pm10)) + geom_contour_filled_ks(colour=1)

The function tidy_kde is a wrapper function for ks::kde, which computes the
tidy density estimate explicitly. This differs from existing layer functions,
e.g., ggplot2::geom_density_2d and ggalt::geom_bkde2d, which compute the
density estimate internally and do not return a user-level R object. The tidy
density estimate from tidy_kde is:

    R> t1
    # A tibble: 22,801  6
         co2  pm10  estimate ks        tks   label
       <dbl> <dbl>     <dbl> <list>    <chr> <chr>
     1  342. -28.8 -5.65e-24 <kde>     kde   Density
     2  346. -28.8 -7.42e-22 <int [1]> kde   Density
     3  350. -28.8  2.38e-22 <int [1]> kde   Density
    ...

This output is a tibble with an added tidy_ks class, which allows for a
ggplot.tidy_ks method to be defined for this object class. Otherwise, it can
be treated as a tibble. The first two columns co2, pm10 (same names as the
input data) are the coordinates of the vertices in the estimation grid, the
third column estimate is the density estimate value at co2, pm10. The fourth
column ks holds the output from ks::kde. This is required for the computation
of probability contours in the new layer function geom_contour_filled_ks to
draw the filled contour plots for tidy_ks objects. The remaining columns
indicate that the output is a density estimate computed from ks::kde, and they
are employed in ggplot.tidy_ks to create default aesthetic mapping and legend
labels. This default aesthetic mapping is ggplot2::aes(x=co2, y=pm10,
z=estimate, weight=ks). Whilst the x, y, z aesthetics are as expected for a
bivariate contour plot, the weight aesthetic is unorthodox, since it is not a
weighting variable: it is a workaround in ggplot2 graphics to mimic the
dynamic display of probability contours in base R graphics.

For the air quality measurements from 11:00 to 12:00, the quartile contour
levels for the optimally smoothed density estimate in the centre panel in
Figure 2 are 3.31e-5, 2.42e-5, 1.22e-5, and for the undersmoothed density
estimate in the right panel are 3.51e-5, 1.22e-5, 2.53e-5. These probability
contour heights are different for each different density estimate, even if the
target contour probabilities remain the same at 25%, 50%, 75%. On the other
hand, it is some times useful to have a set of fixed contour heights for all
density estimates for a direct comparison. A heuristic method consists of
computing the probability contour heights for each density estimate, for a
fixed set of probabilities, which are then aggregated. We compute the
corresponding probabilities for each density estimate for this aggregated set
of contour heights, and we remove any contour levels whose estimated
probability which are too close to each other. This procedure is implemented
in the contour_breaks function, though some trial and error is still likely
required to produce visually appealing contour plots for all density estimates
(Chacón and Duong,, 2018, Section 2.2).

We revisit the density estimates for the air quality measurements from 11:00
to 12:00, this time with the fixed contour heights (3.90e-6, 1.46e-5, 2.48e-5,
3.25e-5, 4.15e-5) in Figure 3. With these fixed contour heights, a direct
comparison of different density estimates is possible. The optimally smoothed
estimate on the right exceeds the two highest contour heights (4.15e-5 dark
pink, 3.25e-5 dark orange) with a unimodal bump, whereas the density estimate
on right exceeds them with a bimodal structure. Since the latter adds some
spurious modal information, it is considered to be undersmoothed in
comparison.

Figure 3: Filled contour plots of density estimates for air quality
measurements from 11:00 to 12:00 ($n=1285$) with fixed contour levels. (Left)
Optimally smoothed. (Right) Undersmoothed.

The code to produce two density estimates with a single set of contour heights
in Figure 3 is

    R> ## fixed contour levels
    R> H2 <- diag(sapply(air2, KernSmooth::dpik)^2)
    R> t2 <- tidy_kde(air2, H=H2)
    R> t3 <- c(t1, t2)
    R> b <- contour_breaks(t3, cont=c(10,30,50,70,90))
    R> ggplot(t3, aes(x=co2, y=pm10)) + geom_contour_filled_ks(colour=1, breaks=b)
    +  facet_wrap(~group)

###  2.2 Geospatial density estimation

To illustrate density estimation for geospatial data, we utilise single
species subsets of the Grevillea locations. Figure 4 compares the density
estimates for the $n=93$ locations of the G. yorkrakinensis species which
result from an optimal plug-in bandwidth [8.84e8, $-8.33$e8; $-8.33$e8,
1.36e9] and a sub-optimal bandwidth [5.43e8, 0; 0, 9.10e8]. The optimally
smoothed density estimate in the centre panel in Figure 4 displays a trimodal
structure with obliquely oriented contours. The oversmoothed density estimate
in the right panel, whilst it also has trimodal structure, it has circular
contours which do not follow closely the orientation of the observed data
points, and so is considered to be oversmoothed.

Figure 4: Filled contour plots of density estimates for G. yorkrakinensis
($n=93$) with quartile probability contour levels. (Left) Scatter plot.
(Centre) Optimally smoothed. (Right) Oversmoothed.

To produce the centre panel in Figure 4 for the G. yorkrakinensis locations,
the commands are

    R> ## geospatial density estimate
    R> yorkr <- filter(grevilleasf, species=="yorkrakinensis")
    R> s1 <- st_kde(yorkr)

The function st_kde is the geospatial equivalent of tidy_kde, and produces an
object of class sf_ks, which is a list of 3 fields: tidy_ks, grid, and sf. The
first field is a summary of the tidy density estimate from tidy_kde, the
second are the rectangular polygons of the estimation grid, and the third are
the 1% to 99% probability contour regions of the density estimate. We focus on
the contour regions.

    R> s1$sf
    Simple feature collection with 99 features and 2 fields
    Geometry type: GEOMETRY
    Dimension:     XY
    Bounding box:  xmin: 429181.6 ymin: 6333015 xmax: 757134 ymax: 6793115
    Projected CRS: GDA2020 / MGA zone 50
       contlabel     estimate                       geometry
    1         99 2.491961e-12 POLYGON ((430965 6755235, 4...
    2         98 3.699348e-12 POLYGON ((437905.1 6746871,...
    3         97 6.714617e-12 POLYGON ((448315.2 6740434,...
    ...

This has 2 attributes: contlabel (label of the probability contour) and
estimate (height of contour probability region). Unlike for tidy_kde where the
probability contour regions are computed dynamically in the layer function
geom_contour_filled_ks, these 1% to 99% regions are converted to
(multi)polygons prior to plotting since the dynamic conversion during plotting
could be computationally heavy. Since we are unable replicate the automatic
selection of the quartile contours 25%, 50%, 75% by default, like in
geom_contour_filled_ks, for the ggplot2::geom_sf layer function, we first
apply st_get_contour to the input of ggplot2::geom_sf. The sf_ks class also
has a ggplot.sf_ks method which computes the default map legend.

    R> ## geospatial density estimate geom_sf plot
    R> ggplot(s1) + geom_sf(data=st_get_contour(s1), aes(fill=contlabel))

The following command produces the equivalent in base R graphics to the
ggplot2 plot in the centre panel in Figure 3.

    R> ## geospatial density estimate base R plot
    R> plot(s1)

This plot.sf_ks method for sf_ks objects method internally calls
st_get_contour to extract the required contour polygons for plotting, so it is
more concise than ggplot2::geom_sf that requires an explicit user-level call
to st_get_contour. The base R and ggplot2 plots are essentially identical
since they comply with the geospatial standard specifications for simple
features.

###  2.3 Optimal bandwidth matrices

Since the bandwidth matrix is the crucial tuning parameter for kernel density
estimates, we explore further their statistical properties. These properties
are the subject of a vast body of research literature, which we do not attempt
to review here, and instead provide a simplified outline of how the optimal
bandwidth matrix in eks is obtained.

We begin with a squared error discrepancy between a density estimate
$\hat{f}_{\bf H}$ and the target density $f$, i.e., $M({\bf
H})=\int\operatorname{\bm{\mathbb{E}}}[\hat{f}_{\bf
H}({\bm{x}})-f({\bm{x}})]^{2}\,\mathsf{d}{\bm{x}}$. Since this expression
involves the unknown target density $f$, it must be estimated for it to be of
practical use. The plug-in bandwidth matrix in ks::Hpi computes the estimate
$\hat{M}({\bf H})=(4\pi)^{-d/2}n^{-1}|{\bf
H}|^{-1/2}+\tfrac{1}{4}\hat{{\bm{m}}}_{4}^{\top}(\operatorname{vec}{\bf
H}\otimes\operatorname{vec}{\bf H})$. We omit to describe this estimate
rigorously since it would require lengthy technical definitions: the
interested reader is encouraged to consult Chacón and Duong, (2018, Chapter 3)
for details. We are content to state that the first term in $\hat{M}$ is
related to the variance of the density estimate, and the second term to the
square of the bias of the density estimate. An optimal bandwidth matrix
$\hat{{\bf H}}$ is defined as

$\hat{{\bf H}}=\underset{{\bf H}}{\mathrm{argmin}}\ \hat{M}({\bf H})$ (2)

where the minimisation is carried out over the space of all symmetric positive
definite matrices. When this minimisation is achieved, then there is an
optimal trade-off between the variance and the squared bias, or equivalently
between over- and under-smoothing. When an optimal bandwidth matrix $\hat{{\bf
H}}$ is substituted into Equation (1), the resulting kernel density estimate
is the closest to the target density $f$ as measured by the discrepancy
$\hat{M}$. Different bandwidth matrices arise from the different ways of
computing $\hat{M}$ and/or from different ways of carrying out the
minimisation. For example, the default bandwidth in ggalt treats the joint
bivariate optimisation in Equation (2) as two separate univariate optimisation
problems. The density estimate functions tidy_kde and st_kde compute
$\hat{{\bf H}}$ in Equation (2) by calling the ks::Hpi function, and then
substitute this $\hat{{\bf H}}$ into Equation (1), to compute an optimal
tidy/geospatial density estimate, as shown in the centre panels in Figures 2
and 3.

Additional bandwidth matrices in the ks package include the normal scale
ks::Hns, unbiased cross validation ks::Hucv and smoothed cross validation
ks::Hscv. The commands are of the type:

    R> ## smoothed cross validation selector
    R> H3 <- ks::Hscv(air2)
    R> t3 <- tidy_kde(air, H=H3)

For most data samples, the plug-in bandwidth ks::Hpi yields fast and robust
kernel estimates, though there remain some cases where other bandwidths are
more suitable. For a review of the performance of these bandwidths, see Chacón
and Duong, (2018, Chapter 3). For brevity, we illustrate kernel estimates only
with the plug-in optimal bandwidths in the sequel.

##  3 Density-based classification (supervised learning)

The goal of classification is to assign future data to one of the known
classes in the current data. That is, this is a supervised learning problem.
The data are $({\bf X}_{1},Y_{1}),\dots,({\bf X}_{n},Y_{n})$, where the ${\bf
X}_{i}$ are the observed attributes, and the $Y_{i}$ is the known class label
from $m$ classes. These data are a random sample from the mixture density
$\pi_{1}f_{1}+\dots+\pi_{m}f_{m}$, where $\pi_{j}$ is the prior probability
and $f_{j}$ is the marginal density function for class $j$, for $j=1,\dots,m$
(Chacón and Duong,, 2018, Section 7.2).

The simplest (and often most robust) classifier is the Bayes classifier. This
assigns a candidate point ${\bm{x}}$ to the class $c$ with the highest density
value at ${\bm{x}}$, i.e.,
$c({\bm{x}})=\mathrm{argmax}_{j=1,\dots,m}\,\pi_{j}f_{j}({\bm{x}})$. The
density-based classifier replaces the prior probability $\pi_{j}$ with the
observed sample class proportion $\hat{\pi}_{j}$, and the marginal density
$f_{j}$ with the marginal density estimate $\hat{f}_{j}$. Each marginal
density estimate is computed with its own optimal bandwidth matrix. The
estimated class label for ${\bm{x}}$ from the kernel density-based classifier
is thus

$\hat{c}({\bm{x}})=\underset{j=1,\dots,m}{\mathrm{argmax}}\
\hat{\pi}_{j}\hat{f}_{j}({\bm{x}}).$

This kernel classifier is more adaptable than the usual linear and quadratic
classifiers. The linear classifier uses Gaussian density fits with a common
variance matrix for all classes, and the quadratic classifier Gaussian density
fits with a different variance matrix for each class.

Our data sample comprises the air quality measurements for three hourly
intervals at six hours apart throughout the day, i.e. 07:00–08:00 (pink
circles, $n_{1}=1282$), 13:00–14:00 (green triangles, $n_{2}=1280$), and
19:00–20:00 (blue squares, $n_{3}=1292$), as shown in the scatter plot in the
left panel of Figure 5. In the centre panel are the quartile probability
contour plots of marginal density estimates $\hat{\pi}_{1}\hat{f}_{1}$ (pink
solid lines), $\hat{\pi}_{2}\hat{f}_{2}$ (green dotted lines), and
$\hat{\pi}_{3}\hat{f}_{3}$ (blue dashed lines), where $\hat{f}_{1}$ is the
density estimate for 07:00–08:00, $\hat{f}_{2}$ for 13:00–14:00, and
$\hat{f}_{3}$ for 19:00–20:00. As the marginal density contours have
considerable overlap in the central regions, it is difficult to decide
visually which marginal density value is higher. This is resolved in the plot
of estimated class labels from the density-based classifier on the right of
Figure 5. The regions where the 07:00–08:00 measurements are more likely are
coloured in pink, the 13:00-14:00 measurements are more likely are in green,
and the 19:00–20:00 measurements are more likely are in blue. The boundaries
of these class label regions are complex and would not be well-estimated by
the linear and quadratic classifiers. The general trend is, as the day
progresses, both levels of CO2 and PM10 increase, with the increase of PM10
being more sustained.

Figure 5: Density-based classifier for air quality measurements at different
hourly intervals. (Left) Scatter plots for 07:00–08:00 ($n_{1}=1282$),
13:00–14:00 ($n_{2}=1280$), 19:00–20:00 ($n_{3}=1292$). (Centre) Quartile
probability contours of marginal density estimates. (Right) Class label
estimates.

The command to compute a tidy kernel classifier is tidy_kda. It requires a
grouped tibble as its input (air_gr), grouped by the class factor variable
(time). To produce the marginal densities plot for the density-based
classifier in the centre panel in Figure 5:

    R> ## tidy density-based classifier contours
    R> air_gr <- filter(air, time %in% c("07:00", "13:00", "19:00"))
    R> air_gr <- mutate(air_gr, time=droplevels(time))
    R> t4 <- tidy_kda(air_gr)
    R> ggplot(t4, aes(x=co2, y=pm10)) + geom_contour_ks(aes(colour=time))

The layer function geom_contour_ks draws the contour lines for tidy_ks
objects. In addition to the columns already present in the density estimate,
the extra columns in the output of a density-based classifier relate to the
classes: prior_prob (class sample proportion), label (estimated class label),
time (same as input class label). The structure of a density-based classifier
is similar to that for a density estimate grouped by a class variable.

Our geospatial data sample comprises the combined G. hakeoides (pink circles,
$n_{1}=207$) and G. paradoxa (green triangles, $n_{2}=358$) locations, as
shown in the scatter plot in the left panel of Figure 6. In the centre panel
are the quartile probability contour plots of marginal density estimates
$\hat{\pi}_{1}\hat{f}_{1}$ (pink solid lines) and $\hat{\pi}_{2}\hat{f}_{2}$
(green dotted lines), where $\hat{f}_{1}$ is the density estimate for G.
hakeoides, and $\hat{f}_{2}$ for G. paradoxa. The regions where G. hakeoides
is more likely are coloured in pink, and where G. paradoxa is more likely are
in green. For display purposes, the class labels have been truncated to the
convex hull of the marginal density estimates so that they remain over the
land area (in grey). We observe again that boundaries of these class label
regions are complex and would not be well-estimated by the linear and
quadratic classifiers.

Figure 6: Density-based classifier for Grevillea locations. (Left) Scatter
plots for G. hakeoides ($n_{1}=207$), G. paradoxa ($n_{2}=358$). (Centre)
Quartile probability contours of marginal density estimates. (Right) Class
label estimates.

A geospatial density-based classifier requires the input (grevilleasf_gr) to
be grouped by the class factor variable (species):

    R> ## geospatial density-based classifier
    R> grevilleasf_gr <- filter(grevilleasf, species %in% c("paradoxa","hakeoides"))
    R> grevilleasf_gr <- mutate(grevilleasf_gr, species=factor(species))
    R> grevilleasf_gr <- group_by(grevilleasf_gr, species)
    R> s4 <- st_kda(grevilleasf_gr)

The estimated class labels are stored in the sf_ks object in the grid field as
a collection of rectangular polygons. To plot these class labels, as in the
right panel in Figure 6, we call ggplot2::geom_sf on the grid field for a
ggplot2 plot, and plot(, which_geometry=="grid") for a base R plot.

    R> ## geospatial density-based classifier geom_sf plot
    R> ggplot(s4) + geom_sf(data=s4$grid, aes(fill=label), alpha=0.2, colour=NA)
    R> ## base R plot
    R> plot(s4, which_geometry="grid", border=NA)

The question of optimal bandwidths for a density-based classifier is more
complicated than that for a density estimate. We opt for a simple and robust
implementation in the eks package, where tidy_kda and st_kda call ks::Hpi for
each class data sub-sample. These class-wise optimal bandwidths are known to
asymptotically minimise the misclassification error, i.e., the probability
that we do not classify a candidate point in class $j$ given that it is drawn
from class $j$, $\operatorname{\bm{\mathbb{P}}}\\{\hat{c}({\bf X})\neq j|{\bf
X}\sim f_{j}\\}$. Whilst there is an intuitive appeal in selecting bandwidths
to exactly minimise the misclassification error, it is not clear how much is
gained in practise with this more complicated approach over the simpler
bandwidths. Moreover, there are currently no efficient algorithms to compute
these more complicated bandwidths. See Chacón and Duong, (2018, Section 7.2)
for a discussion.

##  4 Density derivative estimation

Crucial information about the structure of a data set is not always revealed
by examining solely the density values, and can only be discerned via the
density derivatives. For example, the local minima/maxima of the data density
are characterised as the locations where the first derivative is the
identically zero. A recent example of the utility of density derivative
estimates in data analysis is the segmentation of digital images, which
utilised the first density derivative of pixel colour-locations to guide the
search for similar image segments more efficiently than using only the density
of the pixel colour-locations (Beck et al.,, 2016).

With the same data as for the density estimation case, i.e., ${\bf
X}_{1},\dots,{\bf X}_{n}$ is a random sample drawn from the common density
function $f$, our goal is to estimate the first (gradient) derivative of the
unknown density $f$. For 2-dimensional data, the gradient of a density
function $f$ is comprised of two partial derivatives $\mathsf{D}f=[\partial
f/\partial x_{1},\partial f/\partial x_{2}]$. The kernel estimate of the
density gradient is given by

$\mathsf{D}\hat{f}_{\bf H}({\bm{x}})=n^{-1}\sum_{i=1}^{n}\mathsf{D}K_{\bf
H}({\bm{x}}-{\bf X}_{i})$ (3)

where the gradient kernel function is $\mathsf{D}K_{\bf
H}({\bm{x}})=-(2\pi)^{-1}|{\bf H}|^{-1/2}{\bf
H}^{-1}{\bm{x}}\exp(-\tfrac{1}{2}{\bm{x}}^{\top}{\bf H}^{-1}{\bm{x}})$. Since
there are two components of the density gradient, it can be visualised using
two separate plots, one for each partial derivative. A more concise
alternative is a quiver plot, in which arrows, whose length and direction are
determined by the gradient, are drawn at each point in the estimation grid.
The right panel of Figure 7 is the quiver plot for the density gradient
estimate for the air quality measurements from 13:00 to 14:00, superposed on
the density estimate. The arrows for the density gradient point towards to the
peaks of the modal regions. These arrows are longer where the density gradient
is steeper, and they are shorter in the density tails where the slope is
flatter. These density gradients indicate the rate of change in the data
density, which is not easy to ascertain from the density levels themselves in
the underlying density contour plot.

---
Figure 7: Quiver plots of density gradient estimate, superposed over density
estimates. (Left) Air quality measurements for 13:00–14:00 ($n=1280$). (Right)
G. yorkrakinensis locations ($n=93$).

The command for a tidy density gradient estimate is tidy_kdde(,
deriv_order=1). The function tidy_kquiver converts the output from tidy_kdde
into a format suitable for the quiver plot layer function
ggquiver::geom_quiver (O’Hara-Wild,, 2019). The code to produce a quiver plot
superposed on a density estimate is

    R> ## tidy density gradient estimate
    R> air3 <- ungroup(filter(air, time=="13:00"))
    R> air3 <- select(air3, co2, pm10)
    R> t5 <- tidy_kdde(air3, deriv_order=1)
    R> t6 <- tidy_kquiver(t6, thin=9)
    R> ggplot(t1, aes(x=co2,y=pm10)) + ggquiver::geom_quiver(data=t6, aes(u=u,v=v))

The output from tidy_kdde is a tibble which is grouped by deriv_group. The
columns present in a density estimate are also present in a density derivative
estimate, along with some additional columns relating to the derivative:
deriv_order (derivative order, 1 for the gradient), deriv_ind (partial
derivative enumeration, from 1 to 2), deriv_group (partial derivative indices
(1,0), (0,1) which correspond to $\partial/\partial x_{1},\partial/\partial
x_{2}$ respectively).

For geospatial data, the left panel in Figure 7 is the quiver plot for the
density gradient estimate for G. yorkrakinensis, superposed on the density
estimate. Whilst with st_kquiver we can compute a geospatial output,
ggplot2::geom_sf is not able plot arrows, and it is not possible to overlay a
ggquiver::geom_quiver layer over a geom_sf layer. The current work-around is
to overlay a ggplot2::geom_segment layer over a geom_sf layer, with some trial
and error required in grid::arrow to produce suitable arrows.

    R> ## geospatial density gradient estimate geom_sf plot
    R> s5 <- st_kdde(yorkr, deriv_order=1)
    R> s6 <- st_kquiver(s5, thin=9)
    R> ggplot(s1) + geom_segment(data=s6$sf, aes(x=lon, xend=lon_end, y=lat,
    +  yend=lat_end), arrow=grid::arrow(length=0.1*s6$sf$len))

On the other hand, for a base R plot, the display of geospatial and tidy data
are freely interchangeable, so we can overlay the quiver plot plot(,
display="quiver") for a kernel density gradient estimate from the ks package.

    R> ## geospatial density gradient estimate base R plot
    R> plot(s6$tidy_ks$ks[[1]], display="quiver")

For optimal bandwidth selection for kernel density gradient estimates, it is
crucial to note that the optimal bandwidth matrix for $\mathsf{D}\hat{f}_{\bf
H}$ is not the same as that for $\hat{f}_{\bf H}$. For a density estimate the
optimality criterion is $M({\bf
H})=\int\operatorname{\bm{\mathbb{E}}}[\hat{f}_{\bf
H}({\bm{x}})-f({\bm{x}})]^{2}\,\mathsf{d}{\bm{x}}$, whereas the criterion for
a density gradient estimate is $M_{1}({\bf
H})=\int\operatorname{\bm{\mathbb{E}}}\lVert\mathsf{D}\hat{f}_{\bf
H}({\bm{x}})-\mathsf{D}f({\bm{x}})\lVert^{2}\,\mathsf{d}{\bm{x}}$. Since
$M\neq M_{1}$ then their minimisers are also not equal in general. The default
optimal bandwidth for the density gradient estimate in the eks package is the
plug-in bandwidth (Chacón and Duong,, 2010) obtained from a call to ks::Hpi(,
deriv.order=1). For the air quality measurements for 13:00 to 14:00, this
bandwidth matrix is [441.0, 59.5; 59.5, 305.0]. In comparison, the optimal
bandwidth matrix for the density estimate is [495.0, 88.0; 88.0, 245.0]. For
the G. yorkrakinensis data, the optimal plug-in bandwidth matrix for the
density gradient estimate is [4.39e8, $-4.36$e8; $-4.36$e8, 7.73e8], whereas
for the density estimate, it is [8.84e8, $-8.33$e8; $-8.33$e8, 1.36e9].

##  5 Density-based clustering (unsupervised learning)

The goal of clustering is to discover homogeneous groups within a data set in
a trade-off between similarity/dissimilarity: members of the same cluster are
similar to each other while members of different clusters are dissimilar to
each other. If the $q$ unknown population clusters are
$\\{C_{1},\dots,C_{q}\\}$, then the cluster labelling function is
$c({\bm{x}})=j$ whenever a candidate point ${\bm{x}}$ belongs to cluster
$C_{j}$. Whilst we are able to estimate the cluster labelling function for all
candidate points, for the vast majority of data analysis cases, it is
sufficient to compute $\hat{c}({\bf X}_{1}),\dots,\hat{c}({\bf X}_{n})$ for
the data sample ${\bf X}_{1},\dots,{\bf X}_{n}$. Since the cluster labels are
unknown, then this is an unsupervised learning problem.

Many clustering algorithms have been proposed in the literature. Our chosen
approach is density-based clustering, where a cluster is a data-rich region
(high density values) which is separated from another data-rich region by a
data-poor region (low density values). Thus we associate each data point to
its ‘most representative’ data-rich region. In the eks package, this
association is carried out with a mean shift algorithm (Fukunaga and
Hostetler,, 1975). For a data point ${\bf X}_{i}$, we initialise a sequence
with ${\bf X}_{i,0}={\bf X}_{i}$, then we iterate the recurrence equation

${\bf X}_{i,k+1}={\bf X}_{i,k}+{\bf H}^{-1}\mathsf{D}\hat{f}_{\bf H}({\bf
X}_{k})\big{/}\hat{f}_{\bf H}({\bf X}_{k}),$

where $\hat{f}_{\bf H}$ is a density estimate and $\mathsf{D}\hat{f}_{\bf H}$
is a density gradient estimate. This recurrence equation is closely related to
the well-known gradient ascent algorithm, with the improvement that
accelerates the convergence of the recurrence iterations in regions of low
data density. A more computationally stable form of the mean shift recurrence
equation, since it avoids the explicit computation of the density estimate
$\hat{f}_{\bf H}$ and density derivative estimate $\mathsf{D}\hat{f}_{\bf H}$,
is

${\bf X}_{i,k+1}={\bf X}_{i,k}+{\bm{\beta}}_{\bf H}({\bf
X}_{i,k})=\frac{\sum_{\ell=1}^{n}{\bf X}_{\ell}\,g\big{(}({\bf X}_{i,k}-{\bf
X}_{\ell})^{\top}{\bf H}^{-1}({\bf X}_{i,k}-{\bf
X}_{\ell})\big{)}}{\sum_{\ell=1}^{n}g\big{(}({\bf X}_{i,k}-{\bf
X}_{\ell})^{\top}{\bf H}^{-1}({\bf X}_{i,k}-{\bf X}_{\ell})\big{)}}$ (4)

where $g(x)=x\exp(-\tfrac{1}{2}x)$ and ${\bm{\beta}}_{\bf
H}({\bm{x}})=\frac{\textstyle\sum_{\ell=1}^{n}{\bf X}_{\ell}g(({\bm{x}}-{\bf
X}_{\ell})^{\top}{\bf H}^{-1}({\bm{x}}-{\bf
X}_{\ell}))}{\textstyle\sum_{\ell=1}^{n}g(({\bm{x}}-{\bf X}_{\ell})^{\top}{\bf
H}^{-1}({\bm{x}}-{\bf X}_{\ell}))}-{\bm{x}}$. This ${\bm{\beta}}_{\bf H}$ is
known as the mean shift, since it is the difference between the current
iterate and a weighted mean of all data points. For our stopping rule, we
iterate the recurrence in Equation (4) until either we reach a maximum number
of iterations (400) or that the distance between subsequent iterations is less
than 0.001 times the minimal marginal IQR (interquartile range) of the input
data. This heuristic stopping rule gives sensible results in most cases.

The result is a sequence of points $\\{{\bf X}_{i,0},{\bf X}_{i,1},\dots\\}$
which traces out a path, along the steepest ascent of the density gradient,
from the data point ${\bf X}_{i}$ to the mode of the associated data-rich
region. The data-rich regions are the ‘basins of attraction’ of the density
gradient ascent. If the data points are associated with the same mode, then
they are considered to be members of the same cluster. Thus the number of
clusters is equal to the number of these basins of attraction. For more
details on mean shift and other forms of density-based clustering, see Chacón
and Duong, (2018, Section 6.2).

The result of the mean shift clustering on the $n=1280$ air quality
measurements from 13:00 to 14:00 into 5 clusters is displayed on the left
panel in Figure 8. Observe that we do not need to specify the number of
clusters in advance, and the clusters can be of any arbitrary shape. These
represent two important advantages over $k$-means clustering, which requires
an a priori number of clusters, and whose cluster shapes are restricted to be
(intersections of) ellipses. Cluster #4 (blue crosses) and #5 (magenta boxed
crosses) are the most separate from the other clusters. The points at the
edges of cluster #3 (green squares), cluster #1 (red circles) and cluster #2
(khaki triangles) are close to together, and $k$-means clustering tends to
assign them to the same cluster, whereas the directionality of the mean shift
assigns them to different clusters. Since the mean shift relies on the density
gradient ascent paths, we overlay the arrows of the quiver plot of the density
gradient on the convex hulls of the mean shift clusters on the right of Figure
8. We observe that the gradient ascent arrows within each cluster are oriented
towards the modes.

Figure 8: Mean shift clusters for the air quality measurements 13:00–14:00
($n=1280$). (Left) Cluster members. (Right) Cluster convex hulls, superposed
over the quiver plot of its density gradient estimate.

The command for mean shift clustering for tidy data is tidy_kms. The output is
similar to that for a single density estimate, except that the data points are
returned rather than the estimation grid points, and that estimate indicates
the estimated cluster label rather than the density estimate value.

    R> ## tidy mean shift clusters
    R> t7 <- tidy_kms(air3)
    R> ggplot(t7, aes(x=co2, y=pm10)) + geom_point(aes(colour=estimate))

Since the direction along which the data points are shifted is directly
related to the density gradient, the default bandwidth for mean shift
clustering in tidy_kms is the plug-in bandwidth computed by ks::Hpi(,
deriv.order=1). For the air quality measurements for 13:00 to 14:00, this
bandwidth matrix is [441.0, 59.5; 59.5, 305.0]. The bandwidth choice is made
with the goal of optimal identification of the density gradient ascent paths.
It is also supported by the results that the optimal bandwidth for estimating
the mode of a density is closely related to the optimal bandwidth for density
gradient estimation (Chacón and Duong,, 2018, p. 138).

The result of the mean shift clustering on the $n=93$ G. yorkrakinensis
locations into 4 clusters is displayed on the left panel in Figure 9. Cluster
#4 (magenta crosses) is the most northerly and most separate from the other
clusters. Cluster #2 (green triangles) forms the most southerly cluster and is
also well-separated. The points on the right edge of cluster #1 (red circles)
are close to those on the left edge of cluster #3 (cyan squares), though the
directionality of the mean shift, as indicated by the black arrows of the
density gradient, assigns them to different clusters.

Figure 9: Mean shift clusters for G. yorkrakinensis ($n=93$). (Left) Cluster
members. (Right) Cluster convex hulls, superposed over the quiver plot of its
density gradient estimate.

The command for mean shift clustering for geospatial data is st_kms. To
produce the mean shift clusters, with at least 3 members in each cluster, the
code snippet is:

    R> ## geospatial mean shift clusters geom_sf plot
    R> s7 <- st_kms(yorkr, min.clust.size=3)
    R> ggplot(s7) + geom_sf(aes(colour=estimate))
    R> ## base R plot
    R> plot(s7, pch=16)

The default bandwidth for mean shift clustering in st_kms is the plug-in
bandwidth computed by ks::Hpi(, deriv.order=1). For the G. yorkrakinensis
locations, this is [4.39e8, $-$4.36e8; $-$4.36e8, 7.73e8].

##  6 Other data analysis settings

There are numerous variations on the standard density estimate provided by
tidy_kde for tidy data. For bounded data, there is a boundary density estimate
(tidy_kde_boundary) and a truncated density estimate (tidy_kde_truncate). To
allow for the kernel function and/or the bandwidth matrix to vary, there are
the sample point density estimate (tidy_kde_sp) and balloon density estimate
(tidy_kde_balloon). For data observed with errors, there is the deconvolved
density estimate (tidy_kdcde). There are distribution-based estimates, with a
cumulative distribution estimate (tidy_kcde), a copula estimate
(tidy_kcopula), and a ROC (receiver operating characteristic) curve
(tidy_kroc). Other functionality includes second (curvature) density
derivative estimation (tidy_kdde(, deriv_order=2)), density ridge estimation
(tidy_kdr), significance testing for density differences (tidy_kde_local_test)
and significance testing for modal regions (tidy_kfs). These are all
implemented for geospatial data as st_k*. All of these utilise the appropriate
default bandwidth selector from the ks package. For brevity, we do not
illustrate them here: their usage is demonstrated in their help pages
contained in the eks package, and the details of the statistical framework in
these data analysis settings are provided in Chacón and Duong, (2018). Tidy
kernel smoothers also are applicable to tidy univariate data.

The ability to export the geospatial kernel estimates to standard vectorial
geospatial data formats extends the functionality of the eks package to GIS
software. The commands to export to the geopackage format are:

    R> ## export to vectorial geospatial format
    R> sf::write_sf(yorkr, dsn="grevillea.gpkg", layer="yorkr")
    R> sf::write_sf(st_get_contour(s1), dsn="grevillea.gpkg", layer="yorkr_cont")
    R> sf::write_sf(s1$grid, dsn="grevillea.gpkg", layer="yorkr_grid")
    R> sf::write_sf(s6$sf, dsn="grevillea.gpkg", layer="yorkr_quiver")

The grevillea.gpkg geopackage consists of four layers: yorkr for the point
geometries of the G. yorkrakinensis locations, yorkr_cont for the multi-
polygons of the quartile contour regions of the density estimate, yorkr_grid
for the rectangular polygons of the estimation grid, and yorkr_quiver for the
linestrings of the density gradient flows.

This grevillea.gpkg geopackage can be subsequently employed in QGIS
(QGIS.org,, 2021), which is an industry standard software for GIS
practitioners since it offers features that are not available in R. For
example, it has an interactive point-and-click interface, and it incorporates
fast rendering of the OpenStreetMap base maps. A screenshot from a QGIS
analysis for a quiver plot overlaid on a density estimate is given in the left
panel of Figure 10. Recall that quiver plots can be difficult to produce with
geospatial data in ggplot2 graphics, since the arrows require trial and error
to display suitably with ggplot2::geom_segment. In contrast, quiver plots are
straightforward in QGIS since rescaleable arrows are a native feature.

Figure 10: Screenshots of QGIS analysis for G. yorkrakinensis ($n=93$). (Left)
Contour plot of vectorial density estimate and quiver plot of vectorial
density gradient estimate. (Right) Heat map of raster density estimate.

In addition, QGIS efficiently handles raster geospatial data. Whilst the grid
field of a kernel estimate consists of the rectangular polygons for each pixel
of the estimation grid, it can be converted to a raster via the stars package.
The heat map of the converted raster is displayed in QGIS on the right of
Figure 10.

    R> ## export to raster format
    R> stars::write_stars(stars::st_rasterize(s1$grid), dsn="yorkr.tif")

##  7 Conclusion

We have introduced a new R package eks which serves as a bridge from the
comprehensive suite of kernel smoothers in the ks package to the tidyverse and
geospatial analysis. A wide range of kernel smoothing methods are available,
which (i) improve on the existing kernel density estimates, and (ii) widen the
accessibility to more complex kernel-based data analyses, such as density
gradient estimation, density-based classification (supervised learning) and
mean shift clustering (unsupervised learning). The eks package provides
practitioners with additional tools to create compelling statistical
visualisations from kernel smoothers, whether they are using tidy or
geospatial data, or whether they are using base R or tidyverse graphics.

## References

  * Baíllo and Chacón, (2021) Baíllo, A. and Chacón, J. E. (2021).  Chapter 1 - statistical outline of animal home ranges: An application of set estimation.  In Srinivasa Rao, A. S. and Rao, C., editors, Data Science: Theory and Applications, volume 44 of Handbook of Statistics, pages 3–37. Elsevier.
  * Beck et al., (2016) Beck, G., Duong, T., Azzag, H., and Lebbah, M. (2016).  Distributed mean shift clustering with approximate nearest neighbours.  In Proceedings of the 2016 International Conference on Neural Networks (IJCNN), pages 3110–3115.
  * Béranger et al., (2019) Béranger, B., Duong, T., Perkins-Kirkpatrick, S. E., and Sisson, S. A. (2019).  Tail density estimation for exploratory data analysis using kernel methods.  Journal of Nonparametric Statistics, 31:144–174.
  * Bowman and Foster, (1993) Bowman, A. W. and Foster, P. (1993).  Density based exploration of bivariate data.  Statistics and Computing, 3:171–177.
  * Chacón and Duong, (2010) Chacón, J. E. and Duong, T. (2010).  Multivariate plug-in bandwidth selection with unconstrained bandwidth matrices.  Test, 19:375–398.
  * Chacón and Duong, (2018) Chacón, J. E. and Duong, T. (2018).  Multivariate Kernel Smoothing and Its Applications.  Chapman and Hall/CRC, Boca Raton.
  * Chacón et al., (2011) Chacón, J. E., Duong, T., and Wand, M. P. (2011).  Asymptotics for general multivariate kernel density derivative estimators.  Statistica Sinica, 21:807–840.
  * Duong, (2007) Duong, T. (2007).  ks: Kernel density estimation and kernel discriminant analysis for multivariate data in R.  Journal of Statistical Software, 21(7):1–16.
  * Duong, (2023) Duong, T. (2023).  eks: Tidy and Geospatial Kernel Smoothing.  R package version 1.0.3.
  * Duong and Hazelton, (2003) Duong, T. and Hazelton, M. L. (2003).  Plug-in bandwidth matrices for bivariate kernel density estimation.  Journal of Nonparametric Statistics, 15:17–30.
  * Fukunaga and Hostetler, (1975) Fukunaga, K. and Hostetler, L. (1975).  The estimation of the gradient of a density function, with applications in pattern recognition.  IEEE Transactions on Information Theory, 21:32–40.
  * Hyndman, (1996) Hyndman, R. J. (1996).  Computing and graphing highest density regions.  American Statistician, 50:120–126.
  * Kalair and Connaughton, (2021) Kalair, K. and Connaughton, C. (2021).  Anomaly detection and classification in traffic flow data from fluctuations in the flow–density relationship.  Transportation Research Part C: Emerging Technologies, 127:103178.
  * Myers et al., (2000) Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A. B., and Kent, J. (2000).  Biodiversity hotspots for conservation priorities.  Nature, 403:853–858.
  * OGC, (2021) OGC (2021).  Simple feature access.  Open Geospatial Consortium. https://www.ogc.org/standards/sfa.
  * O’Hara-Wild, (2019) O’Hara-Wild, M. (2019).  ggquiver: Quiver Plots for ‘ggplot2’.  R package version 0.2.0.
  * Otto and Kahle, (2023) Otto, J. and Kahle, D. (2023).  ggdensity: Interpretable Bivariate Density Visualization with ’ggplot2’.  R package version 1.0.0.
  * Pebesma, (2018) Pebesma, E. (2018).  Simple features for R: Standardized support for spatial vector data.  The R Journal, 10:439–446.
  * QGIS.org, (2021) QGIS.org (2021).  QGIS Geographic Information System.  QGIS Association.
  * RATP, (2016) RATP (2016).  Qualité de l’air mesurée dans la station Châtelet.  https://data.iledefrance.fr/explore/dataset/qualite-de-lair-mesuree-dans-la-station-chatelet. Régie autonome des transports parisiens – Département Développement, Innovation et Territoires.  Accessed 2017-09-27.
  * Rudis et al., (2017) Rudis, B., Bolker, B., and Schulz, J. (2017).  ggalt: Extra Coordinate Systems, ‘Geoms’, Statistical Transformations, Scales and Fonts for ‘ggplot2’.  R package version 0.4.0.
  * Schauer et al., (2010) Schauer, K., Duong, T., Bleakley, K., Bardin, S., Bornens, M., and Goud, B. (2010).  Probabilistic density maps to study global endomembrane organization.  Nature Methods, 7:560–568.
  * Wickham, (2014) Wickham, H. (2014).  Tidy data.  Journal of Statistical Software, 59(10):1–23.
  * Wickham, (2016) Wickham, H. (2016).  ggplot2: Elegant Graphics for Data Analysis.  Springer-Verlag, New York.


# Characterising representation dynamics in
recurrent neural networks for object recognition

Sushrut Thorat Adrien Doerig Tim C. Kietzmann
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Institute of Cognitive Science, Osnabrück University, Germany

## Abstract

Recurrent neural networks (RNNs) have yielded promising results for both
recognizing objects in challenging conditions and modeling aspects of primate
vision. However, the representational dynamics of recurrent computations
remain poorly understood, especially in large-scale visual models. Here, we
studied such dynamics in RNNs trained for object classification on MiniEcoset,
a novel subset of ecoset. We report two main insights. First, upon inference,
representations continued to evolve after correct classification, suggesting a
lack of the notion of being “done with classification”. Second, focusing on
“readout zones” as a way to characterize the activation trajectories, we
observe that misclassified representations exhibit activation patterns with
lower L2 norm, and are positioned more peripherally in the readout zones. Such
arrangements help the misclassified representations move into the correct
zones as time progresses. Our findings generalize to networks with lateral and
top-down connections, and include both additive and multiplicative
interactions with the bottom-up sweep. The results therefore contribute to a
general understanding of RNN dynamics in naturalistic tasks. We hope that the
analysis framework will aid future investigations of other types of RNNs,
including understanding of representational dynamics in primate vision111This
article is a revision of the 2023 Conference on Cognitive Computational
Neuroscience (CCN) paper, in which we present a new analysis in the Appendix,
and include suggestions made by the CCN reviewers..

> Keywords: recurrent neural networks, object recognition, neural
> representations, dynamics, naturalistic tasks, readout zones

##  1 Introduction

Feedback connections are ubiquitous in brains Felleman Van Essen (1991). The
resulting recurrent computations are advantageous in challenging conditions
such as recognizing objects in clutter Wyatte . (2014); Kreiman Serre (2020)
and natural scenes Spoerer . (2020). Research into the representation dynamics
underlying recurrent computations is nascent but accelerating Mante . (2013);
Zamir . (2017); Quax van Gerven (2018); Mastrogiuseppe Ostojic (2018); van
Bergen Kriegeskorte (2020); Thorat . (2021); Lindsay . (2022); Driscoll .
(2022). Moving to a more naturalistic setting, this work investigates
representations and their dynamics in a deep recurrent convolutional neural
network (RNN), as they contribute to improving classification responses to
natural images. While we provide novel insights into temporal trajectories of
the RNNs, the developed framework applies more broadly and can be applied to
both, artificial and biological neural network dynamics, and hence contributes
to the toolbox available to researchers interested in modelling vision with
deep neural networks Doerig . (2023).

Figure 1: (A) Architecture of the recurrent neural network. The lateral or
top-down connections interact with bottom-up processing additively or
multiplicatively. (B) Readout vectors capture referential semantic features of
the data.

##  2 Model system and dataset

In our RNN models222The training and evaluation scripts can be found at:
github.com/KietzmannLab/BLT-pytorch-CCN23, lateral or local top-down
connections are included (Fig. 1A). Such RNNs have been used as models of
human neural dynamics and behavior Kietzmann . (2019); Spoerer . (2020);
Doerig . (2022). The lateral and top-down connections interacted with the
bottom-up sweep through either additive or multiplicative interactions. The
RNNs were unrolled for $10$ timesteps. The RNNs were trained to classify the
input images at each timestep (their readouts had no bias terms; see Appendix
6.1.1). The $64\times 64\,$px RGB images were taken from
MiniEcoset333MiniEcoset can be found at: osf.io/msna2/, which is a novel
subset of ecoset Mehrer . (2021) containing $100$ object classes that follow a
hierarchical object structure.

##  3 Analysis

We start our analyses by focusing on an RNN with lateral connections which
interact with the feedforward sweep additively. Please note that these results
generalize across RNN configuration (Fig. 3A).

###  3.1 Learned categorical structure

We start our analysis by asking whether the RNN successfully learns the
hierarchical structure encoded in the dataset statistics. To do so, we
computed the similarities between the readout vectors (rows of the readout
weight matrix, corresponding to connections from the final AvgPool layer to
each of the readout neurons), as they can give us insight into which classes
are considered similar by the RNN. Cosine similarity
($\vec{A}\cdot\vec{B}/|\vec{A}||\vec{B}|$) was computed between each pair of
the readout vectors.

Hierarchical clustering on the pairwise similarities revealed meaningful
clusters (Fig. 1B) resembling the dataset structure and the animacy
organization observed in primate brains Grill-Spector Weiner (2014). This
suggests that our choice of architecture and dataset leads to an interpretable
feature extractor.

###  3.2 Convergent representation dynamics

Next, we moved into analysing the representational dynamics of the RNNs,
asking whether they exhibit a signature of being “done with classification”,
as expected in a stable RNN with attractor dynamics Linsley . (2020).
Additionally, we asked if the changes in pre-readout representations (i.e.,
final AvgPool layer activations) are smaller for images that are already
correctly classified as opposed to images that are not yet correctly
classified. For this analysis, we focused on images that were classified
correctly and consistently starting from a given timestep $t$ (termed stable
classification with $t_{stable}=t$; we only consider these images for
subsequent analyses). To define representational changes, we analysed the
$l^{2}$-norms of the change in representations across time, as a function of
$t_{stable}$.

As seen in Fig. 2A (left), the amount of representational change did not
depend on $t_{stable}$: the changes in representations were not smaller for
images that were classified correctly at earlier timesteps. However, the
change in all representations did decrease with timesteps. These results
indicate that although all representations “settle” across time, the rate of
settling is independent of the correctness of classification. Interestingly,
as seen in Fig. 2A (right), this reduction in the rate of change was also
observed pre-training, suggesting these dynamics are a property of the network
architecture. Finally, note that in contrast to previous findings Linsley .
(2020), these RNNs exhibits stable state dynamics despite them being trained
with backpropagation through time (BPTT), as discussed in Appendix 6.3.

Figure 2: (A) The amount of change in representations does not depend on the
correctness of classification, in both trained and random RNNs. (B) Signatures
of stable correct classification: the norm of the representation and its
cosine similarity to the readout vector of the current class are higher. (C)
Signature of the future correct class: for currently misclassified
representations, the cosine similarity to the correct class readout vector is
higher.

###  3.3 Signatures of stable classification

Originating from the observation that, on average, representations move the
same distance regardless of correct classification, we hypothesized that
representations that are able to transition into another class may initially
be closer to the decision boundary, whereas the ones that do not transition
are initially far from the boundary (and are therefore unable to leave the
current class). As we show in the Appendix 6.1, in networks with linear
readouts and argmax decisions, the “readout zones”, in which representations
are assigned to a given class, resemble conical structures (a 2D schematic is
shown in Fig. 3B). Given this structure, being closer to the decision boundary
either entails having a lower L2 norm or having a lower cosine similarity with
the readout vector (see Appendix 6.1.1 for further explanation). To explore
this hypothesis, we assessed whether currently incorrectly-classified
representations (that will eventually become correct) indeed have lower norms
and/or lower cosine similarities with the readout vector of the current class.

At each timestep $t$, we compared both properties of the representations with
$t_{stable}\leq t$ (i.e., currently correct) and the representations with
$t_{stable}>t$ (i.e., currently incorrect): their norms, and their cosine
similarities to the current readout. As seen in Fig. 2B, both properties were
smaller for $t_{stable}<t$ than for $t_{stable}\leq t$: the norms and cosine
similarities were lower for representations that were incorrectly classified
at a given timestep. As seen in Fig. 3A, these patterns (averaged across
timesteps) are independent of the kind of feedback used or how it interacts
with the bottom-up sweep. This confirms the hypothesis that currently
incorrect representations are closer to the decision boundary.

What constrains incorrect images to be closer to the decision boundary? There
are two main possibilities: either any feedforward sweep, including in a
purely feedforward network, automatically projects them to this position, or
the feedforward sweep is shaped by the fact that recurrent computations move
representations the same distance regardless of correct classification. To
answer this question, we tested if the norms of two feedforward networks could
predict the how fast images are correctly classified by the RNN. They do, as
can be seen in Appendix 6.2. This suggests that the requirement of the
recurrent computations (representations moving out of a class should be closer
to the decision boundary) are satisfied by the representations instantiated by
the feedforward sweep. The reason why incorrect images are closer to the
decision boundary is independent of recurrence, and the effect of recurrent
computations is to move these representations from the incorrect to the
correct class. What properties of the images lead to their representations
being initializated closer to the decision boundary remains to be explored.

###  3.4 Signatures of the correct class

We have now established that currently misclassified objects reside closer to
the decision boundary (in the incorrect readout zone). Do these currently
incorrectly-classified representations exhibit any signatures of their correct
classes? Evidence for this would be provided if the cosine similarity of an
incorrectly-classified representation to its correct class readout vector was
higher than its cosine similarity to the readout vectors corresponding to the
correct class of other incorrectly-classified representations in the same
readout zone (see Fig. 3B for a schematic). As seen in Fig. 2C, the cosine
similarity of the incorrectly-classified representations to the corresponding
correct class readout vector is indeed higher than the cosine similarity to
other correct classes’ readout vectors. Hence, there are signatures of the
correct classes in the incorrectly classified representations. This pattern
(averaged across timesteps) is independent of the kind of feedback (lateral
vs. top-down) and how it interacts with the bottom-up sweep (additive vs.
multiplicative; Fig. 3A).

An intriguing question that arises from this is whether and how recurrent
computations utilize these nascent features to correct the classification.
Future work in understanding these dynamics shall consider: How do the
incorrectly classified representations move through other classes to arrive at
their correct classes? How do the feedback connections hierarchically (given
Fig. 1B) constrain the feedforward sweep to lead to those trajectories? Are
similar dynamics/representations found in biological visual systems?

Figure 3: (A) The results, shown in the previous figures, are independent of
the type of feedback (lateral vs top-down) and its interaction with the
bottom-up sweep (additive or multiplicative). (B) Schematic of the learned
representation arrangements: incorrectly classified representations have a
lower norm and lie closer to the class boundaries which allows for easier
transition to another class.

##  4 Conclusions

In the RNNs studied here, the magnitude of changes in network activations are
surprisingly similar across images and decrease with model timesteps. This
shows that the extent of recurrent dynamics experienced by image
representations does not depend on the correctness of classification. In
addition, we highlight an interesting representation arrangement, presented
schematically in Fig. 3B: image representations that are currently incorrectly
classified (red and blue squares) have lower norms, and are closer to the
current readout zone’s decision boundary. The initial norm of the
representation depends on the alignment of the image features with the
feedforward weights, and can be seen as indicating the certainty of the
network’s inference after the feedforward sweep. For representations where
certainty is low, recurrence can more easily move them towards the correct
readout zone.

This work reported our first advances in deriving a framework for
understanding representational dynamics in RNNs trained on naturalistic
images, which we hope will further clarify how recurrent systems, both
artificial and biological, reach their decisions. Future work should
investigate the representation trajectories in other recurrent systems,
including spatiotemporal data from the primate visual system.

##  5 Acknowledgments

The project was partially funded by the European Union (ERC, TIME, Project
$101039524$). Compute resources were funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation, Project number
$456666331$).

## References

  * Doerig . (2022) doerig2022semanticDoerig, A., Kietzmann, TC., Allen, E., Wu, Y., Naselaris, T., Kay, K. Charest, I.  2022\.  Semantic scene descriptions as an objective of human vision Semantic scene descriptions as an objective of human vision.  arXiv preprint arXiv:2209.11737.
  * Doerig . (2023) doerig2023neuroconnectionistDoerig, A., Sommers, RP., Seeliger, K., Richards, B., Ismael, J., Lindsay, GW.others  2023\.  The neuroconnectionist research programme The neuroconnectionist research programme.  Nature Reviews Neuroscience1–20.
  * Driscoll . (2022) driscoll2022flexibleDriscoll, L., Shenoy, K. Sussillo, D.  2022\.  Flexible multitask computation in recurrent networks utilizes shared dynamical motifs Flexible multitask computation in recurrent networks utilizes shared dynamical motifs.  bioRxiv2022–08.
  * Felleman Van Essen (1991) felleman1991distributedFelleman, DJ. Van Essen, DC.  1991\.  Distributed hierarchical processing in the primate cerebral cortex. Distributed hierarchical processing in the primate cerebral cortex.  Cerebral cortex (New York, NY: 1991)111–47.
  * Grill-Spector Weiner (2014) grill2014functionalGrill-Spector, K. Weiner, KS.  2014\.  The functional architecture of the ventral temporal cortex and its role in categorization The functional architecture of the ventral temporal cortex and its role in categorization.  Nature Reviews Neuroscience158536–548.
  * He . (2016) he2016deepHe, K., Zhang, X., Ren, S. Sun, J.  2016\.  Deep residual learning for image recognition Deep residual learning for image recognition.  Proceedings of the IEEE conference on computer vision and pattern recognition Proceedings of the ieee conference on computer vision and pattern recognition ( 770–778).
  * Kietzmann . (2019) kietzmann2019recurrenceKietzmann, TC., Spoerer, CJ., Sörensen, LK., Cichy, RM., Hauk, O. Kriegeskorte, N.  2019\.  Recurrence is required to capture the representational dynamics of the human visual system Recurrence is required to capture the representational dynamics of the human visual system.  Proceedings of the National Academy of Sciences1164321854–21863.
  * Kreiman Serre (2020) kreiman2020beyondKreiman, G. Serre, T.  2020\.  Beyond the feedforward sweep: feedback computations in the visual cortex Beyond the feedforward sweep: feedback computations in the visual cortex.  Annals of the New York Academy of Sciences14641222–241.
  * Lindsay . (2022) lindsay2022bioLindsay, GW., Mrsic-Flogel, TD. Sahani, M.  2022\.  Bio-inspired neural networks implement different recurrent visual processing strategies than task-trained ones do Bio-inspired neural networks implement different recurrent visual processing strategies than task-trained ones do.  bioRxiv2022–03.
  * Linsley . (2020) linsley2020stableLinsley, D., Karkada Ashok, A., Govindarajan, LN., Liu, R. Serre, T.  2020\.  Stable and expressive recurrent vision models Stable and expressive recurrent vision models.  Advances in Neural Information Processing Systems3310456–10467.
  * Mante . (2013) mante2013contextMante, V., Sussillo, D., Shenoy, KV. Newsome, WT.  2013\.  Context-dependent computation by recurrent dynamics in prefrontal cortex Context-dependent computation by recurrent dynamics in prefrontal cortex.  Nature503747478–84.
  * Mastrogiuseppe Ostojic (2018) mastrogiuseppe2018linkingMastrogiuseppe, F. Ostojic, S.  2018\.  Linking connectivity, dynamics, and computations in low-rank recurrent neural networks Linking connectivity, dynamics, and computations in low-rank recurrent neural networks.  Neuron993609–623.
  * Mehrer . (2021) mehrer2021ecologicallyMehrer, J., Spoerer, CJ., Jones, EC., Kriegeskorte, N. Kietzmann, TC.  2021\.  An ecologically motivated image dataset for deep learning yields better models of human vision An ecologically motivated image dataset for deep learning yields better models of human vision.  Proceedings of the National Academy of Sciences1188e2011417118.
  * Quax van Gerven (2018) quax2018emergentQuax, S. van Gerven, M.  2018\.  Emergent mechanisms of evidence integration in recurrent neural networks Emergent mechanisms of evidence integration in recurrent neural networks.  PLoS One1310e0205676.
  * Spoerer . (2020) spoerer2020recurrentSpoerer, CJ., Kietzmann, TC., Mehrer, J., Charest, I. Kriegeskorte, N.  2020\.  Recurrent neural networks can explain flexible trading of speed and accuracy in biological vision Recurrent neural networks can explain flexible trading of speed and accuracy in biological vision.  PLoS computational biology1610e1008215.
  * Thorat . (2021) thorat2021categoryThorat, S., Aldegheri, G. Kietzmann, TC.  2021\.  Category-orthogonal object features guide information processing in recurrent neural networks trained for object categorization Category-orthogonal object features guide information processing in recurrent neural networks trained for object categorization.  arXiv preprint arXiv:2111.07898.
  * van Bergen Kriegeskorte (2020) van2020goingvan Bergen, RS. Kriegeskorte, N.  2020\.  Going in circles is the way forward: the role of recurrence in visual inference Going in circles is the way forward: the role of recurrence in visual inference.  Current Opinion in Neurobiology65176–193.
  * Wyatte . (2014) wyatte2014earlyWyatte, D., Jilk, DJ. O’Reilly, RC.  2014\.  Early recurrent feedback facilitates visual object recognition under challenging conditions Early recurrent feedback facilitates visual object recognition under challenging conditions.  Frontiers in psychology5674.
  * Zamir . (2017) zamir2017feedbackZamir, AR., Wu, TL., Sun, L., Shen, WB., Shi, BE., Malik, J. Savarese, S.  2017\.  Feedback networks Feedback networks.  Proceedings of the IEEE conference on computer vision and pattern recognition Proceedings of the ieee conference on computer vision and pattern recognition ( 1308–1317).

##  6 Appendix

###  6.1 The shape of the readout zones

Similar to most neural networks used for object recognition, in our RNNs, the
readout layer consists of a linear transformation of the representation
followed by the softmax operation. Given representation
$r\in\mathbb{R}^{d_{r}}$, the readout weight matrix
$M\in\mathbb{R}^{d_{o}\times d_{r}}$, and bias $b\in\mathbb{R}^{d_{o}}$, the
softmaxed output responses $o\in\mathbb{R}^{d_{o}}$ are:

$\begin{split}o^{\prime}&=Mr+b\\\
o_{i}&=e^{o^{\prime}_{i}}/\sum_{j=1}^{d_{o}}e^{o^{\prime}_{j}}\end{split}$ (1)

As $e^{x}$ is a monotonic function, softmax does not change the rank-order of
the responses: $\text{argsort}(o^{\prime})=\text{argsort}(o)$. Therefore, the
output of the network is the class represented by neuron $i$, which is
associated with the readout vector $M[i,:]$ that has the largest dot product
with $r$ (if bias $b=0$, which is the case in our RNNs). The rank-order of the
dot products,
$o^{\prime}_{i}=\left\|M[i,:]\right\|\left\|r\right\|\text{cos}\theta$, is
invariant to the L2 norm of the representation, $\left\|r\right\|$, with the
other factors held constant.

Figure 4: Simulations of $4$ readout zones in a 2D representation space, with
or without readout bias (the readout weights stay constant across a run). The
corresponding readout vectors are also shown. Without readout bias, the zones
resemble triangular areas with one vertex each at origin. Readout bias can add
offsets to the zones, which are better visualized in the regime where the norm
of the representations is much higher than the norm of the readout bias. Note
that, in this low-dimensional readout regime, the readout vectors sometimes do
not fall in the corresponding readout zone. However, with higher dimensional
readouts, this does not happen (see Fig. 5)

For $d_{o}=2$, this implies the regions of representation space corresponding
to the network output decisions, “readout zones”, would resemble triangles,
with one of their vertices at origin. This is confirmed by the simulations of
$4$ readout zones, in a 2D representation space subjected to varying random
$4\times 2$ linear transformations and softmax (Fig. 4 (top)). The addition of
readout bias corresponds to adding offsets to the readout zones as seen in
Fig. 4 (middle). However, these readout zones still look like triangles in the
regime where the norm of the representations is much higher than the norm of
the readout bias (Fig. 4 (bottom)).

####  6.1.1 Implications for our results

As mentioned earlier, our RNNs have no readout bias. Hence, the readout zones
can be interpreted as higher-dimensional conical structures akin to the
triangular regions in 2D, with one of their vertices at origin, supporting the
norm-based representation arrangement in Fig. 3B. However, Fig. 4 suggests
that the readout vectors might not be inside the readout zones, which runs
counter to our cosine similarity analysis. In Fig. 4, the number of readout
vectors is higher than the dimensionality of the representational space,
whereas in our RNNs, the representational space is $512$D while the readout is
$100$D. We checked if, in our setting, the readout vectors are indeed inside
the readout zones.

First, we checked if the readout vectors are classified as the correct class,
which would indicate they lie within their corresponding readout zones. If the
readout matrix is $M\in\mathbb{R}^{100\times 512}$, the class decisions are
given by: $\text{argmax}(M\times M^{T},0)$. We found that $100\%$ of the
readout vectors (from the RNN in Fig. 2) were classified correctly i.e. they
lay in their corresponding readout zones.

Figure 5: The proportion of $100$ readout vectors laying in their
corresponding readout zones increases and quickly reaches 100% with the
dimensionality of the representational space. Trained readouts are obtained
from the readout matrix of the trained RNN (Fig. 2), by sampling the number of
representational dimensions as mentioned on the x-axis. $95\%$ confidence
intervals of the mean of $100$ sampling repetitions are shown in black.

Second, we checked if the representations are in a high-dimensional space and
the readouts span a subspace, the proportion of readout vectors laying in
their corresponding readout zones increases. This might explain the
discrepancy between the “networks” in Fig. 4 and our network (Fig. 2). We
sampled $100$ random readouts from a representational space spanning $1$ to
$50$ dimensions, and asked how many of those lay in their corresponding
readout zones, similar to the analysis above. We did the same for the readout
matrix of our trained RNN - we sampled $1$ to $50$ dimensions of the
representational space. We repeated both analysis $100$ times to get
confidence bounds. As seen in Fig. 5, all the $100$ readout vectors lay in
their corresponding readout zones once the dimensionality of the
representational space was around $30$ (much lower than the dimensionality of
the representational space of our RNNs which is $512$), for both the random
and trained readouts. The proportion of readout vectors being in their
corresponding zones increases faster for the random readouts as they probably
span a larger subspace than the trained readouts (owing to correlations across
classes), leading to a lower proportion of closely located readout vectors.

In summary, in our RNNs, indeed the readout vectors lie in readout zones which
are high-dimensional conical structures anchored to origin, as schematized in
Fig. 3B.

Do these insights generalize to RNNs with readout biases? Adding bias did not
affect the classification accuracy: both, the RNN from Fig. 2 and the same
version trained with readout bias, had a maximum accuracy of $53\%$ on the
“testplus” split of MiniEcoset. Removing the bias after training with it also
did not change the maximum accuracy. Relatedly, in the RNN trained with
readout bias, the norm of the bias was $0.27$ whereas the norms of the
representations lay in the range [$5.1$, $17.9$] with a mean of $8.1$. These
results suggest that the regime in which our RNNs with readout bias operate
would lead to readout zones similar to the ones seen in Fig. 4 (bottom), which
makes our representational arrangement insights (Fig. 3B) generalize to RNNs
with readout biases.

###  6.2 Representational norms in feedforward networks predict $t_{stable}$
in the RNN

We asked whether the arrangement of norms found in the RNN, shown in Fig. 2B,
could be due to the nature of the computations in the feedforward sweep:
images that take longer to be classified by the RNN are the ones that are not
“prototypical” of their category. This would be reflected as them having lower
alignment with the feedforward weights (feature kernels), leading to lower
neural responses i.e. norms. To test this idea, we assessed if the arrangement
of norms found in the feedforward sweep of the RNN (“BNet”; trained without
recurrence on MiniEcoset) predict the time taken by the RNN to stably classify
the images. We also assessed if this arrangement of norms is independent of
network architecture by testing the norms from a ResNet18 He . (2016) trained
on MiniEcoset.

Figure 6: The norms of the feedforward networks’ pre-readout responses to the
stably classified images, considered in Fig. 2, are related to their
$t_{stable}$ values from the RNN. A negative correlation is found for the
trained feedforward networks, in line with previous results: higher norms are
associated with quicker stable classification.

For each of the stably classified images considered in Fig. 2, we extracted
the norms of the responses at the pre-readout layers for each network, and
related them to the corresponding $t_{stable}$ values from the RNN. As seen in
Fig. 6, a negative relationship was found: representations with higher norms
in the feedforward networks corresponded to images which took longer to be
stably classified by the RNN. Crucially, this relationship was not found in a
randomly initialised BNet, indicating that training on the dataset is
essential to acquire this norm arrangement.

This result suggest that while the representational arrangement found in the
RNNs dovetails with the constraint of “equal movement across representations”
(Fig. 2A), the arrangement might not have emerged due to the constraint.
Rather the arrangement is a property of the feedforward network, which is then
used by the recurrent computations to move these non-prototypical
representations into their correct class.

###  6.3 Stability analysis

The overall decrease over time in the change in representations seen in Fig.
2A suggests that these RNNs might have convergent dynamics. However, to
establish convergent dynamics, the change in representations needs to be
assessed past the timesteps the RNNs were trained for. Each of the $4$ RNNs
(analysed in Fig. 3), which were trained with classification losses aggregated
over $10$ timesteps, were unrolled to $80$ timesteps to assess their dynamics
beyond their training. If the RNNs indeed have convergent dynamics i.e. are
stable, the change in representations should decrease with timesteps and be
much smaller than the norm of the representations. Additionally, the
classification accuracies should not drop in the timesteps beyond the training
timesteps. We monitored the relative change in the representation $R$ at
timestep $t$:
$log_{10}(\left\|R_{t}-R_{t-1}\right\|)-log_{10}(\left\|R_{t-1}\right\|)$),
and the classification accuracy on the test set of MiniEcoset, across the $80$
timesteps.

Figure 7: Assessing RNN stability. The RNNs were trained over $10$ timesteps
and were unrolled to $80$ timesteps post-training. (Left) The classification
performance of the $4$ RNNs, analysed in Fig. 3, did not drop substantially
past the point they were trained for. (Middle) The relative change in
representations decreased through time and was orders of magnitude smaller
than the norms of those representations. These results suggest these RNNs have
convergent, stable, dynamics. (Right) These stable dynamics existed pre-
training.

As seen in Fig. 7, the relative change in the representations (averaged across
the test set) decreased with timesteps and was at least $2$ orders of
magnitude smaller than the norms of those representations. This shrinking
change in representations was accompanied by a sustained classification
performance (the worst performance loss was $\sim 6\%$ of the performance at
timestep $10$). Additonally, this signature of decreasing relative change in
representations was also found in these RNNs before training.

These results suggest that these RNNs are indeed exhibiting convergent
dynamics i.e. they are stable. Note that previous research indicated that RNNs
trained with backpropagation through time (BPTT) do not show such convergent
dynamics beyond the training timesteps Linsley . (2020). Future research would
address these discrepancies in the dynamics, by assessing the points of
divergence between the RNNs used in previous research and the current class of
RNNs, in terms of the architecture, the dataset, the classification objective,
and the learning rule.