# Fitting three-dimensional Laguerre tessellations by hierarchical marked
point process models

Filip Seitl<EMAIL_ADDRESS>Jesper Møller  Viktor Beneš  Department
of Probability and Mathematical Statistics, Faculty of Mathematics and
Physics, Charles University, Sokolovská 83, 186 75, Praha 8, Czech Republic
Department of Mathematical Sciences, Aalborg University, Skjernvej 4A, DK-9220
Aalborg Ø, Denmark

###### Abstract

We present a general statistical methodology for analysing a Laguerre
tessellation data set viewed as a realization of a marked point process model.
In the first step, for the point, we use a nested sequence of multiscale
processes which constitute a flexible parametric class of pairwise interaction
point process models. In the second step, for the marks/radii conditioned on
the points, we consider various exponential family models where the canonical
sufficient statistic is based on tessellation characteristics. For each step,
parameter estimation based on maximum pseudolikelihood methods is tractable.
For model selection, we consider maximized log pseudolikelihood functions for
models of the radii conditioned on the points. Model checking is performed
using global envelopes and corresponding tests in both steps and moreover by
comparing observed and simulated tessellation characteristics in the second
step. We apply our methodology for a 3D Laguerre tessellation data set
representing the microstructure of a polycrystalline metallic material, where
simulations under a fitted model may substitute expensive laboratory
experiments.

###### keywords:

exponential family model , global envelope test , multiscale process ,
polycrystalline microstructure , pseudolikelihood

###### MSC:

[2010] 60G55

##  1 Introduction

This paper develops new statistical methodology for analysing Laguerre
tessellation data sets, where we illustrate how the methodology applies for an
example of a 3D Laguerre tessellation data set representing a polycrystalline
material.

###  1.1 Some background and motivation

In brief, a 3D Laguerre tessellation (also called a power diagram or a
generalised Voronoi tessellation) is a flexible way of modelling a subdivision
of 3D space (see Section 2 for the details). In practice the tessellation is
described by a data set $(\mathbf{x}_{n},\mathbf{r}_{n})$ where
$\mathbf{x}_{n}=\\{x_{1},\ldots,x_{n}\\}\subset\mathbb{R}^{3}$ is a finite
point pattern and $\mathbf{r}_{n}=(r_{1},\ldots,r_{n})\in(0,\infty)^{n}$ is an
associated vector of positive numbers called marks, and we identify
$(\mathbf{x}_{n},\mathbf{r}_{n})$ by the marked point pattern
$\\{(x_{1},r_{1}),\ldots,(x_{n},r_{n})\\}$ (see again Section 2 for the
details). Such data sets appear in various contexts of material science. For
example, in 3D X-ray diffraction microscopy of polycrystalline materials the
volumes and centroids of grains are obtained, and representation of such
measurements by a 3D Laguerre tessellation can be produced by various
optimization methods [Lyckegaard et al., 2011, Spettl et al., 2016, Quey and
Renversade, 2018, Kuhn et al., 2020]. It is then useful to develop stochastic
models for the tessellation since simulations can to some extent substitute
expensive laboratory experiments. However, only a few papers [including Spettl
et al., 2015, Seitl et al., 2021] have been dealing with statistical
methodology in the context of polycrystalline materials.

Inference for statistical models of Laguerre tessellation data sets is a
difficult task unless $(\mathbf{x}_{n},\mathbf{r}_{n})$ follows a simple model
such as a (marked) Poisson process which in practice is rarely a reasonable
assumption. Often the points in $\mathbf{x}_{n}$ exhibit regularity and
different marked point process models for $(\mathbf{x}_{n},\mathbf{r}_{n})$
have been suggested: simple random models for packings of hard balls [Chiu et
al., 2013] obtained by iterative procedures such as random sequential
adsorption [Lautensack, 2008, with applications to foam structures] or
variants of collective-rearrangement algorithms [Spettl et al., 2015, with
application to a polycrystalline microstructure]; and inspired by the work in
Dereudre and Lavancier [2011] on Gibbsian models for 2D Voronoi tessellations,
Seitl et al. [2021] introduced Gibbsian models for random 3D Laguerre
tessellations. These Gibbs models were first thought to be tempting to use
because they may incorporate properties of the Laguerre tessellation and
interaction between its cells. However, the models are complicated to use for
statistical inference and they are time-consuming to simulate.

###  1.2 Our contribution

In this paper we introduce hierarchical statistical models consisting of first
a parametric Gibbs point process model for $\mathbf{x}_{n}$ and second a
parametric model for $\mathbf{r}_{n}$ conditioned on $\mathbf{x}_{n}$.
Specifically, we use in the first case a nested sequence of flexible pairwise
interaction points processes called multiscale processes [Penttinen, 1984] and
in the second case various exponential family models where the canonical
sufficient statistic is based on tessellation characteristics such as surface
area or volume of cells or absolute difference in volumes of neighbouring
cells. Apart from reducing the dimension from 4 (when viewing
$(\mathbf{x}_{n},\mathbf{r}_{n})$ as a 4-dimensional point pattern) to 3 (when
considering $\mathbf{x}_{n}$), an advantage is that we specify two much
simpler cases of models with parameters which do not depend on each other.
Hence we can separate between how to simulate and estimate unknown parameters
for $\mathbf{x}_{n}$ and $\mathbf{r}_{n}\mid\mathbf{x}_{n}$, respectively. The
parameters are simply estimated by maximum pseudolikelihood methods and well-
known MCMC algorithms are used for simulations. Thereby estimation and
simulation become much faster than in Seitl et al. [2021] when fitting
specific models given in Section 4 and used for analysing the data from
Section 3 in Section 5.

A further advantage is that the model construction makes it possible to
develop a rather straightforward model selection procedure: For
$\mathbf{x}_{n}$, the procedure starts with the simplest case of a Poisson
process and continues with constructing more and more complex multiscale
processes until a satisfactory fit is obtained when considering global
envelopes and tests [Myllymäki et al., 2017] based on various functional
summary statistics. For $\mathbf{r}_{n}$ conditioned on $\mathbf{x}_{n}$, more
and more complex exponential models are developed, where we demonstrate how to
compare fitted models of the same dimension by considering maximized log
pseudolikelihood functions. Further, we evaluate selected fitted models by
comparing moment properties of tessellation characteristics under simulations
from the model with empirical moments, by considering plots of global
envelopes, and by evaluating values of global envelope tests. This comparison
is not only done by looking at those tessellation characteristics used for
specifying the canonical sufficient statistic of the exponential model but
also for various other tessellation characteristics. To the best of our
knowledge, this is the first time that such a model selection procedure has
been used when analysing polycrystalline materials, cf. Šedivý et al. [2018]
and the references therein.

For the computations of Laguerre tessellations, we use the C++ library Voro++
[Rycroft, 2009]. This open source library enables computation of single
Laguerre cells which is an advantage for the estimation and MCMC simulation
procedures we use. For the calculation of functional summary statistics
(specifically $\hat{L}$, $\hat{F}$ and $\hat{G}$), we use the R-package
spatstat [Baddeley et al., 2015], and for global envelopes and corresponding
$p$-values we use the R-package GET [Myllymäki and Mrkvička, 2019]. We
implemented a code for pseudolikelihood estimation and the MCMC algorithms
used for simulation of the point process models and the models for the radii
given the points. This code is available at
https://github.com/VigoFierry/Lag_mod.

The paper is organised as follows. Section 2 provides details on Laguerre
tessellation data sets. Section 3 describes a particular data set which in
later sections is used to illustrate our statistical methodology. Section 4
specifies our hierarchical model construction and discusses how we used
maximum pseudolikelihood estimation and simulated under fitted models for
$\mathbf{x}_{n}$ and $\mathbf{r}_{n}\mid\mathbf{x}_{n}$. Section 5 details our
model selection procedure and applies it on the Laguerre tessellation data set
from Section 3. Finally, Section 6 contains concluding remarks. The paper is
accompanied by supplementary material available in appendices A-E in this
preprint.

##  2 Laguerre tessellation data sets

A Laguerre tessellation is defined as follows. Consider a marked point pattern
$\\{(x_{j},r_{j});j\in J\\}\subset\mathbb{R}^{3}\times(0,\infty)$ where each
$x_{j}\in\mathbb{R}^{3}$ is a spatial location with an associated mark
$r_{j}>0$ and the index set $J$ is finite or countable. Define the
corresponding point configuration $\mathbf{x}=\\{x_{j};j\in J\\}$ and the
associated (possibly infinite) vector of marks $\mathbf{r}=(r_{j};j\in J)$
(using some arbitrary ordering of $J$). Interpret $(x_{j},r_{j})$ as a closed
ball with center $x_{j}$ and radius $r_{j}$, and define the power distance
from a point $y\in\mathbb{R}^{3}$ to this ball by

$\rho(y,(x_{j},r_{j}))=\lVert y-x_{j}\rVert^{2}-r_{j}^{2}$

where $\|\cdot\|$ denotes Euclidean distance. Assuming that $\min_{j\in
J}\rho(y,(x_{j},r_{j}))$ exists for every $y\in\mathbb{R}^{3}$, the Laguerre
cell with generator $(x_{j},r_{j})$ is

$C(x_{j},r_{j}\mid\mathbf{x},\mathbf{r})=\bigcap_{k\in
J}\\{y\in\mathbb{R}^{3}\mid\rho(y,(x_{j},r_{j}))\leq\rho(y,(x_{k},r_{k}))\\}$

and the Laguerre tessellation generated by $(\mathbf{x},\mathbf{r})$ is the
collection of all non-empty cells $C(x_{j},r_{j}\mid\mathbf{x},\mathbf{r})$.

Notice that $\\{(x_{j},r_{j});j\in J\\}$ and $(\mathbf{x},\mathbf{r})$ are in
one-to-one correspondence. We use this identification without any mentioning
in the following. But it should be kept in mind when we abuse terminology by
calling $(\mathbf{x},\mathbf{r})$ a marked point pattern, and when we abuse
notation and write e.g. $((x,y,z),r)\in(\mathbf{x},\mathbf{r})$ (which means
that $((x,y,z),r)=(x_{j},r_{j})$ for some $j\in J$).

We refer to Lautensack and Zuyev [2008] and the references therein for
detailed studies of deterministic and random Laguerre tessellations. The
reason why Laguerre tessellations are so useful for modelling purposes may be
explained by a fundamental result due to Aurenhammer [1987] and extended by
Lautensack and Zuyev [2008] to infinite cases: any 3D (or higher-dimensional)
normal tessellation with convex cells is a Laguerre tessellation.

In Section 3 and later sections we consider the following setting. Let
$W=[0,a]\times[0,b]\times[0,c]$ be a rectangular parallelepiped and
$(\mathbf{x},\mathbf{r})$ a finite marked point pattern such that
$\mathbf{x}\subset W$. Define a periodic extension of
$(\mathbf{x},\mathbf{r})$ which is in accordance to $W$, i.e., the infinite
marked point pattern

$(\mathbf{x}^{*},\mathbf{r}^{*})=\bigcup_{((x,y,z),r)\in(\mathbf{x},\mathbf{r})}\bigcup_{(i,j,k)\in{\mathbb{Z}}^{3}}\\{((ia+x,jb+y,kc+z),r)\\}$
(1)

where $\mathbb{Z}^{3}$ is the 3D integer lattice. Note that
$\mathbf{x}=\mathbf{x}^{*}\cap W$. Then we consider the Laguerre tessellation
generated by $(\mathbf{x}^{*},\mathbf{r}^{*})$ and keep only those Laguerre
cells $C(x_{j},r_{j}\mid\mathbf{x}^{*},\mathbf{r}^{*})$ which are non-empty
and have $x_{j}\in\mathbf{x}$. Let $n$ be the number of such cells and
$(\mathbf{x}_{n},\mathbf{r}_{n})\subseteq(\mathbf{x},\mathbf{r})$ be the
marked point pattern specifying the generators of these cells. Then
$(\mathbf{x}_{n},\mathbf{r}_{n})$ is the Laguerre tessellation data set used
for our statistical analysis.

