# Deep Parametric Continuous Convolutional Neural Networks

Shenlong Wang${}^{1,3,^{\ast}}$ Simon Suo${}^{2,3,^{\ast}}$ Wei-Chiu Ma3
Andrei Pokrovsky3 Raquel Urtasun1,3
1University of Toronto, 2University of Waterloo, 3Uber Advanced Technologies
Group
{slwang, suo, weichiu, andrei<EMAIL_ADDRESS>
###### Abstract

Standard convolutional neural networks assume a grid structured input is
available and exploit discrete convolutions as their fundamental building
blocks. This limits their applicability to many real-world applications. In
this paper we propose Parametric Continuous Convolution, a new learnable
operator that operates over non-grid structured data. The key idea is to
exploit parameterized kernel functions that span the full continuous vector
space. This generalization allows us to learn over arbitrary data structures
as long as their support relationship is computable. Our experiments show
significant improvement over the state-of-the-art in point cloud segmentation
of indoor and outdoor scenes, and lidar motion estimation of driving scenes.

##  1 Introduction

Discrete convolutions are the most fundamental building block of modern deep
learning architectures. Its efficiency and effectiveness relies on the fact
that the data appears naturally in a dense grid structure (e.g., 2D grid for
images, 3D grid for videos). However, many real world applications such as
visual perception from 3D point clouds, mesh registration and non-rigid shape
correspondences rely on making statistical predictions from non-grid
structured data. Unfortunately, standard convolutional operators cannot be
directly applied in these cases.

Multiple approaches have been proposed to handle non-grid structured data. The
simplest approach is to voxelize the space to form a grid where standard
discrete convolutions can be performed [30, 24]. However, most of the volume
is typically empty, and thus this results in both memory inefficiency and
wasted computation. Geometric deep learning [3, 15] and graph neural network
approaches [25, 16] exploit the graph structure of the data and model the
relationship between nodes. Information is then propagated through the graph
edges. However, they either have difficulties generalizing well or require
strong feature representations as input to perform competitively. End-to-end
learning is typically performed via back-propagation through time, but it is
difficult to learn very deep networks due to the memory limitations of modern
GPUs.

In contrast to the aforementioned approaches, in this paper we propose a new
learnable operator, which we call parametric continuous convolution. The key
idea is a parameterized kernel function that spans the full continuous vector
space. In this way, it can handle arbitrary data structures as long as its
support relationship is computable. This is a natural extension since objects
in the real-world such as point clouds captured from 3D sensors are
distributed unevenly in continuous domain. Based upon this we build a new
family of deep neural networks that can be applied on generic non-grid
structured data. The proposed networks are both expressive and memory
efficient.

We demonstrate the effectiveness of our approach in both semantic labeling and
motion estimation of point clouds. Most importantly, we show that very deep
networks can be learned over raw point clouds in an end-to-end manner. Our
experiments show that the proposed approach outperforms the state-of-the-art
by a large margin in both outdoor and indoor 3D point cloud segmentation
tasks, as well as lidar motion estimation in driving scenes. Importantly, our
outdoor semantic labeling and lidar flow experiments are conducted on a very
large scale dataset, containing 223 billion points captured by a 3D sensor
mounted on the roof of a self-driving car. To our knowledge, this is 2 orders
of magnitude larger than any existing benchmark.

##  2 Related Work

#### Deep Learning for 3D Geometry:

Deep learning approaches that exploit 3D geometric data have recently become
populer in the computer vision community. Early approaches convert the 3D data
into a two-dimensional RGB + depth image [17, 10] and exploit conventional
convolutional neural networks (CNNs). Unfortunately, this representation does
not capture the true geometric relationships between 3D points (i.e.
neighboring pixels could be potentially far away geometrically). Another
popular approach is to conduct 3D convolutions over volumetric representations
[30, 21, 24, 9, 18]. Voxelization is employed to convert point clouds into a
3D grid that encodes the geometric information. These approaches have been
popular in medical imaging and indoor scene understanding, where the volume is
relatively small. However, typical voxelization approaches sacrifice precision
and the 3D volumetric representation is not memory efficient. Sparse
convolutions [9] and advanced data structures such as oct-trees [24] have been
used to overcome these difficulties. Learning directly over point clouds has
only been studied very recently. The pioneer work of PointNet [20], learns an
MLP over individual points and aggregates global information using pooling.
PointNet++ [22], the follow-up, improves the ability to capture local
structures through a multi-scale grouping strategy.

#### Graph Neural Networks:

Graph neural networks (GNNs) [25] are generalizations of neural networks to
graph structured data. Early approaches apply neural networks either over the
hidden representation of each node or the messages passed between adjacent
nodes in the graph, and use back-propagation through time to conduct learning.
Gated graph neural networks (GGNNs) [16] exploit gated recurrent units along
with modern optimization techniques, resulting in improved performance. In
[23], GGNNs are applied to point cloud segmentation, achieving significant
improvements over the state-of-the-art. One of the major difficulties of graph
neural networks is that propagation is conducted in a synchronous manner and
thus it is hard to scale up to graphs with millions of nodes. Inference in
graphical models as well as recurrent neural networks can be seen as special
cases of graph neural networks.

#### Graph Convolution Networks:

An alternative formulation is to learn convolution operations over graphs.
These methods can be categorized into spectral and spatial approaches
depending on which domain the convolutions are applied to. For spectral
methods, convolutions are converted to multiplication by computing the graph
Laplacian in Fourier domain [4, 2, 31]. Parameterized spectral filters can be
incorporated to reduce overfitting [4]. These methods are not feasible for
large scale data due to the expensive computation, since there is no FFT-like
trick over generic graph. Spatial approaches directly propagate information
along the node neighborhoods in the graph. This can be implemented either
through low-order approximation of spectral filtering[6, 15, 7], or diffusion
in a support domain [19, 2, 27, 31, 26]. Our approach generalizes spatial
approaches in two ways: first, we use more expressive convolutional kernel
functions; second, the output of the convolution could be any point in the
whole continuous domain.

Figure 1: Unlike grid convolution, parametric continuous convolution uses
kernel functions that are defined for arbitrary points in the continuous
support domain. As a result, it is possible to output features at points not
seen in the input.

#### Other Approaches:

Edge-conditioned filter networks [27] use a weighting network to communicate
between adjacent nodes on the graph [13] conditioned on edge labels, which is
primarily formulated as relative point locations. In contrast, our approach is
not constrained to a fixed graph structure, and has the flexibility to output
features at arbitrary points over the continuous domain. In a concurrent work,
[26] uses similar parametric function form $f(\mathbf{x}_{i}-\mathbf{x}_{j})$
to aggregate information between points. However, they only use shallow
isotropic gaussian kernels to represent the weights, while we use expressive
deep networks to parameterize the continuous filters.

Figure 2: Detailed Computation Block for the Parametric Continuous Convolution
Layer.

##  3 Deep Parametric Continuous CNNs

###  3.1 Parametric Continuous Convolutions

Standard CNNs use discrete convolutions (i.e., convolutions defined over
discrete domain) as basic operations.

$h[n]=(f\ast g)[n]=\sum_{m=-M}^{M}f[n-m]g[m]$

where $f:\mathcal{G}\rightarrow\mathbb{R}$ and
$g:\mathcal{S}\rightarrow\mathbb{R}$ are functions defined over the support
domain of finite integer set: $\mathcal{G}=\mathcal{Z}^{D}$ and
$\mathcal{S}=\\{-M,-M+1,...,M-1,M\\}^{D}$ respectively.

In contrast, continuous convolutions can be defined as

$h(\mathbf{x})=(f\ast
g)(\mathbf{x})=\int_{-\infty}^{\infty}f(\mathbf{y})g(\mathbf{x}-\mathbf{y})d\mathbf{y}$
(1)

where both the kernel $g:\mathcal{S}\rightarrow\mathbb{R}$ and the feature
$f:\mathcal{G}\rightarrow\mathbb{R}$ are defined as continuous functions over
the support domain $\mathcal{G}=\mathbb{R}^{D}$ and
$\mathcal{S}=\mathbb{R}^{D}$ respectively.

Continuous convolutions require the integration in Eq. (1) to be analytically
tractable. Unfortunately, this is not possible for real-world applications,
where the input features are complicated and non-parametric, and the
observations are sparse points sampled over the continuous domain.

Motivated by monte-carlo integration [5] we derive our continuous convolution
operator. In particular, given continuous functions $f$ and $g$ with a finite
number of input points ${\mathbf{y}_{i}}$ sampled from the domain, the
convolution at an arbitrary point $\mathbf{x}$ can be approximated as:

$h(\mathbf{x})=\int_{-\infty}^{\infty}f(\mathbf{y})g(\mathbf{x}-\mathbf{y})d\mathbf{y}\approx\sum_{i}^{N}\frac{1}{N}f(\mathbf{y}_{i})g(\mathbf{x}-\mathbf{y}_{i})$

The next challenge we need to solve is constructing the continuous
convolutional kernel function $g$. Conventional 2D and 3D discrete convolution
kernels are parameterized in a way that each point in the support domain is
assigned a value (_i.e_. the kernel weight). Such a parameterization is
infeasible for continuous convolutions, since the kernel function $g$ is
defined over an infinite number of points (i.e., has infinite support).
Instead, in this paper we propose to use parametric continuous functions to
model $g$. We name our approach Parametric Continuous Convolutions. In
particular, we use a multi-layer perceptron (MLP) as the approximator. With
reference to the universal approximation theorem of [12], MLPs are expressive
and capable of approximating continuous functions over $\mathbb{R}^{n}$. Thus
we define:

$g(\mathbf{z};\theta)=MLP(\mathbf{z};\theta)$

The kernel function $g(\mathbf{z};\theta):\mathbb{R}^{D}\rightarrow\mathbb{R}$
spans the full continuous support domain while remaining parametrizable by a
finite number of parameters. Note that other choices such as polynomials are
possible, however low-order polynomials are not expressive, whereas learning
high-order polynomials can be numerically unstable for back-propagation.

Figure 3: Architecture of the Deep Parametric Continuous CNNs for Semantic
Labeling Task.

###  3.2 From Convolutions to Deep Networks

In this section, we first design a new convolution layer based on the
parametric continuous convolutions derived in the previous subsection. We then
propose a deep learning architecture using this new convolution layer.

#### Parametric Continuous Convolution Layer:

Note that, unlike standard discrete convolutions which are conducted over the
same point set, the input and output points of our parametric continuous
convolution layer can be different. This is important for many practical
applications, where we want to make dense predictions based on partial
observations. Furthermore, this allow us to abstract information from
redundant input points (i.e., pooling). As a consequence, the input of each
convolution layer contains three parts: the input feature vector
$\mathcal{F}=\\{\mathbf{f}_{\mathrm{in},j}\in\mathbb{R}^{F}\\}$, the
associated locations in the support domain $\mathcal{S}=\\{\mathbf{y}_{j}\\}$,
as well as the output domain locations $\mathcal{O}=\\{\mathbf{x}_{i}\\}$. For
each layer, we first evaluate the kernel function
$g_{d,k}(\mathbf{y}_{i}-\mathbf{x}_{j};\theta)$ for all
$\mathbf{x}_{j}\in\mathcal{S}$ and all $\mathbf{y}_{i}\in\mathcal{O}$, given
the parameters $\theta$. Each element of the output feature vector is then
computed as:

${h}_{k,i}=\sum_{d}^{F}\sum_{j}^{N}g_{d,k}(\mathbf{y}_{i}-\mathbf{x}_{j})f_{d,j}$

Let $N$ be the number of input points, $M$ be the number of output points, and
$D$ the dimensionality of the support domain. Let $F$ and $O$ be predefined
input and output feature dimensions respectively. Note that these are
hyperparameters of the continuous convolution layer analogous to input and
output feature dimensions in standard grid convolution layers. Fig. 1 depicts
our parametric continuous convolutions in comparison with conventional grid
convolution. Two major differences are highlighted: 1) the kernel function is
continuous given the relative location in support domain; 2) the input/ouput
points could be any points in the continuous domain as well and can be
different.

|  |  |  |  |
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
Input | Ground Truth | Ours PCCN | Input | Ground Truth | Ours PCCN
Figure 4: Semenatic Segmentation Results on Stanford Indoor3D Dataset

#### Deep Parametric Continuous CNNs:

Using the parametric continuous convolution layers as building blocks, we can
construct a new family of deep networks which operates on unstructured data
defined in a topological group under addition. In the following discussions,
we will focus on multi-diumensional euclidean space, and note that this is a
special case. The network takes the input features and their associated
positions in the support domain as input. Then the hidden representations are
generated from successive parametric continuous convolution layers. Following
standard CNN architectures, we can add batch normalization, non-linearities
and residual connections between layers. Pooling can also be employed over the
support domain to aggregate information. In practice, we find adding residual
connection between parametric continuous convolution layers is critical to
help convergence. Please refer to Fig. 2 for an example of the computation
graph of a single layer, and to Fig. 3 for an example of the network
architecture employed for our indoor semantic segmentation task.

#### Learning:

All of our building blocks are differentiable, thus our networks can be
learned through back-prop:

$\frac{\partial h}{\partial\theta}=\frac{\partial h}{\partial
g}\cdot\frac{\partial
g}{\partial\theta}=\sum_{d}^{F}\sum_{j}^{N}f_{d,j}\cdot\frac{\partial
g}{\partial\theta}$

|  |  |
---|---|---|---
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
Ground Truth | 3D-FCN | Ours PCCN | Ours 3D-FCN+PCCN
Figure 5: Semenatic Segmentation Results on Driving Scene Dataset; Colored:
correct prediciton; white: wrong prediciton.

###  3.3 Discussions

#### Locality Enforcing Continuous Convolution:

Standard grid convolution are computed over a limited kernel size $M$ to keep
locality. Similarly, locality can be enforced in our parametric continuous
convolutions by constraining the influence of the function $g$ to points close
to $\mathbf{x}$, _i.e_.,

$g(\mathbf{z})=MLP(\mathbf{z})w(\mathbf{z})$

where $w(\cdot)$ is a modulating window function. This can be achieved in
differently. First, we can constrain the cardinality of its local support
domain and only keep non-zero kernel values for its K-nearest neighbors:
$w(\mathbf{z})=\mathbf{1}_{\mathbf{z}\in\mathrm{KNN}(\mathcal{S},\mathbf{x})}$.
Alternatively we can keep non-zero kernel values for points within a fixed
radius $r$: $w(\mathbf{z})=\mathbf{1}_{||\mathbf{z}||_{2}<r}$.

#### Efficient Continuous Convolution:

For each continuous convolution layer, the kernel function is evaluated
$N\times|\mathcal{S}|\times F\times O$ times, where $|\mathcal{S}|$ is the
cardinality of the support domain, and the intermediate weight tensor is
stored for backpropagation. This is expensive in practice, especially when
both the number of points and the feature dimension are large. With the
locality enforcing formulation, we can constrain the cardinality of
$\mathcal{S}$. Furthermore, motivated by the idea of separable filters, we use
the fact that this computation can be factorized if the kernel function value
across different output dimensionality is shared. That is to say, we can
decompose the weight tensor $W\in\mathbb{R}^{N\times|\mathcal{S}|\times
F\times O}$ into two tensors $W_{1}=\mathbb{R}^{F\times O}$ and
$W_{2}=\mathbb{R}^{N\times|\mathcal{S}|\times O}$ , where $W_{1}$ is a linear
weight matrix and $W_{2}$ is evaluated through the MLP. With this
optimization, only $N\times|\mathcal{S}|\times O$ kernel evaluations need to
be computed and stored. Lastly, in inference stage, through merging the
operations of batchnorm and fc layer in MLP, 3x speed boosting can be
achieved.

#### Special Cases:

Many previous convolutional layers are special cases of our approach. For
instance, if the points are sampled over the finite 2D grid we recover
conventional 2D convolutions. If the support domain is defined as
concatenation of the spatial vector and feature vector with a gaussian kernel
$g(\cdot)$, we recover the bilateral filter. If the support domain is defined
as the neighboring vertices of a node we recover the first-order spatial graph
convolution [15].