This construction of $(\mathbf{x}_{n},\mathbf{r}_{n})$ is one way of how
Laguerre tessellation data sets may be produced, and we will consider
$(\mathbf{x}_{n},\mathbf{r}_{n})$ as a realisation of a finite marked point
process [for background material on marked point processes, see e.g. Chiu et
al., 2013]. Note that the statistical methodology introduced in this paper can
immediately be extended to other settings with a finite marked point pattern
data set representing a Laguerre tessellation.

##  3 Application example

The Laguerre tessellation data set used throughout the text to illustrate our
statistical methodology is related to a study in Sedmák et al. [2016] of a
polycrystalline microstructure of a nickel titanium alloy. Using a notation as
in Section 2, we deal with a finite marked point pattern
$(\mathbf{x},\mathbf{r})$ extracted from a larger data set collected in
Petrich et al. [2019] by the so-called cross-entropy method applied to 3D
X-ray diffraction measurements. Specifically, $\mathbf{x}$ consists of the
2009 points which are contained in a 3D rectangular observation window $W$ of
size $40\times 40\times 85$ (the units being micrometers), where $W$ is a cut-
off from the entire material specimen. Some of the cells in the Laguerre
tessellation generated by $(\mathbf{x},\mathbf{r})$ are unbounded and hence
not contained in the material specimen. To account for this and for
computational reasons (i.e., the use of Voro++), we choose to apply a periodic
extension thereby obtaining our Laguerre tessellation data set
$(\mathbf{x}_{n},\mathbf{r}_{n})$, cf. Section 2. Here $n=1965$, so by
applying the periodic extension at most 44 cells are ‘lost’.

The left panel in Figure 1 below depicts $\mathbf{x}_{n}$ together with $W$
and A and B show images of the observed Laguerre tessellation together with
plots of the projections of $\mathbf{x}_{n}$ onto the $xy$, $xz$ and $yz$
planes, respectively. None of these indicate spatial inhomogeneity of
$\mathbf{x}_{n}$. So considering the periodic extension $\mathbf{x}^{*}_{n}$
(defined as in (1) but with $(\mathbf{x},\mathbf{r})$ replaced by
$(\mathbf{x}_{n},\mathbf{r}_{n})$), we find it reasonable to assume that the
distribution underlying $\mathbf{x}^{*}$ is invariant under translations in 3D
space. Equivalently, we assume that the distribution underlying
$\mathbf{x}_{n}$ is invariant under shifts when $W$ is wrapped on a 3D torus.
The right panel in Figure 1 shows the balls defined by
$(\mathbf{x}_{n},\mathbf{r}_{n})$. This indicates that the distribution
underlying the radii $\mathbf{r}_{n}$ conditioned on $\mathbf{x}_{n}$ can be
assumed to be homogeneous, i.e., we have invariance in this conditional
distribution when making shifts of $\mathbf{x}_{n}$ on the torus. This
assumption is further supported by plots in B showing the radii distributions
corresponding to a subdivision of $W$ into sets of equal size and shape.
Moreover, Figure 2 below shows a histogram of the radii together with a kernel
density estimate and the density of a beta distribution fitted by maximum
pseuodolikelihood estimation as discussed later in Sections 4 and 5.

Figure 1: The observation window $W$ and the point pattern $\mathbf{x}_{n}$
(left panel) and the balls defined by $(\mathbf{x}_{n},\mathbf{r}_{n})$ (right
panel). Figure 2: Histogram of the radii together with a kernel density
estimate (solid line) and a fitted density of a beta distribution (dashed
line) obtained by maximum pseudolikelihood estimation.

Figure 3 shows estimated/empirical functional summary statistics
$\hat{L}(t)-t$, $\hat{F}(t)$ and $\hat{G}(t)$ for the point pattern
$\mathbf{x}_{n}$ where we used the R-package spatstat [Baddeley et al., 2015]
for the calculations. For definitions and interpretations of these empirical
functions using edge correction factors, see e.g. Baddeley et al. [2015] (for
$\hat{L}$ we used Ripley’s isotropic edge correction factor and for $\hat{F}$
and $\hat{G}$ we used Kaplan-Meier edge correction factors). Figure 3 also
shows a concatenated 95%-global envelope/confidence region (the grey region)
under a homogeneous Poisson process (with the intensity estimated by $n$
divided by the volume of $W$), meaning that under this Poisson process all
three empirical functions are expected to be within the envelope with
probability 0.95. The envelopes were obtained using the R-package GET
[Myllymäki and Mrkvička, 2019] with 1999 simulations of the Poisson process
(increasing this to 9999 simulations did not change the results). The plots
show that the empirical functions are often outside the envelope, and the
corresponding $p$-value obtained by the global area rank envelope test
[Myllymäki et al., 2017] is below 0.1% [there are several possible choices of
making an envelope test, but the area rank envelope test is recommended in
Myllymäki and Mrkvička, 2020]. The way the functions based on the data differ
from the envelope indicates regularity in the point pattern $\mathbf{x}_{n}$.

Figure 3: From left to right, empirical functional summary statistics
$\hat{L}(t)-t$, $\hat{F}(t)$ and $\hat{G}(t)$ (solid lines) and simulated
95%-global envelope (grey regions) obtained under a fitted homogeneous Poisson
process. Dashed lines are averages of the simulated functional summary
statistics and the dots indicate when the empirical functions are outside the
envelope.

Figure 4 depicts histograms of some cell characteristics: volume of a cell
(abbreviated as ‘vol’); number of faces in a cell (‘nof’); sphericity of a
cell (‘spher’) as defined by

$\text{spher}=\frac{\pi^{1/3}(6\cdot\text{vol})^{2/3}}{\text{surf}}$

where ‘surf’ means surface area of a cell; and absolute difference in volume
‘dvol’ for two neighbouring cells (i.e., they share a face). Note that
$0<\text{spher}\leq 1$ with the value 1 corresponding to a sphere. Figure 4
shows that many cells are small, the distribution of nof is rather symmetric
and ranges from 4 (the smallest possible value) to 35, many cells are rather
spherical and many pairs of neighbouring cells have a range of differences in
cell sizes in the same order as the cell size itself.

Figure 4: From left to right, histograms of cell volume, number of faces,
sphericity and the absolute difference in volumes for neighbouring cells. The
unit for volume and difference in volume is $\mathrm{\mu}m^{3}$.

The upper right triangle in Table 1 shows the empirical correlations for the
cell characteristics vol, surf, nof and spher together with a new cell
characteristic, namely total edge length in a cell (‘tel’). The lower left
triangle in the table shows the empirical correlations for dvol and three
other face characteristics: area of a face (‘farea’), perimeter of a face
(‘fper’) and number of edges in a face (‘fnoe’). We see that all correlations
are high except for dvol which is nearly uncorrelated to any other face
characteristic. In particular vol, surf and tel are highly correlated, and
farea and fper are highly correlated (correlations $>0.9$).

vol | 0.971 | 0.938 | 0.841 | 0.680
---|---|---|---|---
farea | surf | 0.974 | 0.874 | 0.737
0.923 | fper | tel | 0.941 | 0.793
0.751 | 0.751 | fnoe | nof | 0.754
0.074 | 0.062 | 0.025 | dvol | spher
Table 1: Correlations of tessellation characteristics as defined in the text.

To see if there is dependence between $\mathbf{x}_{n}$ and $\mathbf{r}_{n}$,
we follow Stoyan et al. [2021] and consider the empirical mark correlation
function. We perform a permutation test based on the global area rank envelope
test, permuting 1000 times the marks/radii when $\mathbf{x}_{n}$ is fixed and
calculating the empirical mark correlation function each time. Figure 5 shows
the function based on the data: it falls outside the 95%-global envelope on
two intervals, and the corresponding $p$-value obtained by the global area
rank envelope test is 1.4%, so we reject the null hypothesis of independence
between $\mathbf{x}_{n}$ and $\mathbf{r}_{n}$.

Figure 5: The empirical mark correlation function (the solid line) and a
95%-global envelope (the grey region) obtained by permuting the radii when
fixing the points. The dashed line is the average of the simulated mark
correlation functions and the dots indicate when the empirical function is
outside the envelope.

##  4 Hierarchical model

This section details the two steps of our hierarchical model construction
which was briefly discussed in Section 1. To avoid confusion, we reserve the
notation $(\mathbf{x}_{n},\mathbf{r}_{n})$ for the data and use the notation
$(\mathbf{y}_{m},\mathbf{t}_{m})$ when we consider arguments of densities
where $\mathbf{y}_{m}=\\{y_{1},\ldots,y_{m}\\}\subset W$ is a finite point
configuration and $\mathbf{t}_{m}=(t_{1},\ldots,t_{m})\in(0,\infty)^{m}$ are
associated marks (if $m=0$ then $\mathbf{y}_{0}$ is the empty point
configuration and $\mathbf{t}_{0}$ can be ignored).

###  4.1 Point process models for $\mathbf{x}_{n}$

When modelling the point pattern data set $\mathbf{x}_{n}$ as a realization of
a spatial point process $\mathbf{X}$ on $W$, we refer to Møller and
Waagepetersen [2004] for background material on spatial point process models.
In accordance to our observations in Section 3, we assume that the
distribution of $\mathbf{X}$ is invariant under shifts when wrapping $W$ on a
3D torus, and a model for $\mathbf{X}$ should exhibit regularity. Pairwise
interaction point processes constitute a flexible class model for regularity.
Recall that $\mathbf{X}$ is a pairwise interaction point process on $W$ with
homogeneous interaction functions (with respect to shifts on $W$ wrapped on a
3D torus) if its distribution is given by a density with respect to the unit
rate Poisson process on $W$, where the density is of the form

$p(\mathbf{y}_{m})=\frac{1}{Z}\beta^{m}\prod_{i<j}\phi(\|y_{i}-y_{j}\|_{W})$

for $m=0,1,\ldots$ and every point configuration
$\mathbf{y}_{m}=\\{y_{1},\ldots,y_{m}\\}\subset W$. Here $\beta>0$ is a
parameter, $\|\cdot\|_{W}$ denotes shortest distance on $W$ wrapped on a 3D
torus, $\phi\geq 0$ is a so-called interaction function and $Z$ is a
normalising constant which depends on $\beta$ and $\phi$. In order to obtain a
well-defined density, a further condition on the interaction function is
needed; it suffices to assume that $0\leq\phi\leq 1$ in which case
$0<Z<\infty$.

As a parametric approximation of such a pairwise interaction point process we
consider a multiscale process [Penttinen, 1984]: for $d=1,2,\ldots$, let
$\mathcal{M}_{d}$ denote the model class given by densities of the form

$p(\mathbf{y}_{m})\propto\beta^{m}\Pi_{i=1}^{d-1}\gamma_{i}^{\sum_{j<k}\mathbb{I}_{[\delta_{i-1}<\|y_{j}-y_{k}\|_{W}\leq\delta_{i}]}}$
(2)

where $\beta>0$, $0\leq\gamma_{1}\leq 1$, $\ldots$, $0\leq\gamma_{d-1}\leq 1$
and $0<\delta_{1}<\ldots<\delta_{d-1}$ are unknown parameters. Here we have
omitted the normalising constant (which depends on all parameters),
$\mathbb{I}_{[\cdot]}$ denotes the indicator function and we set
$\delta_{0}=0$ and $0^{0}=1$. If $d=1$ we interpret the right hand side in (2)
as $\beta^{m}$. Note that $\mathcal{M}_{1}$ is just a homogeneous Poisson
process on $W$ with intensity $\beta$ and $\mathcal{M}_{2}$ is just a Strauss
process.

Simulation under the Poisson model $\mathcal{M}_{1}$ is well-known (see e.g.
Møller and Waagepetersen [2004]) and for simulation of the models
$\mathcal{M}_{d}$, $d\geq 2$, we used the birth-death-move Metropolis-Hastings
algorithm [Geyer and Møller, 1994, Møller and Waagepetersen, 2004].

###  4.2 Exponential tessellation models for $\mathbf{r}_{n}$ given
$\mathbf{x}_{n}$

When modelling the observed radii $\mathbf{r}_{n}$ as a realisation of an
$n$-dimensional vector, we condition on $\mathbf{X}=\mathbf{x}_{n}$ and
consider a conditional probability density function (pdf)
$p(r_{1},\ldots,r_{n}\mid x_{1},\ldots,x_{n})$ on $(0,\infty)^{n}$. More
precisely we assume $\mathbf{R}$ is a random vector (of random length) which
conditioned on $\mathbf{X}=\mathbf{y}_{m}$ has a conditional pdf
$p(\mathbf{t}_{m}\mid\mathbf{y}_{m})$ which is zero whenever
$C(y_{j},t_{j}\mid\mathbf{y}_{m}^{*},\mathbf{t}_{m}^{*})=\emptyset$ for some
$j\in\\{1,\ldots,m\\}$. Furthermore, in order to work with a well-defined
conditional pdf in (3) below we assume a mark space $\mathbb{M}=[0,6]$ so that
$p(\mathbf{t}_{m}\mid\mathbf{y}_{m})=0$ if
$\mathbf{t}_{m}\not\in\mathbb{M}^{m}$; this is in accordance to our data,
where $\min\mathbf{r}_{n}\approx 1.00$ and $\max\mathbf{r}_{n}\approx 5.15$.

We now give a general exponential family form of the conditional pdf where
$q\in\\{1,2,\ldots\\}$ is the dimension,
$\theta=(\theta_{1},\ldots,\theta_{q})$ denotes the canonical parameter and
$H=(H_{1},\ldots,H_{q})$ the canonical sufficient statistic: for
$\mathbf{t}_{m}\in\mathbb{M}^{m}$,

$p(\mathbf{t}_{m}\mid\mathbf{y}_{m})\propto\mathbb{I}_{[C(y_{j},t_{j}\mid\mathbf{y}_{m}^{*},\mathbf{t}_{m}^{*})\not=\emptyset,\,j=1,\ldots,m]}\exp\left(\sum_{i=1}^{q}\theta_{i}H_{i}(\mathbf{y}_{m},\mathbf{t}_{m})\right).$
(3)

The idea is to let each $H_{i}(\mathbf{y}_{m},\mathbf{t}_{m})$ depend on
either the radii, or tessellation characteristics of the cells
$C(y_{j},t_{j}\mid\mathbf{y}_{m}^{*},\mathbf{t}_{m}^{*})$ with $j=1,\ldots,m$
or interactions between these cells (here $\mathbf{y}_{m}^{*}$ is defined in a
similar way as $\mathbf{x}^{*}$, i.e., in the right hand side of (1)
$(\mathbf{x},\mathbf{r})$ is replaced by $(\mathbf{y}_{m},\mathbf{t}_{m})$).
Specifically, in Section 5 we consider the following cases (a)–(e), using
similar abbreviations as in Section 3, for short writing $H_{i}$ for
$H_{i}(\mathbf{y}_{m},\mathbf{t}_{m})$ and considering in (b)–(d) a sum over
the cells $C(y_{j},t_{j}\mid\mathbf{y}_{m}^{*},\mathbf{t}_{m}^{*})$,
$j=1,\ldots,m$, and in (e) a sum over all unordered pairs of cells sharing a
face:

  1. (a)

including both $H_{i}=\sum_{j=1}^{m}\log{\frac{t_{j}}{6}}$ and
$H_{i^{\prime}}=\sum_{j=1}^{m}\log{\left(1-\frac{t_{j}}{6}\right)}$ (with
$i\neq i^{\prime}$) somehow corresponds to a scaled beta distribution for the
radii if no other terms are included in (3) – ‘somehow’ because it is not
exactly a beta distribution since $p(\mathbf{t}_{m}\mid\mathbf{y}_{m})=0$ if
$C(y_{j},t_{j}\mid\mathbf{y}_{m}^{*},\mathbf{t}_{m}^{*})=\emptyset$ for some
$j\in\\{1,\ldots,m\\}$;

  2. (b)

$H_{i}=\sum{\mathrm{nof}}$ is twice the total number of faces;

  3. (c)

$H_{i}=\sum{\mathrm{surf}}$ is twice the total surface area of faces;

  4. (d)

$H_{i}=\sum{\mathrm{vol}}^{2}$ is the sum of squared volumes of cells;

  5. (e)

$H_{i}=\sum{\mathrm{dvol}}$ is the sum of difference in volumes of two cells
which share a face.

The density in (3) is then well-defined for all $\theta\in\mathbb{R}^{q}$
provided in case of (a) we have $\theta_{i}>-1$ and $\theta_{i^{\prime}}>-1$;
this follows since the mark space is bounded. Moreover, for simulation under
(3) we use a Metropolis within Gibbs algorithm where we alternate between
updating from the conditional densities

$p(t_{j}\mid t_{k},\,k\not=j,\,y_{1},\ldots,y_{m})\propto
p(t_{1},\ldots,t_{m}\mid y_{1},\ldots,y_{m}),\quad j=1,\ldots,m,$ (4)

using a Metropolis algorithm with a normal proposal.

###  4.3 Estimation

We assume that the unknown parameter $(\theta_{1},\ldots,\theta_{q})$ in (3)
varies independently of the parameters in (2). Then parameter estimation can
be done in two steps and the results will not depend on each other.

To avoid calculating the intractable normalising constants appearing in (2)
and (3), we use maximum pseudolikelihood methods [Besag, 1974, 1977, Besag et
al., 1982] in both cases where a few remarks are in order: For the definition
of pseudolikelihood function in connection to (2), we refer to Jensen and
Møller [1991]; see also D. For fixed $(\delta_{1},\ldots,\delta_{d})$, (2) is
an exponential family model with canonical parameter
$\vartheta=(\log\beta,\log\gamma_{1},\ldots,\log\gamma_{d})$ and so the log
pseudolikelihood function is a concave function of $\vartheta$ [Jensen and
Møller, 1991] and we can find a partial maximum pseudolikelihood estimate
(MPLE) of $\vartheta$ using the Newton-Raphson algorithm. Hence we obtain a
profile log pseudolikelihood which is maximised with respect to
$(\delta_{1},\ldots,\delta_{d})$ defined over a $d$-dimensional grid, thereby
providing the final MPLE. Furthermore, when considering (3) the
pseudolikelihood for $\theta$ is defined by the product of the conditional
densities of each $r_{j}$ given the $r_{k}$ with $k\not=j$, cf. (4). Since (3)
is also an exponential family, the log pseudolikelihood is concave and then
again Newton-Raphson can be used for finding the MPLE. Finally, we calculate
various integrals which appear in the pseudolikelihoods (in case of (2), see
Jensen and Møller [1991]; in case of (3), we consider the integral of the
right hand side of (4) with respect to $t_{j}$) by numerical methods.

##  5 Model selection and final joint model

This section details the model selection procedure briefly described in
Section 1.2 and which leads to our final joint model for $\mathbf{x}_{n}$ and
$\mathbf{r}_{n}$. The classical goodness of fit procedures are based on
measures given by maximizing likelihood functions and accounting for model
complexity by including a penalty in terms of the number of parameters (e.g.
AIC and BIC criteria). In Šedivý et al. [2018] the Vapnik-Chervonenkis theory
is used instead and a criterion based on residual sum of squares is developed
for partially ordered sets of deterministic tessellation models. Since we deal
with stochastic tessellation models but our likelihood functions are
intractable to maximize, we present another approach that corresponds to
modern trends in spatial statistics. It is based on how well geometrical
tessellation characteristics are described when comparing fitted models, where
we account for mutual correlations among the characteristics and we consider
the maximum of pseudolikelihood functions when comparing models with the same
number of parameters.

###  5.1 The fitted model for $\mathbf{x}_{n}$

For the point pattern $\mathbf{x}_{n}$ and $d=1,2,\ldots$, we select the first
model $\mathcal{M}_{d}$ which provides a satisfactory fit when considering
global envelopes and global area rank envelope test in the same way as in
Section 3. As noticed in Section 3, the Poisson model $\mathcal{M}_{1}$ is not
providing a satisfactory fit and as discussed in C this is also the case for
the Strauss model $\mathcal{M}_{2}$. The first model which is not significant
at level 5% is $\mathcal{M}_{3}$: Figure 6 shows that the empirical functional
summary statistics are within the 95%-global envelope produced in a similar
way as in Figure 3 but using simulations under the fitted multiscale process
with $d=3$. The corresponding $p$-value obtained by the global area rank
envelope test is $18.8\%$, and the maximum pseudolikelihood estimates (see
Section 4.3) are $\hat{\beta}=0.0168$, $\hat{\gamma}_{1}=0.5328$,
$\hat{\gamma}_{2}=0.8432$, $\hat{\delta}_{1}=1.25$ and
$\hat{\delta}_{2}=2.25$.

Figure 6: From left to right, empirical functional summary statistics
$\hat{L}(t)-t$, $\hat{F}(t)$ and $\hat{G}(t)$ (solid lines) and simulated
95%-global envelope (grey regions) obtained under a fitted multiscale process
with $d=3$. Dashed lines are averages of the simulated functional summary
statistics.

###  5.2 The fitted model for $\mathbf{r}_{n}$ conditioned on $\mathbf{x}_{n}$

Now, having fitted the model $\mathcal{M}_{3}$ for $\mathbf{x}_{n}$ it remains
to obtain a model for $\mathbf{r}_{n}$ conditioned on $\mathbf{x}_{n}$ where
we use a model selection procedure as follows. We consider a list of
tessellation characteristics

${\cal
L}=\\{\text{nof},\text{vol},\text{surf},\text{tel},\text{spher},\text{dvol}\\}$
(5)

when creating and evaluating more and more complex models as given by (3) and
(a)–(e) in Section 4.2. When comparing fitted models of the same dimension
$q$, we select the one with the highest value of the maximized log
pseudolikelihood function. The selected model is then evaluated by a global
area rank envelope test based on kernel smoothed densities of the empirical
distributions for the six tessellation characteristics in $\mathcal{L}$, using
$499$ simulations of the joint model for $\mathbf{x}_{n}$ and
$\mathbf{r}_{n}$. Here we concatenate the six densities and use the R-package
GET [Myllymäki and Mrkvička, 2019], and we simulate from the joint model
rather than from the conditional model of $\mathbf{r}_{n}$ given
$\mathbf{x}_{n}$, since our aim is to replace expensive laboratory experiments
with simulations from the joint model. Therefore, in Table 2, we also evaluate
our fitted models by comparing empirical (column ‘data’) and simulated means
and standard deviations of the tessellation characteristics in $\mathcal{L}$,
using 100 simulations of first $\mathbf{y}_{m}$ and second under the different
fitted models of $\mathbf{t}_{m}$ conditioned on $\mathbf{y}_{m}$. The table
also shows the (signed) difference between the empirical and the simulated
values divided by the empirical values (the values given in percentages) – we
refer to such a value as a deviation.

First, we consider the empirical distribution of the radii: As the histogram
in Figure 2 looks like a beta distribution, we first propose in (3) only to
include the terms in (a) so that $q=2$,
$H_{1}=\sum_{j=1}^{m}\log\tfrac{t_{j}}{6}$ and
$H_{2}=\sum_{j=1}^{m}\log\left(1-\tfrac{t_{j}}{6}\right)$ – we refer to this
model as ‘beta’. The column ‘beta’ in Table 2 shows that except for dvol the
means match the data well, and the deviations of standard deviations vary by
$12$ to $24$ percent. However, the fitted ‘beta’ model was highly significant
when evaluated by the global area rank envelope test.