|  |
---|---|---
Figure 6: Semantic Labeling on KITTI Dataset without Retraining

##  4 Experimental Evaluation

We demonstrate the effectiveness of our approach in the tasks of semantic
labeling and motion estimation of 3D point clouds, and show state-of-the-art
performance. We conduct point-wise semantic labeling experiments over two
datasets: a very large-scale outdoor lidar semantic segmentation dataset that
we collected and labeled in house and a large indoor semantic labeling
dataset. To our knowledge, these are the largest real-world outdoor and indoor
datasets that are available for this task. The datasets are fully labeled and
contain 137 billion and 629 million points respectively. The lidar flow
experiment is also conducted on this dataset with ground-truth 3D motion label
for each point.

Method |  blue!25 mIOU |  blue!25 mAcc | ceiling | floor | wall | beam | column | window | door | chair | table | bookcase | sofa | board | clutter
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
PointNet [20] |  blue!2541.09 |  blue!2548.98 | 88.80 | 97.33 | 69.80 | 0.05 | 3.92 | 46.26 | 10.76 | 52.61 | 58.93 | 40.28 | 5.85 | 26.38 | 33.22
3D-FCN-TI [29] |  blue!2547.46 |  blue!25 54.91 | 90.17 | 96.48 | 70.16 | 0.00 | 11.40 | 33.36 | 21.12 | 76.12 | 70.07 | 57.89 | 37.46 | 11.16 | 41.61
SEGCloud [29] |  blue!25 48.92 |  blue!2557.35 | 90.06 | 96.05 | 69.86 | 0.00 | 18.37 | 38.35 | 23.12 | 75.89 | 70.40 | 58.42 | 40.88 | 12.96 | 41.60
Ours PCCN |  blue!25 58.27 |  blue!25 67.01 | 92.26 | 96.20 | 75.89 | 0.27 | 5.98 | 69.49 | 63.45 | 66.87 | 65.63 | 47.28 | 68.91 | 59.10 | 46.22
Table 1: Semantic Segmentation Results on Stanford Large-Scale 3D Indoor Scene
Dataset

###  4.1 Semantic Segmentation of Indoor Scenes

#### Dataset:

We use the Stanford large-scale 3D indoor scene dataset [1] and follow the
training and testing procedure used in [29]. We report the same metrics, i.e.,
mean-IOU, mean class accuracy (TP / (TP + FN)) and class-wise IOU. The input
is six dimensional and is composed of the xyz coordinates and RGB color
intensity. Each point is labeled with one of 13 classes shown in Tab. 1.

#### Competing Algorithms:

We compare our approach to PointNet [20] and SegCloud [29]. We evaluate the
proposed end-to-end continuous convnet with eight continuous convolution
layers (Ours PCCN). The kernels are defined over the continuous support domain
of 3D Euclidean space. Each intermediate layer except the last has 32
dimensional hidden features followed by batchnorm and ReLU nonlinearity. The
dimension of the last layer is 128. We observe that the distribution of
semantic labels within a room is highly correlated with the room type (_e.g_.
office, hallway, conference room, _etc_.). Motivated by this, we apply max
pooling over all the points in the last layer to obtain a global feature,
which is then concatenated to the output feature of each points in the last
layer, resulting in a 256 dimensional feature. A fully connected layer with
softmax activation is used to produce the final logits. Our network is trained
end-to-end with cross entropy loss, using Adam optimizer.

#### Results:

As shown in Tab. 1 our approach outperforms the state-of-the-art by 9.3% mIOU
and 9.6% mACC. Fig. 4 shows qualitative results. Despite the diversity of
geometric structures, our approach works very well. Confusion mainly occurs
between columns vs walls and window vs bookcase. It is also worth noting that
our approach captures visual information encoded in RGB channels. The last row
shows two failure cases. In the first one, the door in the washroom is labeled
as clutter whearas our algorithm thinks is door. In the second one, the board
on the right has a window-like texture, which makes the algorithm predict the
wrong label.

|
---|---
Flow Field | Overlay of Target and Warped Source
Figure 7: Right: purple shows target frame, yellow shows source frame warped
to target frame using ground truth flow

###  4.2 Semantic Segmentation of Driving Scenes

#### Dataset:

We first conduct experiments on the task of point cloud segmentation in the
context of autonomous driving. Each point cloud is produced by a full sweep of
a roof-mounted Velodyne-64 lidar sensor driving in several cities in North
America. The dataset is composed of snippets each having 300 consecutive
frames. The training and validation set contains 11,337 snippets in total
while the test set contains 1,644 snippets. We report metrics on a subset of
the test set which is generated by sampling 10 frames from each snippet to
avoid bias brought due to scenes where the ego-car is static (e.g., when
waiting at a traffic light). Each point is labeled with one of seven classes
defined in Tab. 2. We adopt mean intersection-over-union (meanIOU) and point-
wise accuracy (pointAcc) as our evaluation metrics.

Method |  blue!25 pACC |  blue!25 mIOU | vehicle | bicyclist | pedestrian | motorcycle | animal | background | road | params size
---|---|---|---|---|---|---|---|---|---|---
PointNet [20] |  blue!25 91.96 |  blue!25 38.05 | 76.73 | 2.85 | 6.62 | 8.02 | 0.0 | 89.83 | 91.96 | 20.34MB
3D-FCN [11] |  blue!25 94.31 |  blue!25 49.28 | 86.74 | 22.30 | 38.26 | 17.22 | 0.98 | 86.91 | 92.56 | 74.66MB
Ours PCCN |  blue!25 94.56 |  blue!25 46.35 | 86.62 | 8.31 | 41.84 | 7.24 | 0.00 | 87.27 | 93.20 | 9.34MB
Ours 3D-FCN+PCCN |  blue!25 95.45 |  blue!25 58.06 | 91.83 | 40.23 | 47.74 | 42.91 | 1.25 | 89.27 | 93.18 | 74.67MB
Table 2: Semenatic Segmentation Results on Driving Scenes Dataset

#### Baselines:

We compare our approach to the point cloud segmentation network (PointNet)
[20] and a 3D fully convolutional network (3D-FCN) conducted over a 3D
occupancy grid. We use a resolution of 0.2m for each voxel over a
160mx80mx6.4m range. This results in an occupancy grid encoded as a tensor of
size 800x400x32. We define a voxel to be occupied if it contains at least one
point. We use ResNet-50 as the backbone and replace the last average pooling
and fully connected layer with two fully convolutional layers and a trilinear
upsampling layer to obtain dense voxel predictions. The model is trained from
scratch with the Adam optimizer[14] to minimize the class-reweighted cross-
entropy loss. Finally, the voxel-wise predictions are mapped back to the
original points and metrics are computed over points. We adapted the open-
sourced PointNet model onto our dataset and trained from scratch. The
architecture and loss function remain the same with the original paper, except
that we removed the point rotation layer since it negatively impacts
validation performance on this dataset.

|  |  |
---|---|---|---
|  |  |
Ground Truth | Ours 3D-FCN+PCCN | Ground Truth | Ours 3D-FCN+PCCN
Figure 8: Lidar Flow Results on Driving Scene Dataset

#### Our Approaches:

We evaluate two versions of our approach. Our first instance conducts
continuous convolutions directly over the raw xyz-intensity lidar points (Ours
PCCN). Our second version (Ours 3D-FCN+PCCN) performs continuous convolutions
over the features extracted from 3D-FCN. Ours PCCN has 16 continuous conv
layers with residual connections, batchnorm and ReLU non-linearities. We use
the spatial support in $\mathbb{R}^{3}$ to define our kernel. We train the
network with point-wise cross-entropy loss and Adam [14] optimizer. In
contrast, Ours 3D-FCN+PCCN model has 7 residual continuous convolutional
layers on top of the trained 3D-FCN model and performs end-to-end fine-tuning
using Adam optimizer.

#### Results:

As shown in Tab. 2, by exploiting sophisticated feature via 3D convolutions,
3D-FCN+PCCN results in the best performance. Fig. 5 shows qualitative
comparison between models. As shown in the figure, all models produce good
results. Performance differences often result from ambiguous regions. In
particular, we can see that the 3D-FCN model oversegements the scene: it
mislabels a background pole as vehicle (red above egocar), nearby spurirous
points as bicyclist (green above egocar), and a wall as pedestrian (purple
near left edge). This is reflected in the confidence map (as bright regions).
We observe a significant improvement in our 3D-CNN + PCCN model, with all of
the above corrected with high confidence. For more results and videos please
refer to the supplementary material.