|  | models |  |  | data
---|---|---|---|---|---
|  | beta | ${\mathrm{beta}}+{\mathrm{dvol}}$ | ${\mathrm{beta}}+{\mathrm{nof}}+{\mathrm{dvol}}$ |
nof | mean | 14.82 | 14.87 | 14.95 | 14.98
|  | -1% | -1% |  $<$-1% |
| sd | 5.58 | 5.38 | 5.20 | 4.92
|  | +13% | +9% | +6% |
vol | mean | 70.65 | 70.65 | 70.65 | 69.21
|  | +2% | +2% | +2% |
| sd | 65.65 | 63.42 | 62.01 | 58.89
|  | +12% | +8% | +5% |
surf | mean | 91.14 | 91.53 | 93.58 | 92.46
|  | -1% | -1% | +1% |
| sd | 57.26 | 54.07 | 51.50 | 47.86
|  | +20% | +13% | +8% |
tel | mean | 66.16 | 66.12 | 67.69 | 67.92
|  | -3% | -3% |  $<$-1% |
| sd | 31.99 | 30.86 | 29.52 | 27.38
|  | +17% | +13% | +8% |
spher | mean | 0.76 | 0.76 | 0.77 | 0.78
|  | -3% | -3% | -1% |
| sd | 0.11 | 0.10 | 0.10 | 0.087
|  | +24% | +19% | +16% |
dvol | mean | 79.87 | 73.72 | 69.82 | 68.89
|  | +16% | +7% | +1% |
| sd | 72.75 | 68.57 | 62.78 | 65.04
|  | +12% | +5% | -4% |
Table 2: Means and standard deviations of the tessellation characteristics
given in (5) and as obtained by simulations under various joint models and
from the data. Deviations are given in percentages. See the text for details.

Second, we expand the model by including one of the terms
$H_{3}=\sum{\mathrm{vol}}^{2}$, $\sum{\mathrm{nof}}$, $\sum{\mathrm{surf}}$ or
$\sum{\mathrm{dvol}}$ so that $q=3$. Here, we do not include
$\sum{\mathrm{vol}}$, since this is a constant; or $\sum{\mathrm{tel}}$, since
the correlation coefficient between vol and tel for each cell is close to 1;
or $\sum{\mathrm{spher}}$, since by definition spher is given by vol and surf
for each cell. Table 3 shows that the model with $H_{3}=\sum{\mathrm{dvol}}$
has the largest maximized log pseudolikelihood function. For this model,
comparing the columns ‘beta’ and ‘beta+dvol’ in Table 2 we obtain now better
results for dvol. Moreover, all standard deviations are now reduced, and the
$p$-value based on global area rank envelope test is 9.8%.

Third, we investigate the effect of expanding the model ‘beta’ with any two of
the four terms $\sum{\mathrm{vol}}^{2}$, $\sum{\mathrm{nof}}$,
$\sum{\mathrm{surf}}$ and $\sum{\mathrm{dvol}}$ so that $q=4$. Table 3 shows
that ‘beta+nof+dvol’ provides the best fit according to the maximized log
pseudolikelihood function. For this model, the $p$-value based on global area
rank envelope test is 10.6% which is slightly larger than the $p$-value of
9.8% for the model ‘beta+dvol’. Table 2 shows an improvement for both the mean
values and the standard deviations when comparing ‘beta+nof+dvol’ with
‘beta+dvol’. Figure 7 shows empirical kernel estimates of the densities for
the six characteristics in $\mathcal{L}$ together with a $95\%$-global
envelope obtained when concatenating all six empirical densities. The
empirical functions are completely covered by the envelope. Moreover, the
maximum pseudolikehood estimates are $\hat{\theta}_{1}=4.709$,
$\hat{\theta}_{2}=5.982$, $\hat{\theta}_{3}=-2.376\times 10^{-1}$ and
$\hat{\theta}_{4}=3.021\times 10^{-2}$. Thus under this fitted model
realisations become more likely as the total number of faces decreases or the
sum of differences in volumes between neighbouring cells increases (when all
other terms $H_{i}$ in (3) are fixed).

In a similar way, we also tried fitting models without including the terms in
(a), thereby obtaining ‘nof+surf+dvol’ as the final model for $q=3$ and where
the $p$-value is 7.7%. A plot similar to Figure 7 but for the fitted
‘nof+surf+dvol’ model is given in E.

model |  $\log$PL$(\hat{\theta})$
---|---
${\mathrm{beta}}$ | -2532.45
${\mathrm{beta}}+{\mathrm{vol}}^{2}$ | -2866.39
${\mathrm{beta}}+{\mathrm{nof}}$ | -2594.72
${\mathrm{beta}}+{\mathrm{surf}}$ | -2483.98
${\mathrm{beta}}+{\mathrm{dvol}}$ | -2468.76
${\mathrm{beta}}+{\mathrm{vol}}^{2}+{\mathrm{nof}}$ | -2513.89
${\mathrm{beta}}+{\mathrm{vol}}^{2}+{\mathrm{surf}}$ | -2823.71
${\mathrm{beta}}+{\mathrm{vol}}^{2}+{\mathrm{dvol}}$ | -2714.10
${\mathrm{beta}}+{\mathrm{nof}}+{\mathrm{surf}}$ | -2498.58
${\mathrm{beta}}+{\mathrm{nof}}+{\mathrm{dvol}}$ | -2465.12
${\mathrm{beta}}+{\mathrm{surf}}+{\mathrm{dvol}}$ | -2477.63
Table 3:  Maximized log pseudolikelihood functions for radii models
conditioned on the points. The value for the ‘beta’ model is given as a
reference.  Figure 7:  Estimated densities of tessellation characteristics,
namely nof (a), vol (b), surf (c), tel (d), spher (e) and dvol (f). The solid
lines are the functions based on the data $(\mathbf{x}_{n},\mathbf{r}_{n})$
and the grey regions are 95%-global envelopes under the fitted model
‘beta+nof+dvol’. Dashed lines are averages of the simulated densities.

##  6 Concluding remarks

Our hierarchical model construction for Laguerre tessellation data sets, using
a Gibbs model for the points and a conditional model for the radii given the
points, may be useful in several other contexts, including 2D cases. As
compared to the use of Gibbs models in Dereudre and Lavancier [2011] and Seitl
et al. [2021], our hierarchical approach is less time demanding: To simulate a
vast number of standard point process models is a routine task which is not
too much time demanding. The more exhaustive part is to simulate the radii
since Laguerre cells have to be recomputed, but still in comparison with Seitl
et al. [2021] it becomes faster in our case since we condition on the points.

For the Gibbs point process model, we considered multiscale models which
approximate the class of pairwise interaction point processes and lead to a
natural model selection procedure by considering an increasing number of
interaction parameters. We might consider point process models with higher
order interactions such as Geyer’s triplet interaction process [Geyer, 1999]
or the area interaction process involving interactions of all orders [Baddeley
and Van Lieshout, 1995]. However, the computational complexity increases with
more complex models.

For the conditional distribution of radii given points, we have demonstrated
the usefulness of exponential family models with a canonical sufficient
statistic based on tessellation characteristics. As there are several
possibilities for selecting tessellation characteristics and there is no
natural ordering (as there is in the case of the multiscale model), a model
selection procedure is less straightforward. However, we demonstrated how to
compare fitted models of the same dimension by considering maximized log
pseudolikelihood functions and how to evaluate selected fitted models by
partly comparing moment properties of tessellation characteristics under
simulations from the model with empirical moments and partly by considering
plots of global envelopes and by evaluating their corresponding $p$-values.

When calculating pseudolikelihood functions, various integrals have to be
evaluated, cf. Section 4.3. We used simple grid-based approximation methods
for both the pseudolikelihood based on the points and that based on the radii
given the points. A small simulation study indicated that this may cause some
bias in the MPLEs (also for MPLEs of the Gibbs models in Dereudre and
Lavancier [2011] and Seitl et al. [2021] some bias appeared). Baddeley and
Turner [2014] noticed that bias of MPLEs for spatial point process models
usually does not reflect weaknesses of the MPLE method, but is probably due to
the effects of discretization of the observation window. For 2D spatial point
process models, efficient numerical methods have been developed [Baddeley and
Turner, 2014, Baddeley et al., 2014] and spatstat provides useful software,
but to the best of our knowledge it remains to make progress for in general 3D
point processes and for the special case of our model. We leave this important
and huge task for future research.

## Acknowledgements

The research was supported by the Czech Science Foundation, project 19-04412S,
and by The Danish Council for Independent Research — Natural Sciences, grant
DFF – 7014-00074 ‘Statistics for point processes in space and beyond’. We
thank the referees for useful comments and suggestions.

##  Appendix A Visualization of the observed Laguerre tessellation

In order to visualize the Laguerre tessellation data set described in Section
3, Figure 8 shows four equidistant slices perpendicular to the $z$ axis.

Figure 8: Four slices in $xy$ plane of the observed Laguerre tessellation.

##  Appendix B Homogeneity of points and radii

The key assumption of homogeneity of $\mathbf{x}_{n}$ and $\mathbf{r}_{n}$
conditioned on $\mathbf{x}_{n}$ is investigated in this section. Figure 9
depicts $\mathbf{x}_{n}$ together with the observation window and the
projections of $\mathbf{x}_{n}$ onto $xy$, $xz$ and $yz$ planes. The figure
indicates spatial homogeneity of $\mathbf{x}_{n}$. In Figure 10, the left
panel shows a kernel density estimate of the radii distribution and, for
comparison purposes, the right panel shows eight kernel density estimates of
the radii distributions corresponding to a subdivision of $W$ into the
$2^{3}=8$ subsets obtained by dividing the three sides of the rectangular
region $W$ into half sides. For the number of points per subset, the mean is
$232$, the minimum is $219$ and the maximum is $270$. The similarity of the
eight density estimates is in accordance with the assumption that the
distribution of $\mathbf{r}_{n}$ conditioned on $\mathbf{x}_{n}$ is
homogeneous.

.

Figure 9: The spatial point pattern $\mathbf{x}_{n}$ of the generators of the
Laguerre tessellation after applying a periodical boundary condition together
with its projections onto the $xy$, $xz$ and $yz$ planes.

Figure 10: A kernel density estimate based on the radii $\mathbf{r}_{n}$ (left
panel) and eight kernel density estimates based on the radii associated to the
points in the eight sets for the subdivision of $W$ obtained by dividing its
sides into half parts (right panel).

##  Appendix C Strauss process

Figure 11 shows three empirical functional summary statistics together with
$95\%$-global envelopes obtained by simulations under the Strauss process
model $\mathcal{M}_{2}$ (the multiscale process with $d=2$) fitted by maximum
pseudolikelihood. The corresponding $p$-value obtained by the global area rank
envelope test is $4.4\%$.

Figure 11: From left to right, empirical functional summary statistics
$\hat{L}(t)-t$, $\hat{F}(t)$ and $\hat{G}(t)$ (solid lines) and simulated
$95\%$-global envelope (grey region) obtained under a fitted Strauss process.
Dashed lines are averages of the simulated functional summary statistics.

##  Appendix D Pseudolikelihood functions

###  D.1 Pseudolikelihood for the points

As in the paper, consider a point process with a given realization
$\mathbf{y}_{m}$ in the observation window $W$ and modelled by a parametric
density $p_{\vartheta}$ (with respect to the unit rate Poisson process) where
$\vartheta$ is an unknown real parameter vector. Assume $p_{\vartheta}$ is
hereditary, i.e., for $u\in W\setminus\mathbf{y}_{m}$ we have
$p(\\{u\\}\cup\mathbf{y}_{m})>0$ whenever $p(\mathbf{y}_{m})>0$. Then the
pseudolikelihood function for the points is

${\cal{PL}}(\vartheta)=\exp{\left(|W|-\int_{W}\lambda^{\star}_{\vartheta}(u,\mathbf{y}_{m})\,\mathrm{d}u\right)}\prod_{y_{j}\in\mathbf{y}_{m}}\lambda^{\star}_{\vartheta}(y_{j},\mathbf{y}_{m}\setminus
y_{j})$