#### Model Sizes:

We also compare the model sizes of the competing algorithms in Tab. 2. In
comparison to the 3D-FCN approach, the end-to-end continuous convolution
network’s model size is eight times smaller , while achieving comparable
results. And the 3D-FCN+PCCN is just 0.01MB larger than 3D-FCN, but the
performance is improved by a large margin in terms of mean IOU.

#### Complexity and Runtime

We benchmark the proposed model’s runtime over a GTX 1080 Ti GPU and Xeon
E5-2687W CPU with 32 GB Memory. The forward pass of a 8-layer PCCN model (32
feature dim in each layer with 50 neighbours) takes 33ms. The KD-Tree
neighbour search takes 28 ms. The end-to-end computation takes 61ms. The
number of operations of each layer is 1.32GFLOPs.

#### Generalization:

To demonstrate the generalization ability of our approach, we evaluate our
model, trained with only North American scenes, on the KITTI dataset [8],
which was captured in Europe. As shown in Fig. 6, the model achieves good
results, with well segmented dynamic objects, such as vehicles and
pedestrians.

###  4.3 Lidar Flow

#### Dataset:

We also validate our proposed method over the task of lidar based motion
estimation, refered to as lidar flow. In this task, the input is two
consecutive frames of lidar sweep. The goal is to estimation the 3D motion
field for each point in the first frame, to undo both ego-motion and the
motion of dynamic objects. The ground-truth ego-motion is computed through a
comprehensive filters that take GPS, IMU as well as ICP based lidar alignment
against pre-scaned 3D geometry of the scene as input. And the ground-truth
6DOF dynamics object motion is estimated from the temporal coherent 3D object
tracklet, labeled by in-house annotators. Combining both we are able to get
the ground-truth motion field. Fig. 7 shows the colormapped flow field and the
overlay between two frames after undoing per-point motion. This task is
crucial for many applications, such as multi-rigid transform alignment, object
tracking, global pose estimation, _etc_. The training and validation set
contains 11,337 snippets while the test set contains 1,644 snippets. We use
110k frame pairs for training and validation, and 16440 frame pairs for
testing. End-point error, and outlier percentage at 10 cm and 20 cm are used
as metric.

#### Competing Algorithms:

We compare against the 3D-FCN baseline using the same architecture and
volumetric representation as used in Sec. 4.2. We also adopt a similar 3D-FCN
+ PCCN architecture with 7 residual continuous convolution layers added as a
polishing network. In this task, we remove the ReLU nonlinearity and supervise
the PCCN layers with MSE loss at every layer. The training objective function
is mean square error loss between the ground-truth flow vector and the
prediction.

#### Results:

Tab. 3 reports the quantitative results. As shown in the table, our
3D-FCN+PCCN model outperforms the 3D-FCN by 0.351cm in end-point error and our
method reduces approximately $20\%$ of the outliers. Fig. 18 shows sample flow
predictions compared with ground truth labels. As shown in the figure, our
algorithm is able to capture both global motion of the ego-car including self
rotation, and the motion of each dynamic objects in the scene. For more
results please refer to our supplementary material.

Method | EPE (cm) |  Outlier$\%_{10}$ |  Outlier$\%_{20}$
---|---|---|---
3D-FCN | 8.161 | 25.92% | 7.12 %
Ours 3D-FCN+PCCN | 7.810 | 19.84% | 5.97%
Table 3: Lidar Flow Results on Driving Scenes Dataset

##  5 Conclusions

We have presented a new learnable convolution layer that operates over non-
grid structured data. Our convolution kernel function is parameterized by
multi-layer perceptrons and spans the full continuous domain. This allows us
to design a new deep learning architecture that can be applied to arbitrary
structured data, as long as the support relationships between elements are
computable. We validate the performance on point cloud segmentation and motion
estimation tasks, over very large-scale datasets with up to 200 bilion points.
The proposed network achieves state-of-the-art performance on all the tasks
and datasets.

## References

  * [1] I. Armeni, O. Sener, A. R. Zamir, H. Jiang, I. Brilakis, M. Fischer, and S. Savarese.  3d semantic parsing of large-scale indoor spaces.  In CVPR, 2016.
  * [2] D. Boscaini, J. Masci, E. Rodolà, and M. Bronstein.  Learning shape correspondence with anisotropic convolutional neural networks.  In NIPS, 2016.
  * [3] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst.  Geometric deep learning: going beyond euclidean data.  IEEE Signal Processing Magazine, 2017.
  * [4] J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun.  Spectral networks and locally connected networks on graphs.  ICLR, 2014.
  * [5] R. E. Caflisch.  Monte carlo and quasi-monte carlo methods.  Acta numerica, 1998.
  * [6] M. Defferrard, X. Bresson, and P. Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  In NIPS, 2016.
  * [7] D. K. Duvenaud, D. Maclaurin, J. Iparraguirre, R. Bombarell, T. Hirzel, A. Aspuru-Guzik, and R. P. Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In NIPS, 2015.
  * [8] A. Geiger, P. Lenz, and R. Urtasun.  Are we ready for autonomous driving? the kitti vision benchmark suite.  In CVPR, 2012.
  * [9] B. Graham and L. van der Maaten.  Submanifold sparse convolutional networks.  arXiv, 2017.
  * [10] S. Gupta, R. Girshick, P. Arbeláez, and J. Malik.  Learning rich features from rgb-d images for object detection and segmentation.  In ECCV, 2014.
  * [11] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In CVPR, 2016.
  * [12] K. Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural networks, 1991.
  * [13] X. Jia, B. De Brabandere, T. Tuytelaars, and L. V. Gool.  Dynamic filter networks.  In NIPS. 2016.
  * [14] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  ICLR, 2015.
  * [15] T. N. Kipf and M. Welling.  Semi-supervised classification with graph convolutional networks.  arXiv, 2016.
  * [16] Y. Li, D. Tarlow, M. Brockschmidt, and R. Zemel.  Gated graph sequence neural networks.  arXiv, 2015.
  * [17] J. Long, E. Shelhamer, and T. Darrell.  Fully convolutional networks for semantic segmentation.  In CVPR, 2015.
  * [18] D. Maturana and S. Scherer.  Voxnet: A 3d convolutional neural network for real-time object recognition.  In IROS, 2015.
  * [19] F. Monti, D. Boscaini, J. Masci, E. Rodolà, J. Svoboda, and M. M. Bronstein.  Geometric deep learning on graphs and manifolds using mixture model cnns.  CVPR, 2017.
  * [20] C. R. Qi, H. Su, K. Mo, and L. J. Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  2016\.
  * [21] C. R. Qi, H. Su, M. Nießner, A. Dai, M. Yan, and L. J. Guibas.  Volumetric and multi-view cnns for object classification on 3d data.  In CVPR, 2016.
  * [22] C. R. Qi, L. Yi, H. Su, and L. J. Guibas.  Pointnet++: Deep hierarchical feature learning on point sets in a metric space.  2017\.
  * [23] X. Qi, R. Liao, J. Jia, S. Fidler, and R. Urtasun.  3d graph neural networks for rgbd semantic segmentation.  In CVPR, 2017.
  * [24] G. Riegler, A. O. Ulusoys, and A. Geiger.  Octnet: Learning deep 3d representations at high resolutions.  2017\.
  * [25] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini.  The graph neural network model.  TNN, 2009.
  * [26] K. T. Schütt, P. Kindermans, H. Sauceda, S. Chmiela, A. Tkatchenko, and K. Müller.  Schnet: A continuous-filter convolutional neural network for modeling quantum interactions.  arXiv, 2017.
  * [27] M. Simonovsky and N. Komodakis.  Dynamic edge-conditioned filters in convolutional neural networks on graphs.  CVPR, 2017.
  * [28] H. Su, S. Maji, E. Kalogerakis, and E. Learned-Miller.  Multi-view convolutional neural networks for 3d shape recognition.  In ICCV, 2015.
  * [29] L. P. Tchapmi, C. B. Choy, I. Armeni, J. Gwak, and S. Savarese.  Segcloud: Semantic segmentation of 3d point clouds.  arXiv preprint arXiv:1710.07563, 2017.
  * [30] Z. Wu, S. Song, A. Khosla, F. Yu, L. Zhang, X. Tang, and J. Xiao.  3d shapenets: A deep representation for volumetric shapes.  In CVPR, 2015.
  * [31] L. Yi, H. Su, X. Guo, and L. Guibas.  Syncspeccnn: Synchronized spectral cnn for 3d shape segmentation.  CVPR, 2017.

Appendix

In this supplementary material, we first validate our method’s generalization
ability by testing the proposed models on KITTI dataset without re-training.
We then provide more details about the lidar flow task. We also show
additional results for all tasks, with analysis of failure modes. Lastly, we
conduct a point cloud classification experiment using the deep parameteric
continuous convolutional networks. In the supplementary video, we show our
results on semantic segmentation and lidar flow results over a sequential
data. We also show the generalization ability of the proposed method by
training on our dataset and test it on KITTI as well as a truck sequence.

##  Appendix A Generalization

We show the generalization ability of our proposed model by training over one
dataset and test it over another in our supplementary video. To be specific,
we have used the following configurations:

  * •

Train our proposed semantic labeling network on the driving scene data
(several north America cities), and test it on KITTI (Europe).

  * •

Train our proposed semantic labeling network on the driving scene data (non-
highway road), and test it on a lidar sequence mounted on top of a high-way
driving truck (highway road).

  * •

Train our proposed lidar flow network on the driving scene data (several north
America cities), and test it on KITTI (Europe).

|  |
---|---|---
Figure 9: Lidar Flow Results on KITTI Dataset |  |
---|---|---
|  |
Figure 10: Semantic Segmentation on KITTI Dataset

Under all the settings, our algorithm is able to generalize well. Fig. 9 shows
our lidar flow model’s performance on KITTI. As shown in Fig. 9, our Lidar
flow model generalizes well to the unseen KITTI dataset. From left to right,
the figures show the most common scenarios with moving, turning, and
stationary ego-car. The model produces covincing flow predictions in all three
cases. Fig. 10 includes some additional segmentation results on KITTI.

For more results over several sequences, please refer to our supplementary
video.

##  Appendix B Lidar Flow Data and Analysis

#### Ground-truth Generation

In this paragraph, we describe how we generate the ground-truth lidar flow
data in detail. For each consecutive frame, we first get the global vehicle
pose transform from frame 0 to frame 1:
$\mathbf{R}_{\mathrm{ego}},\mathbf{t}_{\mathrm{ego}}$ with the help of
additional sensors and prior information. This global vehicle pose transform
represents how far away the vehicle moves and how the vehicle turns. This
localization accuracy is at centi-meter scale. Therefore, the motion per each
static point is:

$\mathbf{f}_{\mathrm{static-
gt}}^{(0)}=\mathbf{R}_{\mathrm{ego}}^{T}(\mathbf{x}^{(0)}_{static}-\mathbf{t}_{\mathrm{ego}})-\mathbf{x}^{(0)}_{static}$

where $\mathbf{f}_{\mathrm{gt}}^{(k)}$ is the ground-truth flow at the frame
$k$ under the ego-car centered coordinate, $\mathbf{x}^{(k)}$ is the point’s
location at the frame $k$ under the ego-car centered coordinate.

For dynamic objects in the scene, e.g. other vehicles and pedestrians, the
motion between each lidar frame in the vehicle coordinate is not only due to
the self-driving car’s ego-motion. The movement of the dynamic objects
themselves are also contributing to the motion. In this project, we assume
rigid motion for all the objects. The labeling of the dynamics objects include
two steps. Firstly, using the global pose that we get, we visualize the point
cloud of 3D objects from two frames in the same reference coordinate and label
the pose changes $\mathbf{R}_{\mathrm{obj}},\mathbf{t}_{\mathrm{obj}}$ between
the objects at the different time. Secondly, both ego-motion and object motion
are considered in order to generate the ground-truth flow vector:

$\mathbf{f}_{\mathrm{dynamic-
gt}}^{(0)}=\mathbf{R}_{\mathrm{obj}}^{T}(\mathbf{R}_{\mathrm{ego}}^{T}(\mathbf{x}^{(0)}_{dynamic}-\mathbf{t}_{\mathrm{ego}})-\mathbf{t}_{\mathrm{obj}})-\mathbf{x}^{(0)}_{dynamic}$

Please refer to Fig. 11 for an illustration.

Figure 11: Flow data generation. The source of motion comes from two
components: motion of the ego-car and motion of the dynamic objects.

#### Ground-truth Motion Analysis

We also conduct an analysis over the ground-truth motion distribution. In Fig.
12 we show the 2D histogram of the GT 3D translation component along $x$ and
$y$ axis respectively. We also show the motion distribution across different
object types, e.g. static background, vehicle and pedestrian. As we can see,
different semantic types have different motion patterns. And the heaviest
density of distribution is on the y-axis, which suggests the forward motion is
the major motion pattern of our ego-car.

Figure 12: Ground-truth Motion Distribution of the Lidar Flow Dataset (unit in
meters)

#### Ground-truth Validation and Visualization

We validate the quality of our ground-truth motion labels by overlaying target
frame points ($\mathbf{x}^{(1)}$) and source frame points warped with ground-
truth motion ($\mathbf{x}^{(0)}+\mathbf{f}_{\mathrm{gt}}^{(0)}$). Fig. 13
shows overlays of entire scenes and Fig. 14 shows overlays of individual
dynamic objects. For vehicles, points align perfectly across frames. For
pedestrians, the correspondence is also near perfect: the only discrepancy is
caused by non-rigid motion (e.g. changing posture).

|  |
---|---|---
|  |
Figure 13: Flow Ground Truth Overlay of Entire Scene. Yellow: target frame, purple: warped source frame. |  |  |  |
---|---|---|---|---
|  |  |  |
Figure 14: Flow Ground Truth Overlay of Individual Dynamic Objects. Green:
target frame, red: warped source frame.

##  Appendix C More Results

In this section, we show additional qualitative results of the proposed
algorithm over all the tasks.

###  C.1 Semantic Segmentation for Indoor Scenes

Fig. 15 and Fig. 16 show more qualitatitive results over the stanford dataset.
As the figure shown, in most cases our model is able to predict the semantic
labels correctly.

|  |
---|---|---
|  |
|  |
|  |
|  |
|  |
Input | GT | Ours
Figure 15: Semantic Segmentation Results on Stanford Indoor Dataset |  |
---|---|---
|  |
|  |
|  |
|  |
|  |
Input | GT | Ours
Figure 16: Semantic Segmentation Results on Stanford Indoor Dataset

###  C.2 Semantic Segmentation for Driving Scenes

Fig. 17 shows additional results for semantic labeling in driving scenes. As
shown, the results capture very small dynamics, _e.g_. pedestrians and
bicyclists. This suggests our model’s potential in object detection and
tracking. The model is also able to distinguish between road and non-road
through lidar intensity and subtle geometry structure such as road curbs. This
validates our model’s potential in map automation.

More specifically, we see that most error occur on road boundaries (bright
curves in error map).

|  |
---|---|---
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
Ground Truth | Ours | Error Map
Figure 17: Semantic Segmentation Results on Driving Scenes Dataset

###  C.3 Lidar Flow

We show additional results on Lidar flow estimation in Fig. 18. Unlike the
visualization in the main submission, we visualize the colored vector in order
to better depicts the magnitudes of the motion vector. As shown in the figure,
our model is able to capture majority flow field. The majority of the error
happens at the object boundary. This suggests that a better support domain
that includes both space and intensity features could be potentially used to
boost performance.

|  |
---|---|---
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
Ground Truth | Ours | Error Map
Figure 18: Lidar Flow Results on Driving Scenes Dataset

##  Appendix D Activations