where
$\lambda^{\star}_{\vartheta}(u,\mathbf{y}_{m})={p_{\vartheta}(u\cup\mathbf{y}_{m})}/{p_{\vartheta}(\mathbf{y}_{m})}$
is the Papangelou conditional intensity (setting $0/0=0$). Denoting the
natural logarithm by $\log$, the log pseudolikelihood function in the case of
the multiscale process ${\cal M}_{d}$ (see (2) with $d>1$) becomes

$\begin{split}\log{\cal{PL}}(\beta,\gamma)=&|W|-\int_{W}\beta\prod_{i=1}^{d-1}\gamma_{i}^{t_{\delta_{i}}(u,\mathbf{y}_{m})-t_{\delta_{i-1}}(u,\mathbf{y}_{m})}\,\mathrm{d}u\\\
&+m\log\beta+2\sum_{i=1}^{d-1}\left(S_{\delta_{i}}(\mathbf{y}_{m})-S_{\delta_{i-1}}(\mathbf{y}_{m})\right)\log\gamma_{i},\\\
\end{split}$

where
$t_{\delta}(u,\mathbf{y}_{m})=\sum_{y_{j}\in\mathbf{y}_{m}}\mathbb{I}_{\\{0<\|u,y_{j}\|\leq\delta\\}}$
and
$S_{\delta}(\mathbf{y}_{m})=\sum_{u\in\mathbf{y}_{m}}t_{\delta}(u,\mathbf{y}_{m})/2$.

###  D.2 Pseudolikelihood for the radii given the points

Denoting the conditional density (3) by $p_{\theta}$, the pseudolikelihood
function for the radii $\mathbf{t}_{m}$ given the points $\mathbf{y}_{m}$ is

${\cal{PL}}(\theta)=\prod_{j=1}^{n}p_{\theta}(t_{j}\mid
y_{j},(y_{k},t_{k})\text{ with }k\neq
j)=\prod_{j=1}^{n}\frac{h_{\theta}(\mathbf{t}_{m}\mid\mathbf{y}_{m})}{\int_{0}^{6}h_{\theta}(\mathbf{t}_{m}\mid\mathbf{y}_{m})\,\mathrm{d}t_{j}}$

where $p_{\theta}(y_{m},\mathbf{t}_{m})\propto
h_{\theta}(\mathbf{t}_{m}\mid\mathbf{y}_{m})$. Assume that
$(\mathbf{y}_{m},\mathbf{t}_{m})$ is feasible, i.e., all Laguerre cells
$C(y_{j},t_{j}\mid\mathbf{y}_{m}^{*},\mathbf{t}_{m}^{*})$ are nonempty; we
denote this property by
$(\mathbf{y}_{m},\mathbf{t}_{m})\notin{\cal{L}}(\emptyset)$. Then the log
pseudolikelihood function for the radii distribution given the points is

$\begin{split}\log{\cal{PL}}(\theta)=&n\sum_{i=1}^{q}\theta_{i}H_{i}(\mathbf{y}_{m},\mathbf{t}_{m})\\\
&-\sum_{j=1}^{m}\log\int_{0}^{6}\mathbb{I}_{\\{(\mathbf{y}_{m},\mathbf{t}_{m}^{j,u})\notin{\cal{L}}(\emptyset)\\}}\exp\left(\sum_{i=1}^{q}\theta_{i}H_{i}(\mathbf{y}_{m},\mathbf{t}_{m}^{j,u})\right)\,\mathrm{d}u\\\
\end{split}$

where $\mathbf{t}_{m}^{j,u}=(t_{1},\ldots,t_{j-1},u,t_{j+1},\ldots,t_{m})$.

##  Appendix E A model for the radii given the points where the ‘beta’ term
does not appear

In order to see the effect of omitting the ‘beta’ term in the density (3) for
the radii conditioned on the points, we fitted the model ‘nof$+$surf$+$dvol’
corresponding to
$H=(\sum{\mathrm{nof}},\sum{\mathrm{surf}},\sum{\mathrm{dvol}})$ in (3). Since
this is a three-dimensional model, we compare with the results for the fitted
three-dimensional model ‘beta+dvol’ in the paper.

Figure 12 shows empirical kernel estimates of the densities for the six
characteristics in $\mathcal{L}$ (see (5)) together with a $95\%$-global
envelope obtained when concatenating all six empirical densities and using the
R-package GET with $499$ simulations of the fitted joint model for
$\mathbf{x}_{n}$ and $\mathbf{r}_{n}$. The plots show that the functions based
on the data are completely covered by the envelope. Moreover, the
corresponding $p$-value obtained by the global area rank envelope test is
7.7%, which is less than the $p$-value 9.8% obtained for the fitted
‘beta+dvol’ model. Finally, the maximum pseudolikehood estimates are
$\hat{\theta}_{1}=3.209\times 10^{-2}$, $\hat{\theta}_{2}=-2.993\times
10^{-2}$ and $\hat{\theta}_{3}=2.049\times 10^{-4}$, and the value of the
maximized log pseudolikelihood function is $-2470.12$, which is smaller than
the corresponding value of $-2468.76$ for the fitted ‘beta+dvol’ model.

Figure 12:  Estimated densities of tessellation characteristics, namely nof
(a), vol (b), surf (c), tel (d), spher (e) and dvol (f). The solid lines are
the functions based on the data $(\mathbf{x}_{n},\mathbf{r}_{n})$ and the grey
regions are 95%-global envelopes under the fitted model ‘nof+surf+dvol’.
Dashed lines are averages of the simulated densities.

## References

  * Aurenhammer [1987] Aurenhammer, F., 1987.  A criterion for the affine equivalence of cell complexes in $\mathbb{R}^{d}$ and convex polyhedra in $\mathbb{R}^{d+1}$.  Discrete Computational Geometry 2, 49–64.
  * Baddeley et al. [2014] Baddeley, A., Coeurjolly, J.F., Rubak, E., Waagepetersen, R., 2014\.  Logistic regression for spatial Gibbs point processes.  Biometrika 101, 377–392.
  * Baddeley et al. [2015] Baddeley, A., Rubak, E., Turner, R., 2015.  Spatial Point Patterns: Methodology and Applications with R.  Chapman & Hall/CRC, Boca Raton.
  * Baddeley and Turner [2014] Baddeley, A., Turner, R., 2014\.  Bias correction for parameter estimates of spatial point process models.  Journal of Statistical Computation and Simulation 84, 1621–1643.
  * Baddeley and Van Lieshout [1995] Baddeley, A.J., Van Lieshout, M.N.M., 1995\.  Area-interaction point processes.  Annals of the Institute of Statistical Mathematics 47, 601–619.
  * Besag [1974] Besag, J.E., 1974.  Spatial interaction and the statistical analysis of lattice systems (with discussion).  Journal of the Royal Statistical Society Series B 36, 192–236.
  * Besag [1977] Besag, J.E., 1977.  Some methods of statistical analysis for spatial data.  Bulletin of the International Statistical Institute 47, 77–92.
  * Besag et al. [1982] Besag, J.E., Milne, R.K., Zachary, S., 1982.  Point process limits of lattice processes.  Journal of Applied Probability 19, 210–216.
  * Chiu et al. [2013] Chiu, S.N., Stoyan, D., Kendall, W.S., Mecke, J., 2013\.  Stochastic Geometry and Its Applications.  John Wiley & Sons, Chichester.
  * Dereudre and Lavancier [2011] Dereudre, D., Lavancier, F., 2011\.  Practical simulation and estimation for Gibbs Delaunay–Voronoi tessellations with geometric hardcore interaction.  Computational Statistics & Data Analysis 55, 498–519.
  * Geyer [1999] Geyer, C.J., 1999.  Likelihood inference for spatial point processes, in: Barndorff-Nielsen, O.E., Kendall, W.S., van Lieshout, M.N.M. (Eds.), Stochastic Geometry: Likelihood and Computation. Chapman & Hall/CRC, Boca Raton, pp. 141–172.
  * Geyer and Møller [1994] Geyer, C.J., Møller, J., 1994\.  Simulation procedures and likelihood inference for spatial point processes.  Scandinavian Journal of Statistics 21, 359–373.
  * Jensen and Møller [1991] Jensen, J.L., Møller, J., 1991\.  Pseudolikelihood for exponential family models of spatial point processes.  The Annals of Applied Probability 1, 445–461.
  * Kuhn et al. [2020] Kuhn, J., Schneider, M., Sonnweber-Ribic, P., Boelke, T., 2020\.  Fast methods for computing centroidal Laguerre tessellations for prescribed volume fractions with applications to microstructure generation to polycrystalline materials.  Computer Methods in Applied Mechanics and Engineering 369, 113–175.
  * Lautensack [2008] Lautensack, C., 2008.  Fitting three-dimensional Laguerre tessellations to foam structures.  Journal of Applied Statistics 35, 985–995.
  * Lautensack and Zuyev [2008] Lautensack, C., Zuyev, S., 2008\.  Random Laguerre tessellations.  Advances in Applied Probability 40, 630–650.
  * Lyckegaard et al. [2011] Lyckegaard, A., Lauridsen, E., Ludwig, W., Fonda, R., Poulsen, H., 2011.  On the use of Laguerre tessellations for representations of 3D grain structures.  Advanced Engineering Materials 13, 165–170.
  * Møller and Waagepetersen [2004] Møller, J., Waagepetersen, R.P., 2004\.  Statistical Inference and Simulation for Spatial Point Processes.  Chapman & Hall/CRC, Boca Raton.
  * Myllymäki and Mrkvička [2019] Myllymäki, M., Mrkvička, T., 2019\.  GET: Global Envelopes in R.  Avaliable at arXiv:1911.06583.
  * Myllymäki and Mrkvička [2020] Myllymäki, M., Mrkvička, T., 2020\.  Comparison of non-parametric global envelopes.  Avaliable at arXiv:2008.09650.
  * Myllymäki et al. [2017] Myllymäki, M., Mrkvička, T., Seijo, H., Grabarnik, P., Hahn, U., 2017.  Global envelope tests for spatial processes.  Journal of the Royal Society of Statistics Series B 79, 381–404.
  * Penttinen [1984] Penttinen, A., 1984.  Modelling Interactions in Spatial Point Patterns: Parameter Estimation by the Maximum Likelihood Method.  Number 7 in Jyväskylä Studies in Computer Science, Economics and Statistics, University of Jyväskylä.
  * Petrich et al. [2019] Petrich, L., Staněk, J., Wang, M., Westhoff, D., Heller, L., Šittner, P., Krill, III, C., Beneš, V., Schmidt, V., 2019.  Reconstruction of grains in polycrystalline materials from incomplete data using Laguerre tessellations.  Microscopy and Microanalysis 25, 743–752.
  * Quey and Renversade [2018] Quey, R., Renversade, L., 2018\.  Optimal polyhedral description of 3D polycrystals: Method and application to statistical and synchrotron X-ray diffraction data.  Computational Methods in Applied Mechanical Enginering 330, 308–333.
  * Rycroft [2009] Rycroft, C.H., 2009.  A three-dimensional Voronoi cell library in C++.  Chaos 19, 041111\.
  * Sedmák et al. [2016] Sedmák, P., Pilch, J., Heller, L., Kopeček, J., Wright, J., Sedlák, P., Frost, M., Šittner, P., 2016\.  Grain-resolved analysis of localized deformation in nickel-titanium wire under tensile load.  Science 353, 559–562.
  * Seitl et al. [2021] Seitl, F., Petrich, L., Staněk, J., Krill, C.E., Schmidt, V., Beneš, V., 2021\.  Exploration of Gibbs-Laguerre tessellations for three-dimensional stochastic modeling.  Methodology and Computing in Applied Probability, 23, 669–693.
  * Spettl et al. [2016] Spettl, A., Brereton, T., Duan, Q., Werz, T., Krill III, C.E., Kroese, D.P., Schmidt, V., 2016.  Fitting Laguerre tessellation approximations to tomographic image data.  Philosophical Magazine 96, 166–189.
  * Spettl et al. [2015] Spettl, A., Wimmer, R., Werz, T., Heinze, M., Odenbach, S., Krill III, C.E., Schmidt, V., 2015.  Stochastic 3D modeling of Ostwald ripening at ultra-high volume fractions of the coarsening phase.  Modelling and Simulation in Materials Science and Engineering 23, 065001.
  * Stoyan et al. [2021] Stoyan, D., Beneš, V., Seitl, F., 2021.  Dependent radius marks of Laguerre tessellations: a case study.  Australian & New Zealand Journal of Statistics 63, 19–32.
  * Šedivý et al. [2018] Šedivý, O., Westhoff, D., Kopeček, J., Krill III, C.E., Schmidt, V., 2018.  Data-driven selection of tessellation models describing polycrystalline microstructures.  Journal of Statistical Physics 172, 1223–1246.