We visualize the activation maps of the trained PCCN network over a single
lidar frame from the driving scene dataset for segmentation. Fig. 19 depicts
the activation map at layer 1 of PCCN. As we can see, at the early conv layer
the method mainly captures low-level geometry details, _e.g_. the z
coordinate, the intensity peak, _etc_. Fig. 20 shows the activation map at
layer 8 of PCCN. The conv layers begin to capture information with more
semantic meaning, e.g. the road curb and the dynamic objects.

|  |
---|---|---
|  |
Figure 19: Activation Map of PCCN at Layer 1 |  |
---|---|---
|  |
Figure 20: Activation Map of PCCN at Layer 8

##  Appendix E Point Cloud Classification

To verify the applicability of the proposed parameteric continuous convolution
over global prediction task, we conduct a simple point cloud classification
task on the ModelNet40 benchmark. This dataset contains CAD models from 40
categories. The state-of-the-art and most representative algorithms conducted
on ModelNet40 are compared [28]. We randomly sampled 2048 points for each
training and testing sample over the 3D meshes and feed the point cloud into
our neural network. The architecture contains 6 continuous convolution layers
with 32-dimensional hidden features, followed by two layers with
128-dimensions and 512 dimensions respectively. The output of the last
continuous convolution layer is fed into a max pooling layer to generate the
global 512-dimensional feature, followed by two fc layers to output the final
logits. Tab. 4 reports the classification performance. As we can see in the
table, the performance is comparable with PointNet and slightly below
PointNet++. Here we use a naive global max pooling to aggregate global
information for our method. We expect to achieve better results with more
comprehensive and hierachical pooling strategies.

Table 4: ModelNet40 Point Cloud Classification Method | Input | Accuracy
---|---|---
MVCNN [28] | Multi-view Image | 90.1%
3DShapeNet [30] | Volume | 84.7%
VoxNet [18] | Volume | 85.9%
Subvolume [21] | Volume | 89.2%
ECC [27] | Point | 87.4%
PointNet vanilla [20] | Point | 87.2%
PointNet [20] | Point | 89.2%
PointNet++ [22] | Point | 91.9%
Ours | Point | 88.9%


# Population dynamics in the triplet annihilation model with a mutating
reproduction rate

Ronald Dickman,111email<EMAIL_ADDRESS>Departamento de Física and
National Institute of Science and Technology for Complex Systems,
ICEx, Universidade Federal de Minas Gerais,
C. P. 702, 30123-970 Belo Horizonte, Minas Gerais - Brazil

###### Abstract

I study a population model in which the reproduction rate $\lambda$ is
inherited with mutation, favoring fast reproducers in the short term, but
conflicting with a process that eliminates agglomerations of individuals. The
model is a variant of the triplet annihilation model introduced several
decades ago [R. Dickman, Phys. Rev. B 40, 7005 (1989)] in which organisms
(“particles”) reproduce and diffuse on a lattice, subject to annihilation when
(and only when) occupying three consecutive sites. For diffusion rates below a
certain value, the population possesses two “survival strategies”: (i) rare
reproduction ($0<\lambda<\lambda_{c,1}$), in which a low density of diffusing
particles renders triplets exceedingly rare, and (ii) frequent reproduction
($\lambda>\lambda_{c,2}$). For $\lambda$ between $\lambda_{c,1}$ and
$\lambda_{c,2}$ there is no active steady state. In the rare-reproduction
regime, a mutating $\lambda$ leads to stochastic boom-and-bust cycles in which
the reproduction rate fluctuates upward in certain regions, only to lead to
extinction as the local value of $\lambda$ becomes excessive. The global
population can nevertheless survive due to the presence of other regions, with
reproduction rates that have yet to drift upward.

##  I Introduction

My regrettably few encounters with Dietrich Stauffer left an important
imprint. He always impressed with his unconventional, even provocative
attitude, his honesty, generosity, and sense of humor. I first met Prof. Dr.
Stauffer at a Rutgers Meeting in the early ’80s. I’d been enumerating lattice
animals as part of my PhD project on nucleation in lattice gases. When he
asked what computer I used for the task, I simply gestured to indicate that
this was hand enumeration. But a year later, having exhausted the
possibilities of this approach (and/or my patience with filling notebooks with
strange little diagrams), Stauffer’s question led me to develop Monte Carlo
methods for counting lattice animalsanimals .

Toward the end of the 1980s, on a visit to Germany at the kind invitation of
Prof. Kurt Binder, I met Stauffer at his lab in Jülich, and described my
ongoing study of the triplet annihilation model, the subject of the present
work. Hearing my explanation that particles could only be annihilated in
groups of three, he turned and said, “Oh, you’re studying guerrilla warriors
who spread out to evade detection from the air.” This colorful and somewhat
poignant interpretation stuck with me, as did his advice (in reply to a letter
sent shortly before defending my PhD, in which I confessed to having become
“addicted” to Monte Carlo simulation), that I simply “enjoy my addiction.” And
so I have. In a handful of encounters in Europe and Brazil, we chatted about
science, politics and other matters. (I’m afraid my answers were sometimes
disappointing: Seen any good movies recently? - My tastes were too art-house.
Your favorite football team? - I had none.)

Dietrich Stauffer made major contributions to percolation theory, so it is
fitting that the triplet annihilation model studied here belongs to the
(rather broad) class of models with a directed percolation (DP)-like phase
transition. In the stochastic particle systems of interest here, the directed
spatial axis of DP corresponds to time; the (continuous) phase transition is
between an active and an absorbing state. The simplest example is the contact
process (CP) harris74 ; liggett ; marro ; odor07 ; henkel , in which each site
of a lattice is either vacant or occupied by a single particle. The latter die
at unit rate, and attempt to reproduce at rate $\lambda$, choosing a
neighboring site at random and adding a particle if the chosen neighbor is
vacant. In the infinite-size limit there is a transition between an active
phase (for $\lambda>\lambda_{c}$) with a nonzero particle density and an
absorbing phase (all sites vacant). To the CP dynamics, the triplet
annihilation model tam89 ; tam90  adds diffusion (nearest-neighbor hopping)
and a modified death process, in which only groups of three particles are
removed. The result is that (for not too rapid diffusion) an active phase
exists for small $\lambda$ (with a very low abundance of triplets) as well as
for large $\lambda$, but not for intermediate values. Here, I explore the
effect of a reproduction rate, $\lambda_{i}$ for each particle $i$, that is
inherited from its parent with some random variation. Depending on the degree
of variability, this can lead to a population with minimal reproduction, to
extinction, or to long-time survival punctuated by a series of local
population explosions and collapses.

The remainder of this paper is organized as follows. In the next section I
define the model and describe the simulation algorithm. Sec. III reports
simulation results, followed by a summary and discussion in Sec. IV.

##  II Model

###  II.1 Triplet annihilation model

The triplet annihilation model (TAM) tam89 ; tam90  is a stochastic particle
system on a lattice, a ring of $L$ sites in the present study dimnote . Each
site can be either vacant or occupied by a (single) particle. Transitions
between configurations occur through three basic processes:

  * •

Creation: Each particle attempts independently, at rate $2\lambda$, to create
a new particle at its nearest-neighbor to the left or right. The attempt is
successful if the chosen neighbor is unoccupied.

  * •

Hopping: Each particle attempts to hop, at rate $2d$, to its nearest-neighbor
to the left or right. The attempt is successful if the chosen neighbor is
unoccupied.

  * •

Annihilation: If sites $i-1$, $i$ and $i+1$ are all occupied, we say there is
a triplet centered at site $i$. Triplets are annihilated a rate of unity. (All
three particles are removed in this event.)

From this definition it is evident that the configuration with all sites
vacant is absorbing. (The parameterization of the rates used here differs from
that used in the earlier studies, although it has no qualitative effect on the
phase diagram.)

Note that hopping tends to destroy triplets, whereas creation tends to
generate them, raising the possibility of subtle interactions between these
processes. The original motivation for studying the model was to investigate
what effect such conflicting tendencies might have on the phase diagram and
critical behavior.aside

The phase diagram, shown schematically in Fig. 1 is indeed intriguing. In
simpler lattice birth-and-death models such as the contact process, sustained
activity (i.e., a nonzero particle density, $\rho>0$), requires a reproduction
rate greater than some value $\lambda_{c}>0$. The presence, in the TAM, of
activity for arbitrarily small $\lambda$ is due, of course, to the death rate
being proportional to the density of triplets, $\rho_{3}$, rather than of
particles, $\rho$. In the presence of hopping $(d>0)$, there is an interval
$(0,\lambda_{-}(d))$ of reproduction rates such that $\rho_{3}<\\!\\!<\rho$,
permitting survival of activity. Of particular interest in the present work is
the portion of the active phase bounded by $\lambda=0$ and the small-$\lambda$
branch of the phase boundary, $\lambda_{-}(d)$, because in this region,
paradoxically, increasing the reproduction rate leads to extinction, due to
the increased abundance of triplets. An example of this behavior is shown in
Fig. 2. The quasistationary particle density QS  peaks near $\lambda\simeq
0.005$ and falls rapidly for higher values. For this hopping rate
($d=1.2346$), the upper critical reproduction rate, $\lambda_{+}\simeq 1.17$,
much larger than the values considered here; preliminary studies yield
$\lambda_{-}\simeq 0.03$.

Although critical behavior is not the principal subject of the present study,
it is worth noting that the entire phase boundary appears to be a line of
critical points in the directed percolation universality class, and that the
particle density $\rho\sim\lambda^{1/2}$ as $\lambda\to 0$. This property, as
well as the qualitative form of the phase diagram, is captured by mean-field
theory at the pair level tam89 ; tam90 . Several works have investigated the
phase diagram and closely related models. Katori and Konno katori  established
some general properties of systems like the TAM, while Poland poland  used
time-series expansions to investigate population dynamics in the TAM and
related models. A model with annihilation of quadruplets was analyzed by Ódor
odor2004 ; odor2004a .

Figure 1: Schematic of the triplet annihilation model phase diagram in the
reproduction rate-hopping rate ($\lambda-d$) plane. The phase boundary is a
line of continuous phase transitions between the active and absorbing phases.
For hopping rates $d>d^{*}$, there is activity for any $\lambda>0$, while for
$d<d^{*}$, the active phase exists for large $\lambda$, but also for very
small $\lambda$. The present study focuses on the region bounded by
$\lambda=0$ and the small-$\lambda$ branch of the phase boundary,
$\lambda_{-}(d)$.  Figure 2: Quasistationary particle density $\rho$ versus
reproduction rate $\lambda$ in the TAM, for hopping parameter $d=1.2346$,
system size $L=10^{4}$. Error bars are smaller than the symbols. Inset: the
same data plotted versus $\sqrt{\lambda}$.

###  II.2 TAM with mutating reproduction rate (MTAM)

In the present study, the reproduction rate is no longer a fixed parameter.
Instead, each particle $i$ is assigned a value $\lambda_{i}$ when created, and
carries this value during its lifetime, as it hops on the lattice. If particle
$i$ creates a new particle $j$, then $\lambda_{j}$ is inherited from the
parent with random variation:

$\lambda_{j}=\max\left[\lambda_{i}\\!+\\!\Delta,\,\lambda_{min}\right],$ (1)

where $\Delta$ is a random variable with zero mean, and $\lambda_{min}$ is a
small basal reproduction rate. I consider both additive and multiplicative
schemes. In the former, $\Delta$ is simply chosen from a distribution uniform
on $[-\gamma,\gamma]$. In the multiplicative case, $\Delta=\xi\lambda_{i}$,
where $\xi$ is uniform on $[-\gamma,\gamma]$. Each realization of the process
begins with all sites occupied, and $\lambda_{i}=\lambda_{0},\,\forall
i=1,...,L$. In what follows I refer to this as the mutating triplet
annihilation model (MTAM).

###  II.3 Simulation algorithm

I developed an event-driven (rejection-free) Monte Carlo simulation algorithm
for the present study. The basic requirements this type of simulation are
lists of the possible transitions and pointers to facilitate updating the
lists. In the TAM simulation one must maintain a list of triplets and another
list of nearest-neighbor particle-hole pairs (PHPs). One maintains, in
addition, vectors carrying (i) the occupation status of each site; (ii)
whether sites $i$ and $i+1$ belong to a PHP, and if so, the position on the
PHP list, and which of the two sites is occupied; (iii) whether site $i$ is
the central site of a triplet, and if so, the corresponding position on the
triplet list. In the case of mutating reproduction rates, a fourth vector
carries $\lambda_{i}$ associated with each occupied site $i$.

Given the numbers, $N_{trp}$ and $N_{php}$ of triplets and PHPs, respectively,
the total transition rate out of the current configuration is,

$R=N_{trp}+(d+\lambda)N_{php}\equiv R_{ann}+R_{hop}+R_{rep},$ (2)

where in the second equality we define the total rates of annihilation,
hopping and reproduction. The next transition is taken to be of type $j$
(annihilation, hopping or reproduction), with probability $P_{j}=R_{j}/R$.
Then the specific event on the corresponding list is selected (with uniform
probability over the list), the new configuration constructed, the lists
updated, and the time advanced so: $t\to t+1/R$.

In the case of mutating reproduction rates, i.e., the MTAM, we must alter the
expression for the total creation rate to $R_{rep}=\Lambda$, where

$\Lambda=\sum_{i=1}^{N_{php}}\lambda_{j(i)},$ (3)

where the sum is over PHPs and $j(i)$ denotes the occupied site in the $i$th
PHP. In this case, of course, the PHP for creation must be chosen with a
probability proportional to the associated $\lambda$ value. This is done by
letting $\tilde{\Lambda}=\Lambda z$, where $z$ is a random number uniform on
$[0,1]$, and then evaluating the partial sums $\Lambda_{k}$, with
$\Lambda_{1}=\lambda_{j(1)}$ and,

$\Lambda_{k}=\Lambda_{k-1}+\lambda_{j(k)},$ (4)

until $\Lambda_{k}>\tilde{\Lambda}$, at which point creation is realized in
the $k$th PHP. While more efficient methods for choosing the PHP exist, this
approach is simple and unlikely to increase run times significantly since
reproduction counts for only a small fraction of events and the number of PHPs
is not very large in the regime of interest.

##  III Results

Given the sensitivity of the TAM to small changes in the reproduction rate,
and the counterintuitive nature of the response, it is of interest to examine
how a population evolves under a mutating rate.

In this section I report preliminary results for the time evolution and
spatial structure of the MTAM. The fixed simulation parameters are $d=1.2346$,
initial creation rate $\lambda_{0}=0.005$, basal creation rate
$\lambda_{min}=0.001$, and system size $L=2\,000$. Studies run to a maximum
time of $10^{6}$ unless the absorbing state is attained prior to this time.
(For these parameters, but without mutation, the probability of a realization
surviving to time $10^{6}$ is about 0.86, and the quasistationary particle
density is $\rho=0.0331(1)$.)

To begin I consider the effect of multiplicative noise in the inherited
reproduction rate, varying the noise intensity $\gamma$ between $0.005$ and
0.25. Most realizations survive until the maximum time: For $\gamma=0.025$,
for example, the survival probability is about 0.87. Time series for the
particle density show significant fluctuations, with or without mutation, as
shown in Fig. 3. This figure also shows the time series for
$\bar{\lambda}(t)$, the instantaneous average of $\lambda$ over all particles,
as well as the instantaneous maximum, $\lambda_{max}(t)$. These quantities
exhibit a pattern of slow growth followed by a sudden crash. The mean
reproduction rate remains “pinned” near its basal rate of $10^{-3}$ for long
periods. The inset shows that the variations in the standard deviation of
$\lambda$ mirror those of $\bar{\lambda}(t)$. The particle density also
exhibits crashes (followed by recovery, at least in this example), but the
relation between crashes in population and in $\lambda$ is unclear, since one
may occur without the other. Similar examples of time series with mutating
reproduction rates are shown in Figs. 4 and 5.