# Representation Based Complexity Measures for Predicting Generalization in
Deep Learning

Parth Natekar
Department of Engineering Design
Indian Institute of Technology, Madras
<EMAIL_ADDRESS>
Manik Sharma
Department of Engineering Design
Indian Institute of Technology, Madras
<EMAIL_ADDRESS>

###### Abstract

Deep Neural Networks can generalize despite being significantly
overparametrized. Recent research has tried to examine this phenomenon from
various viewpoints and to provide bounds on the generalization error or
measures predictive of the generalization gap based on these viewpoints, such
as norm-based, PAC-bayes based, and margin-based analysis. In this work, we
provide an interpretation of generalization from the perspective of quality of
internal representations of deep neural networks, based on neuroscientific
theories of how the human visual system creates invariant and untangled object
representations. Instead of providing theoretical bounds, we demonstrate
practical complexity measures which can be computed ad-hoc to uncover
generalization behaviour in deep models. We also provide a detailed
description of our solution that won the NeurIPS competition on Predicting
Generalization in Deep Learning held at NeurIPS 2020. An implementation of our
solution is available at https://github.com/parthnatekar/pgdl.

##  1 Introduction

The generalization ability of deep neural networks despite having more
parameters than training data is a property of great interest from both
theoretical and application perpectives. Understanding the factors behind
generalization and possibly creating practical measures which are predictive
of generalization can lead to better designs as well as a better understanding
of this phenomanon.

Multiple hypothesis have been presented before based on various properties of
deep neural networks, from optimization based theories, to empirical measures
based on sharpness, to theoretically motivated generalization bounds on the
test error based on some intrinsic property of deep networks.

Theoretically motivated techniques which aim to relate generalization bounds
to measures of capacity such as VC Dimension and Rademacher Complexity are
vacuous when the network falls under the ‘over-parameterized’ regime. Bartlett
[3] provide a bound based on the norm of the weights, and Neyshabur et al.
(2015) [20] extend this to other norm based measures. Multiple studies have
related generalization to PAC-Bayes based bounds [8, 16]. Optimization based
theories rely on the presence of implicit regularization via SGD or explicit
regularization via dropout or weight decay [19, 23]. However, Zhang et al.
[25] showed that despite such implicit and explicit regularization, deep
neural networks can achieve negligible training error on randomly labelled
samples.

Empirical measures have also been proposed before. Keskar et al. [15] connect
generalization to flatness of minima, Arora et al. [2] derive a bound based on
the noise attenuation properties of deep networks, while Morcos et al. [18]
hypothesize that networks’ reliance on single directions is a good predictor
of generalization performance, and demonstrate this via ablation studies and
addition of noise. Jiang et al. (2018) [12] demonstrate that the margin
distribution based measures show a high correlation with generalization gap.

Many of the theoretical measures discussed above are upper bounds and are not
tight enough to use in an empirical setting. Measures based on the weight
norms cannot be used when there are variations in network topology. Other
measures like PAC-Bayes or empirical measures like sharpness are stochastic
and it is not clear how these scale with changes in network size and other
network parameters. Thus, most of these measures do not always provide
intuitive support to practitioners in terms of designing models which
generalize better or comparing generalization ability of already trained
networks with limited access to their hyperparameters and design choices.

In this work, we provide an interpretation of generalization ability in terms
of the quality of representations deep neural networks. We propose simple and
intuitive complexity measures that can be computed in a post-hoc setting, i.e.
when comparing models after they have already been trained. Our measures are
predictive of generalization irrespective of changes in hyperparameters and
model architectures, and achieve high scores on the final and development
tasks of the NeurIPS competition on Predicting Generalization in Deep
Learning.

##  2 Quality of Representations and Neuroscientific Parallels

Representation learning ascribes the effectiveness of deep learning to the
ability of deep networks to learn representations that are invariant to
nuisances in the input such as translations, rotations, occlusions, and also
“disentangled”, that is separating out independent factors in the high-
dimensional space of data [1]. Bengio et al. [4] discuss, broadly, the
properties of good feature representations and discuss how general priors
about the world around us can be reflected through these. On a high level, it
makes sense that deep models which build ’better’ feature representations are
more likely to generalize better. We posit that certain properties of learned
representations in deep neural networks improve generalization ability.
However, not all properties that are hypothesized to be linked with
representation quality may be positively correlated with generalization, and
if the are, may not be easily computable. Moreover, while it is intuitively
understood and often taken for granted that representation quality and
generalization are linked, there have been limited studies which explicitly
relate the two and empirically investigate this relationship.

Neuroscientific and psychological studies have shown that humans rely on
invariant and consistent feature representations irrespective of controlled
naturalistic variations of objects [14, 21]. Studies have shown that one of
the key aspects of invariant object detection is the ability of the visual
hierarchy to create easily separable representation manifolds [7]. Robustness
is another related and well-known quality of human vision that helps in
invariant object detection. Human vision has been shown to be more robust than
deep neural networks on image manipulations like contrast reduction, additive
noise or eidolon-distortions [10].

While it is unclear if deep feed-forward convolutional neural networks use
similar strategies, studies have shown a hierarchical correspondence between
deep network representations and spatio-temporal human visual brain
representations [5]. Thus, creating measurable metrics that quantify the
properties of representations discussed above can potentially predict
generalization in deep learning.

Based on such neuroscientific arguments and the above ideas from
representation learning, we choose three metrics of quality of representations
that have the potential to predict generalization in a post-hoc diagnostic
setting. These are: (i) Consistency, (ii) Robustness (to controlled
variations), and (iii) Separability.

In the next sections, we detail simple complexity measures based on these
metrics of quality of representations that allow us to compare models that
have been trained with a range of different hyperparameters and architectures.
We report the performance of our complexity measures on a range of datasets
and models which are part of the NeurIPS competition on Predicting
Generalization in Deep Learning (PGDL). These complexity measures achieve top
scores on the development and final phases of PGDL. We report our scores on
the public and private datasets of PGDL, along with scores of the different
techniques we implemented that were not submitted for evaluation.

##  3 Preliminaries and Notation

First, we define some notation. A feedforward neural network is denoted by
$f_{w}:X\rightarrow\mathbb{R}^{\kappa}$ and its weight parameters by $w$. The
weight tensor of the $i^{th}$ layer of the network is denoted by $W_{i}$, so
that $w=$vec$(\mathcal{W}_{1},...,\mathcal{W}_{d})$, where d is the depth of
the network, and vec represents the vectorization operator. Let
$\mathcal{A}=(\mathcal{A}_{1},...,\mathcal{A}_{d})$ be the feature maps of the
network at each layer, and $f_{w,k}:A_{k}\rightarrow\mathbb{R}^{\kappa}$ the
intermediate model that maps representations at layer $k$ to the network
output. Let $\mathcal{D}$ be the data distribution over inputs and outputs.
Let $S=\\{X,y\\}$ be the training set, where $X=\\{X_{i}\\}_{i=1}^{N}$ are the
data points and $y=\\{y_{i}\\}_{i=1}^{N}\in\\{1,...,\kappa\\}$ are the
corresponding training labels sampled randomly from $\mathcal{D}$. We denote
by $f_{w}(X)[j]$ be the $j^{th}$ output of the network. The empirical 0-1
classification loss $L$ is defined as
$\hat{L}=\frac{1}{m}\sum_{i=1}^{m}I\big{(}f_{w}(\tilde{X}_{i})[y_{i}]\leq\max_{j\neq\tilde{y}_{i}}f_{w}(\tilde{X}_{i})[j]\big{)}$.
Finally, the generalization gap is defined as the difference between the
training and test accuracy of the network.

Our task is to find a complexity measure $\mathcal{C}$ that is predictive of
the generalization gap using only the network weights and the input dataset in
a post-hoc setting, i.e. after models have been trained. Such a measure should
allow a relative comparison between models trained using a range of different
hyperparameters and architectures, thus capturing properties of deep neural
networks that are responsible of generalization irrespective of how they are
trained.

In order to better capture relationships between complexity measures and the
generalization gap and to separate these from spurious correlations, a
conditional mutual information based metric is used as in Jiang et al. (2019)
[13].

##  4 Representation Based Complexity Measures

###  4.1 Consistency of Representations

It is logical that networks which build consistent representations of similar
objects are good at generalizing. In such a case, it is likely that the
network is recognizing the underlying causal factors that define that object,
which allows it to create consistent representations in-spite of variations.

Clustering of feature representations to encourage compactness and
interpretability and to create parsimonious representations has been explored
before [17]. Clustering is also used in the context of network compression
[22]. In the context of generalization, measuring how clustered
representations of a particular class are tell us about their consistency and
compactness. Models where feature representations at layer $L_{i}$ for each
class lie in their own cluster and are easily distinguishable from other class
clusters in high dimensional feature spaces are likely to have similar
representations at subsequent layers $L_{j|j>i}$, and hence likely to be
predicted to be in the same class. Figure 1 provides a visual depiction of
this idea.

Therefore, we require a measure that tells us about the quality of the
clustering or compactness of feature representations in intermediate layers of
the model, given the ground-truth label as the cluster index. However, simply
measuring the compactness of representations using distance-based measures
does not suffice when comparing models with different topology. For example,
one possible metric could be the expectation over all classes of the average
distance of representations of a class to the cluster centroid of that class.
The smaller the average distance to cluster centroids for each class, the more
consistent the representations. However, this is likely to fail when models
have different shapes, and comparisons happen across different vector spaces.
We use the Davies-Bouldin Index to measure clustering quality. The Davies-
Bouldin Index measures the ratio of within-cluster scatter to between cluster
separation. Being a ratio, it is easily comparable across models with
different architectures. Measuring the Davies-Bouldin index on intermediate
feature representations tells us how consistent representations are within a
class and how different they are from other classes. Mathematically, for a
given layer $L_{k}$ and its activations $A_{k}$,

$\tilde{\mathcal{A}}_{k}=\Phi(\mathcal{A}_{k})$ (1)
$\mathcal{S}_{i}=\Big{(}\frac{1}{n_{i}}\sum_{1}^{n_{i}}|\tilde{\mathcal{A}}_{k}^{i}-\mu_{\tilde{\mathcal{A}}_{k}^{i}}|^{p}\Big{)}^{1/p}$
(2)
$\mathcal{M}_{i,j}=||\mu_{\tilde{\mathcal{A}}_{k}^{i}}-\mu_{\tilde{\mathcal{A}}_{k}^{j}}||_{p}$
(3)

Where, $\Phi$ indicates a dimensionality reduction operation, $i$ and $j$ are
two different classes, $A_{k}^{i}$ are the feature representations of samples
belonging to class $i$, $\mu_{\tilde{\mathcal{A}}_{k}^{i}}$ is the cluster
centroid of the representations of class $i$, $S_{i}$ is a measure of scatter
within representations of class $i$, and $\mathcal{M}_{i,j}$ is a measure of
separation between representations of classes $i$ and $j$. Then, the
complexity measure based on the Davies Bouldin Index for representations at
layer $k$ is calculated as:

$\mathcal{C}_{k}=\frac{1}{\kappa}\sum_{i=1}^{\kappa}\max_{i\neq
j}\frac{\mathcal{S}_{i}+\mathcal{S}_{i}}{\mathcal{M}_{i,j}}$ (4)

In practice, we compute the Davies-Bouldin Index over minibatches of data and
average the result to get the final complexity measure. We choose a suitable
batch size based on the number of output classes such that there are a
reasonable number of samples from each class in a batch. We also find that
taking the mean as opposed to the maximum over clusters gives better results
on most tasks, possibly because this quantifies the separation of a class from
all other classes as opposed to only the closest cluster.

Figure 1: A visual depiction of how consistency of representations could aid
in generalization. Models in both Figure 1(a) and Figure 1(b) are saturated,
and thus correctly classify the training data $X$. However, the model in
Figure 1(b) has more consistent representations. Calculating the DB Index of
intermediate representations $A_{k}$ can then be predictive of generalization
gap.

We explore PCA and Max-Pooling for dimensionality reduction. Also, the Davies-
Bouldin Index could be computed for any intermediate layer. It makes sense
that middle and deeper layers of the network are most suitable for computing
the DB Index; these are likely to capture high level concepts in a compact
fashion. Initial layers are likely to detect edges and simple shapes and thus
have high variance, while the penultimate layer before the output may have
converged representations that are consistent irrespective of generalization
gap. We find that computing the DB Index on middle and deeper layers before
the penultimate layer gives highest scores for the public tasks. However, from
the perspective of comparing models with different topologies, we find that
the DB Index on the first layer representations is most suitable. We also
explore other measures of clustering quality, such as Silhouette Coefficient
and Calinski-Harabasz Index but find that the Davies Bouldin Index works best.

###  4.2 Robustness of Representations

Multiple studies have found the human visual system to be more robust to most
types of noise than deep neural networks [10, 11]. However, it has also been
shown that networks trained on perturbed images perform better than humans on
the exact perturbation types they were trained on, but generalize poorly on
other perturbation types [11]. Robustness of internal network representations
is thus closely related to generalization ability and a metric which
quantifies this could be predictive of generalization.

Supervised Deep Learning approaches are usually trained on the principle of
Empirical Risk Minimization, where the average error over the training samples
is minimized. Vicinal Risk Minimization, the formal name for data augmentation
and related techniques, aims to reduce the error on samples in the vicinity of
the training samples.

Measuring error on a vicinal dataset is one seemingly straightforward way to
predict performance of models on test data. This would allow us to determine
how robust a model is to nuisance or perturbations in features. However,
multiple problems arise here. One, data augmentation or perturbation
techniques are often dataset dependent. Creating a vicinal distribution
requires visualization and expert understanding of the underlying factors of
variation of the data, which is possible for natural image datasets but may
not be for the range and types of data deep neural networks are being used for
today. Secondly, such a vicinal distribution would only allow us to check
robustness in the input space. Thirdly, in a post-hoc diagnostic setting, it
is necessary to know whether networks are trained on a vicinal distribution,
and if they are, what sort of perturbations have they been trained to handle.
Checking performance of a model on a vicinal distribution may not tell us
anything about the test performance if our chosen perturbations take the data
into a neighbourhood that the model has never seen; all models are likely to
give out random results here, irrespective of their test performance. This is
especially likely in higher dimensional spaces, where it is easy to be close
to the training distribution and yet be in an unseen region. Moreover,
comparing models trained on vicinal distributions and those trained on the
original training dataset becomes difficult here.

Thus, while predicting generalization in a post-hoc setting, perturbations to
the training dataset or the intermediate representations have to be (i) data-
agnostic, and (ii) controlled, to ensure that they do not take representations
into a space that the model has never seen before.

Mixup [26] and manifold mixup [24] provide a potential solution that meets our
requirements. Mixup is data agnostic, and perturbs representations in
directions that our model has already seen, i.e. towards another training
sample.

However, mixup by itself is an ante-hoc procedure which trains models to
behave linearly between samples irrespective of class. If models have not been
trained on mixup, it is likely that models will perform poorly on mixed-up
samples irrespective of how well they generalize. We use a label-wise version
of mixup where we check performance on linear combinations of training samples
restricted to a class. Models are more likely to have linear behaviour within
samples of a particular class as opposed to samples of different classes,
especially if they haven’t been trained with mixup. We find that our variation
of mixup is especially useful in a post-hoc setting and achieves significantly
higher scores on all tasks of PGDL. Mathematically, we create a vicinal
distribution $\mathcal{V}$ made up of tuples $\\{\tilde{X},\tilde{y}\\}$, such
that

$\tilde{X}^{i}=\lambda X_{1}^{i}+(1-\lambda)X_{2}^{i}$ (5)
$\tilde{y}^{i}=y^{i}$ (6)

More generally, we create mixed-up representations $\tilde{A}$ such that,

$\tilde{A}^{i}=\lambda A_{1}^{i}+(1-\lambda)A_{2}^{i}$ (7)
$\tilde{y}^{i}=y^{i}$ (8)

Where inputs $X$ are denoted by $A_{0}$. Computing model performance on mixed-
up representations can then tell us about its generalization ability. The
mixup complexity measure based on representations at layer $k$ is then defined
as:

$\mathcal{C}_{k}=\sum_{i=1}^{\kappa}\frac{1}{N_{i}}\sum_{n=1}^{N_{i}}I\big{(}f_{w,k}(\tilde{A}_{n}^{i})[y_{i}]\leq\max_{j\neq\tilde{y}_{n}}f_{w,k}(\tilde{A}_{n}^{i})[j]\big{)}$
(9) $=\sum_{i=1}^{\kappa}\hat{L}\big{(}f_{w,k}(\tilde{A}^{i})\big{)}$ (10)

Where $f_{w,k}(\tilde{A})$ denotes the model mapping intermediate
representations to the output. We report scores for both mixup and manifold
mixup on various layers. However, we find that mixup on the input data is most
predictive of generalization. Mixup also has the added advantage of accounting
for random labelling. Networks are likely to misclassify linear combinations
of input samples where atleast one sample has been randomly labelled. This is
shown graphically in Figure 2. Thus, performance on mixed-up samples can also
be used as a proxy for the percentage of random labels the network has been
trained on.

Figure 2: A visualization of how mixup can help quantify random labelling.
$X^{i}$ depicts samples of class $i$, $X^{j(i)}$ depicts samples of class $j$
wrongly labelled as class $i$, $\tilde{X}$ depicts the mixed-up samples, and
$D$ shows the hypothetical margin (which overfits on the wrongly labelled
samples). The mixed-up points $\tilde{X}$ are likely to be misclassified when
either one or both samples in the linear combination are randomly labelled.

###  4.3 Separability of Representations

The entanglement of object manifolds is hypothesized to be one reason why
object recognition is difficult for humans. It has been proposed that the
human visual hierarchy untangles object manifolds to transform inseparable
manifolds into those which can be linearly or easily separated [7], behaviour
which has also been seen in deep neural networks [6].

Approximate formulations for the margin, the minimum distance to a decision
boundary between two classes, have been proposed before for deep neural
networks. Margin distributions provide a reasonable proxy for measuring
separability of representation manifolds in high dimensional spaces [9]. The
margin distribution has also been shown to be empirically predictive of
generalization in deep neural networks [12]. We use the margin distribution as
a measure of separability of representations in intermediate layers of the
network. Jiang et al. (2018) [12] show that the normalized margin distribution
becomes heavier tailed and shifts to the right as generalization gap
increases, and train a linear classifier to predict generalization gap from a
signature of the margin distribution. It may not always be possible to train
such a classifier in a real world setting. We replace the classifier simply by
a measure of central tendency of the margin distribution as a summary measure
and find that this gives good results. Moreover, we argue that the margin
distribution computed on a perturbed representation manifold or the vicinal
distribution is more useful as a predictor of generalization. Mathematically,
Let $D_{(i,j)}$ be the decision boundary between two classes $i$ and $j$, and
$X^{\prime}$ be a perturbed input sample. Let $A^{\prime}_{k}$ be the
corresponding intermediate layer representation. Then, the margin is defined
as:

$D_{(i,j),k}=\\{A^{\prime}_{k}\;|\;f_{w,k}(A^{\prime}_{k})[i]=f_{w,k}(A^{\prime}_{k})[j]\\}$
(11)

And the margin distance is approximated as:

$d_{f_{w,k},(i,j)}(A^{\prime}_{k})=\frac{f_{w,k}(A^{\prime}_{k})[i]-f_{w,k}(A^{\prime}_{k})[j]}{||\nabla_{A^{\prime}_{k}}f_{w,k}(A^{\prime}_{k})[i]-\nabla_{A^{\prime}_{k}}f_{w,k}(A^{\prime}_{k})[j]||_{2}}$
(12)

The margin complexity measure at layer $k$ is then:

$\mathcal{C}_{k}=-\theta(d_{f_{w,k},(i,j)}(A^{\prime}_{k}))$ (13)

Where $\theta$ is a distribution summary measure. We normalize the margin by
the total variation of the representation as in Jiang et al. (2018) [12]. The
intuition behind calculating the margin on a perturbed distribution is shown
in Figure 3. Figure 3(a) and Figure 3(b) show two models that correctly
classify samples $X$ from the original distribution as well as samples
$X^{\prime}$ from a perturbed distribution. The original data points $X$ are
at a similar distance from the margin for both models. However, the model in
Figure 3(b) also has a high margin for perturbed samples $X^{\prime}$. Thus,
it is expected that this model will better handle a wider range of data and
therefore generalizes better. Note that simply calculating the margin on
samples from the original distribution will not help here (since this will be
similar for both models), neither will trivially comparing classification
performance or model confidence on perturbed samples (since both models
correctly classify all points; we also show this empirically by comparing
Augment Margin and simple Augment Performance, shown in Table 1). We explore
simple data augmentation, mixup, and adverserial attacks for creating
perturbed distributions. We find that the margin distribution summary on
mixed-up samples is most predictive of generalization.

Figure 3: Calculating the margin on the perturbed samples tells us more about
generalization ability than the naive margin. The margin is shown in red, $X$
depicts samples from the original distribution, and $X^{\prime}$ are samples
from the vicinal distribution.

##  5 Experiments

We evaluate our methods on the 8 datasets provided as part of the PGDL
Competition. The public datasets of this competition contain 150 deep CNN
models trained on the CIFAR-10 and SVHN datasets. The CIFAR-10 models are VGG-
like, while the SVHN models follow a Network-in Network type of architecture.
These networks are trained with various hyperparameter choices - convolutional
width, number of layers, batch size, dropout, depth, and weight decay. All
models are trained till saturation, i.e. nearly perfect training accuracy,
while their test performance varies.

Complexity Measure |  | CIFAR-10 |  | SVHN
---|---|---|---|---
Davies Bouldin Index |  | 1.76 (25.22) |  | 7.98 (22.19)
Label-Wise Mixup |  | 0.00 |  | 22.75
Margin Summary |  | 1.87 |  | 39.95
Manifold Mixup |  | 2.24 |  | 12.11
DBI * Label-Wise Mixup |  | 0.00 |  | 32.05
Augment Margin Summary |  | 5.73 (15.66) |  | 44.60 (48.34)
Mixup Margin Summary |  | 1.11 |  | 47.33
prod-of-spec-over-margin |  | 1.54 |  | 5.13
prod-of-fro-over-margin |  | 2.84 |  | 4.92
Sharpness |  | 1.77 |  | 15.85
Augment Performance |  | 1.54 |  | 6.61
VC-Dimension |  | 0.00 |  | 0.00
Table 1: Mutual Information Scores for various representation based complexity measures as well as some baseline measures on CIFAR-10 and SVHN. Specific details for the measures as well as design choices for the task-specialized scores in brackets are given in the Appendix. Complexity Measure | CIFAR-10 | SVHN | PGDL Task 4 | PGDL Task 5
---|---|---|---|---
DBI * Label-Wise Mixup | 0.00 | 32.05 | 31.79 | 15.92
Augment Margin Summary | 5.73 | 44.60 | 47.22 | 22.82
Mixup Margin Summary | 1.11 | 47.33 | 43.22 | 34.57
Table 2: Mutual Information scores of our final complexity measures on public and development tasks of PGDL Complexity Measure | PGDL Task 6 | PGDL Task 7 | PGDL Task 8 | PGDL Task 9
---|---|---|---|---
DBI * Label-Wise Mixup | 43.99 | 12.59 | 9.24 | 25.86
Augment Margin Summary | 8.67 | 11.97 | 1.28 | 15.25
Mixup Margin Summary | 11.46 | 21.98 | 1.48 | 20.78
Table 3: Mutual Information scores of our final complexity measures on final
tasks of PGDL