Figure 3: Time series of the particle density $\rho(t)$ and reproduction rate
in the TAM and MTAM with hopping parameter $d=1.2346$, and system size
$L=2000$. Lower panel: $\lambda=0.005$, without mutation (TAM). Middle panel:
$\rho(t)$ for the MTAM with mutation parameter $\gamma=0.1$ (multiplicative
case). Upper panel: $\bar{\lambda}(t)$(lower) and $\lambda_{max}(t)$ (upper)
in the same realization as the middle panel. Inset: comparison of time series
for $\bar{\lambda}(t)$ (upper curve) and the standard deviation of the set of
$\lambda_{i}$ values at time $t$ (lower).  Figure 4: MTAM: Time series of the
particle density $\rho(t)$ and reproduction rate as in Fig. 3, except
$\gamma=0.025$.  Figure 5: MTAM: Time series of the particle density $\rho(t)$
and reproduction rate, parameters as in Fig. 4.  Figure 6:  MTAM: Time series
of the particle density $\rho(t)$ and mean reproduction rate in with
multiplicative noise, hopping parameter $d=1.2346$, initial reproduction rate
$\lambda_{0}=1.15$, $\gamma=5\times 10^{-4}$, and system size $L=2000$.

For an initial reproduction rate $\lambda_{0}$ close to the upper phase
boundary, $\lambda_{+}(d)$, mutations lead, as expected, to steady growth in
$\bar{\lambda}$ and the particle density. An example is shown in Fig. 6, for
$\lambda_{0}=1.15$, close to $\lambda_{+}=1.17$.

Time series for additive noise are generally similar to those for the
multiplicative case, if we use noise intensities $\gamma$ about two orders of
magnitude smaller. (Recall that in the multiplicative case, the noise term is
scaled by $\lambda\sim 0.005$.) A typical example, for $\gamma=5\times
10^{-4}$, is shown in Fig. 7. “Pinning” of $\bar{\lambda}$ near
$\lambda_{min}$ is not observed under additive noise, since the magnitude of
mutations (i.e., of $\Delta\lambda$) is independent of $\lambda_{i}$. For
$\gamma\geq 10^{-3}$ most realizations with additive noise do not survive
until time $10^{6}$.

An example of a population that goes extinct rapidly is shown in Fig. 8.
Inspection of the time series shows that even the exceptionally high values of
$\bar{\lambda}$ for $t\simeq 1.7\times 10^{4}$ do not lead to an immediate
population crash. Partial extinction events near $t=1.8$ and $2.1\times
10^{4}$ appear to cull the rapidly reproducing subpopulations. The population
and reproduction rate continue to fluctuate, until, at around $t=2.8\times
10^{4}$, the system attains a regime with small $\rho$ and small
$\bar{\lambda}$. This combination makes appearance of individuals with larger
$\lambda$ values unlikely, and the population goes extinct some time later.
This example shows that although high reproduction rates play a role in
extinction, they may not the direct or immediately cause. Further analyses
will be required to assess the generality of this scenario.

Figure 7: MTAM: Time series of the particle density $\rho(t)$ and reproduction
rate, but with additive noise of intersity $\gamma=5\times 10^{-4}$. Other
parameters as in Fig. 3.  Figure 8: MTAM: Time series of the particle density
$\rho(t)$ and reproduction rate, but with additive noise of intersity
$\gamma=0.005$. Other parameters as in Fig. 3. The population goes extinct at
$t=3.192\times 10^{4}$.

###  III.1 Spatiotemporal evolution

An interesting aspect of the population dynamics is the appearance of
localized subpopulations with a more or less uniform reproduction rate. The
average reproduction rate within such a group may fluctuate upward until it
exceeds a “safe” value, and the group dies out due to overproduction or
triplets. For parameter sets allowing long-time survival, it appears that when
such a group goes extinct, there remain several other subpopulations,
typically with small reproduction rates. An example of such an evolution is
shown in Fig. 9, which covers the period $t=4.0-4.8\times 10^{5}$ shown in
Fig. 5, during which the population grows and then crashes. Localized groups
with restricted variations in $\lambda$ are evident; one does not find
individuals with markedly different $\lambda$ values in the same region at the
same time. In the region extending from $x\simeq 500-1000$, the mean value of
$\lambda$ initially drifts upward. The region then becomes nearly empty, to be
reoccupied at the final time shown, presumably by particles that have diffused
in from nearby groups. By contrast, the region between $x=1500$ and 200 (i.e.,
2200 under the periodic boundary) exhibits only modest variations in density
and reproduction rate. One may conjecture that this dynamic homogeneity is
what allows for long-term survival in a population undergoing “boom-and-bust”
cycles in localized regions. (Survival, one might say, lies in diversity.)
Further studies, including the effect of varying system sizes, will be needed
to test these ideas.

Figure 9: MTAM: Spatiotemporal evolution during a period of population growth
and decay. Each type of symbol (see legend) corresponds to a specific time in
the time series of Fig. 5. Each symbol represents the position ($x$) and
reproduction rate ($\lambda$) of a single particle.

##  IV Summary and discussion

I study the population dynamics of the triplet annihilation model (TAM) in a
regime in which only a narrow range of (quite small) reproduction rates assure
survival. The the mutating version, MTAM, individual reproduction rates are
inherited with random mutation, leading to complex spatiotemporal dynamics.
Much remains to be analyzed following this preliminary study, including
analytic approximations such as deterministic and/or stochastic partial
differential equations.

There are further variations to be considered. In the present context, for
example, individuals with an optimal reproduction rate (in the sense of Fig.
2) may be seen as cooperators, and rapid reproducers as cheaters or betrayers.
This raises the possibility of introducing mechanisms that foster cooperation
in an evolutionary-game context.

Stochastic particle systems offer fascinating challenges for physicists and
mathematicians. Many of their interpretations, as models of populations,
epidemics, social systems, or even guerrilla warriors, are imperfect, serving
nevertheless as intuitively attractive metaphors for processes in the real
world. I wonder what fanciful interpretation Dietrich Stauffer would have
offered for the model considered here.

Acknowledgment

This work was supported by CNPq, Brazil, under grant No. 303766/2016-6.

## References

  * (1) R. Dickman and W. C. Schieve, “Lattice animal specific heats and the collapse of branched polymers,” J. Phys. France 45, 1727 (1984).
  * (2) T. E. Harris, “Contact Interactions on a Lattice,” Ann. Probab. 2, 969 (1974).
  * (3) T. M. Liggett, Interacting particle systems Springer-Verlag, Berlin (2005).
  * (4) J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999), Ch. 6.
  * (5) G. Ódor, Universality In Nonequilibrium Lattice Systems: Theoretical Foundations (World Scientific,Singapore, 2007)
  * (6) M. Henkel, H. Hinrichsen and S. Lubeck, Non-Equilibrium Phase Transitions Volume I: Absorbing Phase Transitions (Springer-Verlag, The Netherlands, 2008).
  * (7) R. Dickman, “Universality and diffusion in nonequilibrium critical phenomena,” Phys. Rev. B 40, 7005 (1989).
  * (8) R. Dickman, “Nonequilibrium critical behavior of the triplet annihilation model,” Phys. Rev. A 42, 6985 (1990).
  * (9) Thus the present study, like those of Refs. tam89 , tam90 , and poland , is restricted to one dimension. The behavior of the model in higher dimensions, or on fractal substrates, remains as an interesting topic for future study.
  * (10) A suitable interpretation of the TAM might be that of a mobile population subject to predation when its local density exceeds a certain threshold. See J. D. Murray, Mathematical Biology I. (Springer, New York, 1993) §1.7 and references therein.
  * (11) M. Katori and N. Konno, “On the extinction of Dickman’s reaction- diffusion processes,” Physica A 186, 578 (1992).
  * (12) D. Poland, “Time-power series study of the triplet annihilation model and other cooperative growth models,” Physica A 193, 1 (1993).
  * (13) G. Ódor, “Phase transitions of the binary production 2A $\to$ 3A, 4A $\to$ X model,” Phys. Rev. E 69, 036112 (2004).
  * (14) G. Ódor, “Critical behavior of the two-dimensional 2A $\to$ 3A, 4A $\to$ $\not{\\!0}$ binary system,” Phys. Rev. E 70, 026119 (2004).
  * (15) A quasistationary mean is conditioned on long-term survival, or more formally, the limiting ($t\to\infty$) value of a mean taken over the set of realizations that survive until time $t$ or longer. For models with an absorbing state, the quasistationary mean is the closest analog to the usual notion of a stationary mean.