Tables 1, 2 and 3 show the performance of our methods on the different tasks
of PGDL. We report scores of individual complexity measures as well as
combinations which achieve highest scores. Trivially combining measures is
difficult without access to the distribution of measures over the entire task;
scale differences between non-normalized measures across models makes it
difficult to create meaningful combinations. Our highest scoring complexity
measure on the final task combines two measures - DB Index and Label-Wise
Mixup, both of which have the advantage of being relative quantities and are
hence scale invariant, allowing us to combine them by simply taking a product.
We also report scores of some baseline measures taken from Jiang et al. (2019)
[13]. For some complexity measures, various design choices are possible and
lead to high scores on individual tasks. In such cases we show the scores of
generic/minimalist variants of each measure which we expect to run across
tasks irrespective of type of data and model architectures, as well as the
highest scores obtained by making intelligent design choices based on
knowledge of the data and model architectures (shown in brackets in Table 1),
which is achievable in a real-world setting. We find that calculating our
complexity measures over as little as 1% of the training samples is enough to
give reasonable estimates of generalization ability, underlining their
computational efficiency. Details about design choices are detailed in the
appendix. We show the scores of all complexity measures we tested on the
public tasks; however, due to limited access to the private tasks we only show
the scores of our final 3 solutions on these.

##  6 Discussion

A number of methods have been proposed before for understanding and predicting
generalization performance in deep learning. However, many of these are in the
form of upper bounds which are unsuitable to use in a practical setting, while
others are unreliable when dealing with a wide range of hyperparameters and
topological differences in models.

In this work, we have introduced complexity measures based on the quality of
representations of deep neural networks that are especially useful in a post-
hoc setting. We discuss our motivation and approach that led to a top score in
the NeurIPS Competition on Predicting Generalization in Deep Learning. We
demonstrate the effectiveness and ubiquity of our measures on a wide range of
tasks part of the PGDL Competition. We show that our measures can lead to
better mutual information scores than baseline measures such as spectral norm
and sharpness on CIFAR-10 and SVHN datasets. Our main observations and
conclusions are as follows:

  * •

Simple measures based on the quality of internal representation are predictive
of generalization gap. Specifically, we demonstrate how quantifying three
qualities of representations - consistency, robustness, and separability, can
predict generalization across a wide range of hyperparameter and model
architectures.

  * •

Our hypothesis aligns with neuroscientific theories of how the human visual
system creates invariant object representations in the face of controlled
naturalistic variations and has the ability to untangle inseparable
representation manifolds into easily separable ones.

  * •

We measure consistency of internal representations using the Davies Bouldin
Index. This indicates that clustered and compact representations of objects in
internal layers aid in generalization.

  * •

To measure robustness, we use a label-wise variation of Mixup. Mixup as a
technique is uniquely qualified to probe the invariance of representations
inside the model in a post-hoc setting since it is both data-agnostic and
controlled, i.e. it does not introduce external elements which might take the
model into unseen regions.

  * •

We measure separability of representations using the Margin Distribution,
which has been explored before. However, we find that the perturbed margin is
more expressive of generalization than just the margin.

While it is not surprising that better representations lead to better
generalization, there have been limited studies to explicitly understand the
relationship between these. We provide a way to think of generalization from
the perspective of the parallels of deep neural network representations with
the human visual system. We show that interpreted this way, existing
regularization techniques can be modified to be more suitable for predicting
generalization across a wide range of tasks in a post-hoc setting. We hope our
methods and analysis can provoke further research on understanding qualities
of representation and other intermediate measures that are responsible for
better generalization, as well as more critical evaluations of the parallels
between deep neural networks and the human visual system.

##  7 Acknowledgements

We’d like to thank the organizers of the NeurIPS 2020 Competition on
Predicting Generalization in Deep Learning for hosting this competition and
for providing a platform for us to test our hypotheses.

## References

  * Achille and Soatto [2018] Alessandro Achille and Stefano Soatto.  Emergence of invariance and disentanglement in deep representations.  _The Journal of Machine Learning Research_ , 19(1):1947–1980, 2018.
  * Arora et al. [2018] Sanjeev Arora, Rong Ge, Behnam Neyshabur, and Yi Zhang.  Stronger generalization bounds for deep nets via a compression approach.  _arXiv preprint arXiv:1802.05296_ , 2018.
  * Bartlett [1998] Peter L Bartlett.  The sample complexity of pattern classification with neural networks: the size of the weights is more important than the size of the network.  _IEEE transactions on Information Theory_ , 44(2):525–536, 1998.
  * Bengio et al. [2013] Yoshua Bengio, Aaron Courville, and Pascal Vincent.  Representation learning: A review and new perspectives.  _IEEE transactions on pattern analysis and machine intelligence_ , 35(8):1798–1828, 2013.
  * Cichy et al. [2016] Radoslaw Martin Cichy, Aditya Khosla, Dimitrios Pantazis, Antonio Torralba, and Aude Oliva.  Comparison of deep neural networks to spatio-temporal cortical dynamics of human visual object recognition reveals hierarchical correspondence.  _Scientific reports_ , 6:27755, 2016.
  * Cohen et al. [2020] Uri Cohen, SueYeon Chung, Daniel D Lee, and Haim Sompolinsky.  Separability and geometry of object manifolds in deep neural networks.  _Nature communications_ , 11(1):1–13, 2020.
  * DiCarlo and Cox [2007] James J DiCarlo and David D Cox.  Untangling invariant object recognition.  _Trends in cognitive sciences_ , 11(8):333–341, 2007.
  * Dziugaite and Roy [2017] Gintare Karolina Dziugaite and Daniel M Roy.  Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data.  _arXiv preprint arXiv:1703.11008_ , 2017.
  * Elsayed et al. [2018] Gamaleldin Elsayed, Dilip Krishnan, Hossein Mobahi, Kevin Regan, and Samy Bengio.  Large margin deep networks for classification.  In _Advances in neural information processing systems_ , pages 842–852, 2018.
  * Geirhos et al. [2017] Robert Geirhos, David HJ Janssen, Heiko H Schütt, Jonas Rauber, Matthias Bethge, and Felix A Wichmann.  Comparing deep neural networks against humans: object recognition when the signal gets weaker.  _arXiv preprint arXiv:1706.06969_ , 2017.
  * Geirhos et al. [2018] Robert Geirhos, Carlos RM Temme, Jonas Rauber, Heiko H Schütt, Matthias Bethge, and Felix A Wichmann.  Generalisation in humans and deep neural networks.  In _Advances in neural information processing systems_ , pages 7538–7550, 2018.
  * Jiang et al. [2018] Yiding Jiang, Dilip Krishnan, Hossein Mobahi, and Samy Bengio.  Predicting the generalization gap in deep networks with margin distributions.  _arXiv preprint arXiv:1810.00113_ , 2018.
  * Jiang et al. [2019] Yiding Jiang, Behnam Neyshabur, Hossein Mobahi, Dilip Krishnan, and Samy Bengio.  Fantastic generalization measures and where to find them.  In _International Conference on Learning Representations_ , 2019.
  * Karimi-Rouzbahani et al. [2017] Hamid Karimi-Rouzbahani, Nasour Bagheri, and Reza Ebrahimpour.  Invariant object recognition is a personalized selection of invariant features in humans, not simply explained by hierarchical feed-forward vision models.  _Scientific reports_ , 7(1):1–24, 2017.
  * Keskar et al. [2016] Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang.  On large-batch training for deep learning: Generalization gap and sharp minima.  _arXiv preprint arXiv:1609.04836_ , 2016.
  * Langford and Caruana [2001] John Langford and Rich Caruana.  (not) bounding the true error.  _Advances in Neural Information Processing Systems_ , 14:809–816, 2001.
  * Liao et al. [2016] Renjie Liao, Alex Schwing, Richard Zemel, and Raquel Urtasun.  Learning deep parsimonious representations.  In _Advances in Neural Information Processing Systems_ , pages 5076–5084, 2016.
  * Morcos et al. [2018] Ari S Morcos, David GT Barrett, Neil C Rabinowitz, and Matthew Botvinick.  On the importance of single directions for generalization.  _arXiv preprint arXiv:1803.06959_ , 2018.
  * Neyshabur et al. [2014] Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  In search of the real inductive bias: On the role of implicit regularization in deep learning.  _arXiv preprint arXiv:1412.6614_ , 2014.
  * Neyshabur et al. [2015] Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  Norm-based capacity control in neural networks.  In _Conference on Learning Theory_ , pages 1376–1401, 2015.
  * Nielsen et al. [2008] Kristina J Nielsen, Nikos K Logothetis, and Gregor Rainer.  Object features used by humans and monkeys to identify rotated shapes.  _Journal of Vision_ , 8(2):9–9, 2008.
  * Son et al. [2018] Sanghyun Son, Seungjun Nah, and Kyoung Mu Lee.  Clustering convolutional kernels to compress deep neural networks.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , pages 216–232, 2018.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * Verma et al. [2019] Vikas Verma, Alex Lamb, Christopher Beckham, Amir Najafi, Ioannis Mitliagkas, David Lopez-Paz, and Yoshua Bengio.  Manifold mixup: Better representations by interpolating hidden states.  In _International Conference on Machine Learning_ , pages 6438–6447. PMLR, 2019.
  * Zhang et al. [2016] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  _arXiv preprint arXiv:1611.03530_ , 2016.
  * Zhang et al. [2017] Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization.  _arXiv preprint arXiv:1710.09412_ , 2017.

##  8 Appendix

###  8.1 Details for Table 1

Here, we detail the various design choices for complexity measures in Table 1.

  1. 1.

For the generic version of the Davies Bouldin Index, we use max-pooling with a
filter size of 4 and stride of 1 for spatial aggregation/dimensionality
reduction for better clustering. We compute the DB Index on the first layer,
since this is most suitable when model architectures are unknown. For the
task-specialized version for CIFAR-10 and SVHN, we compute DB Index on the
third-from-last layer. This may be a max-pool layer or a convolutional layer
for different models; we find that comparing across these layer types still
gives reasonable mutual information scores.

  2. 2.

For Label-Wise Mixup, we hold $\lambda$ constant at 0.5.

  3. 3.

For the Margin Summary, we use the mean of the quantile signature in Jiang et
al. (2018) [12]. We also explore simply the mean of the margin distribution
and find that this gives similar results.

  4. 4.

Manifold Mixup scores are reported on the first convolutional layer since we
are comparing models with varying topologies; we do not conduct a study to
explore the suitability of different intermediate layers for this measure.

  5. 5.

For the Augment Margin Summary, we use the following augmentations: Random Hue
with a maximum value of 0.5, random saturation with a saturation value sampled
from [0.6, 1.2], random brightness with a maximum value of 0.5, random
contrast sampled from [0.7, 1.0], random clipping/zoom [1 to 15 %], and left-
right flips. For the generic version, we only use contrast, left-right flips,
and zoom. Random seeds are fixed so that all models are tested on the same
vicinal distribution.

  6. 6.

Both norm-based measures, prod-of-spec-over margin and prod-of-fro-over-
margin, are computed as in Jiang et al. (2019) [13].

  7. 7.

For Augment Performance, we simply check the categorical cross-entropy of the
model on augmented samples. Augmentations are same as those for the Augmented
Margin.