# Concave-Convex PDMP-based sampling

Matthew Sutton
Centre for Data Science, Queensland University of Technology
and
Paul Fearnhead
Department of Mathematics and Statistics, Lancaster University This research
was supported by EPSRC grant EP/R018561

###### Abstract

Recently non-reversible samplers based on simulating piecewise deterministic
Markov processes (PDMPs) have shown potential for efficient sampling in
Bayesian inference problems. However, there remains a lack of guidance on how
to best implement these algorithms. If implemented poorly, the computational
costs of simulating event times can out-weigh the statistical efficiency of
the non-reversible dynamics. Drawing on the adaptive rejection literature, we
propose the concave-convex adaptive thinning approach for simulating a
piecewise deterministic Markov process (CC-PDMP). This approach provides a
general guide for constructing bounds that may be used to facilitate PDMP-
based sampling. A key advantage of this method is its additive structure -
adding concave-convex decompositions yields a concave-convex decomposition.
This facilitates swapping priors, simple implementation and computationally
efficient thinning. In particular, our approach is well suited to local PDMP
simulation where known conditional independence of the target can be exploited
for potentially huge computational gains. We provide an R package for
implementing the CC-PDMP approach and illustrate how our method outperforms
existing approaches to simulating events in the PDMP literature.

##  1 Introduction

Monte Carlo methods based on continuous-time Markov processes have shown
potential for efficient sampling in Bayesian inference challenges (Goldman and
Singh, 2021; Fearnhead et al., 2018). These new sampling algorithms are based
on simulating piecewise-deterministic Markov processes (PDMPs). Such processes
evolve deterministically between random event times with possibly random
dynamics at event times (see Davis, 1993, for an introduction to PDMPs). Monte
Carlo methods based on simulating a PDMP with a desired target distribution of
interest were first proposed in the computational physics literature where
they were motivated as examples of non-reversible Markov processes (Peters and
de With, 2012; Michel et al., 2014). These ideas transitioned to the
statistics community as an alternative to traditional Markov chain Monte Carlo
(MCMC) samplers. Popular algorithms in this development include the Bouncy
Particle Sampler (Bouchard-Côté et al., 2018) and the Zig-Zag sampler
(Bierkens et al., 2019) amongst others (Vanetti et al., 2017; Wu and Robert,
2020; Michel et al., 2020; Bierkens et al., 2020).

There is substantial theoretical and empirical evidence to support the claim
that non-reversible MCMC samplers offer more efficient sampling as they reduce
the chance of the sampler moving back to areas of the state space that have
recently been visited (Diaconis et al., 2000; Bierkens and Roberts, 2017;
Bouchard-Côté et al., 2018). Samplers based on PDMPs simulate from a density
$\pi(\bm{\theta})$ by introducing a velocity component $\bm{v}$ which is used
to evolve the position of the state, $\bm{\theta}$, with deterministic
dynamics. At stochastic event times the velocity component is updated using a
Markov kernel after which the process continues with the new velocity. The
main challenge to implementing a PDMP-based sampler is the simulation of the
random event times. This involves simulating the next event time in a time-
inhomogeneous Poisson process with rate function that depends on the current
position of the sampler.

Simulation of event times in a Poisson process is often facilitated through a
method known as thinning (Lewis and Shedler, 1979). The aim of this paper is
to provide a general framework to aid practitioners in implementing PDMP-based
samplers. Specifically we introduce the concave-convex thinning approach to
facilitate simulating a PDMP-based sampler (CC-PDMP). This method can be
applied whenever the rate function of the PDMP can be decomposed into the sum
of concave and convex functions, and can be used to facilitate thinning from
polynomial rate functions, allowing for broad applicability of the method. We
discuss the efficiency of the method compared to alternative approaches in the
literature for exact sampling of event times.

Related to our work is the concave-convex adaptive rejection sampler of Görür
and Teh (2011) which uses the same upper-bounding approach to construct bounds
on the log density within a rejection sampler that draws independent samples
from univariate density functions. In contrast, by applying this technique
within PDMP sampling, the bounds are constructed on the gradient of the log
density (i.e. the rate) and used within the non-reversible Markov Process
framework to generate samples from multivariate densities. Moreover, the CC-
PDMP framework allows for both subsampling and local updating schemes to
facilitate these methods for high dimensional sampling. A short comparison
between concave-convex adaptive rejection sampling and CC-PDMP for sampling
from a univariate density is given in Appendix C.

Some authors have proposed using automatic but approximate methods for
simulating a PDMP-based sampler. Among these approaches, Cotter et al. (2020)
propose numerical integration and root-finding algorithms to facilitate
simulation. Others have considered using approximate local bounds which can be
simulated exactly (Goldman and Singh, 2021; Goan and Fookes, 2021; Pakman et
al., 2017). Both of these approaches sacrifice exact sampling of the posterior
and involve a trade-off between the computational cost and the approximation
of the sampling distribution.

The paper is organised as follows. In Section 2, we introduce the technical
details of simulating from a PDMP-based sampler and the details of some
popular samplers. Section 3 reviews the literature on thinning for PDMP-based
samplers. We introduce the concave-convex PDMP approach for adaptive thinning
in Section 4. Empirical evaluation of CC-PDMP, and comparison to existing
methods, is given in Section 5. We conclude with discussion of limitations and
extensions that are possible for the method. Code for implementing our method
and replicating our results is available at https://github.com/matt-
sutton/ccpdmp.

##  2 Sampling using a PDMP

###  2.1 Piecewise deterministic Markov processes

Before we explore PDMP-based samplers, we first review PDMPs. A PDMP is a
stochastic process, and we will denote its state at time $t$ by $\bm{z}_{t}$.
There are three key components that define the dynamics of a PDMP:

  1. 1.

A set of deterministic dynamics defined as $\bm{z}_{t+s}=\psi(\bm{z}_{t},s)$.

  2. 2.

Event times that are driven by an inhomogeneous Poisson process with rate
$\lambda(\bm{z}_{t})$.

  3. 3.

A Markov transition kernel $q(\bm{z}|\bm{z}_{t-})$ where
$\bm{z}_{t-}:=\lim_{s\uparrow t}\bm{z}_{s}$.

A PDMP with these specifications will evolve according to its deterministic
dynamics between event times and at event times will update according to the
Markov transition kernel $q(\bm{z}|\bm{z}_{t-})$. The event times must be
simulated from a Poisson process with rate $\lambda(\bm{z}_{t})$, which
depends on the position of the process at time $t$. The result of simulating
this process is a PDMP skeleton $\\{(t_{k},\bm{z}_{t_{k}})\\}_{k=1}^{n}$ where
for $k=1,\ldots,n$, $t_{k}$ denotes the $k$th event time and we have stored
the state of the process immediately after this event time, $\bm{z}_{t_{k}}$.
Given this skeleton it is possible to fill in the values of the state between
the events using the deterministic dynamics.

###  2.2 PDMP-based samplers

Assume we wish to construct a PDMP process to sample from a target
distribution of interest, $\pi(\bm{\theta})$. Current PDMP-based samplers are
defined on an augmented space $\bm{z}\in\mathcal{E}$, where
$\mathcal{E}=\mathcal{X}\times\mathcal{V}$, that can be viewed as having a
position, $\bm{\theta}_{t}$, and a velocity, $\bm{v}_{t}$. As described below,
the dynamics of the PDMP can be chosen so that the PDMP’s invariant
distribution has the form $\nu(\bm{z})=\pi(\bm{\theta})p(\bm{v}|\bm{\theta})$,
for some conditional distribution of the velocities, $p(\bm{v}|\bm{\theta})$.
As a result, the $\bm{\theta}$-marginal of the invariant distribution is the
target distribution we wish to sample from. We will consider target
distributions of the form

$\pi(\bm{\theta})\propto\exp(-U(\bm{\theta}))$

where $\bm{\theta}\in\mathbb{R}^{p}$ and $U(\bm{\theta})$ is the associated
potential. An important feature of PDMP samplers is that they need to know the
target distribution only up to a constant of proportionality.

If we simulate the PDMP, and the process converges to stationarity, then we
can use the simulated skeleton of the PDMP to estimate expectations with
respect to the target distribution in two ways. Assume we are interested in a
Monte Carlo estimate of $\int g(\bm{z})\nu(\bm{z})\mbox{d}\bm{z}$. The first
estimator averages over the continuous time path of the PDMP

$\frac{1}{t_{n}-t_{1}}\sum_{k=1}^{n-1}\int_{t_{k}}^{t_{k+1}}g(\psi(\bm{z}_{t_{k}},s))ds.$

This requires that the integral of $g(\psi(\bm{z}_{t_{k}},s))$ with respect to
$s$ may be computed. The second, simpler, approach is to “discretise” the PDMP
path. This involves taking an integer $M>0$, defining $s=(t_{n}-t_{1})/M$,
calculating the value of the state at discrete-time points,
$\bm{z}_{t_{1}+s},\ldots,\bm{z}_{t_{1}+Ms}$, and then using the standard Monte
Carlo estimator

$\frac{1}{M}\sum_{j=1}^{M}g(\bm{z}_{t_{1}+js}).$

The points $\bm{z}_{t_{1}+js}$ may be found by evolving the process along its
deterministic dynamics for the appropriate event time interval.

###  2.3 Bouncy particle sampler

The Bouncy Particle Sampler (Bouchard-Côté et al., 2018) takes the velocity
space $\mathcal{V}$ to be either $\mathbb{R}^{p}$ or the unit sphere
$\mathcal{S}^{p-1}$ and targets an invariant distribution
$\nu(\bm{z})=\pi(\bm{\theta})p(\bm{v})$ where $p(\bm{v})$ is either the
standard $p$-dimensional Normal distribution, or a uniform distribution on
$\mathcal{S}^{p-1}$. This sampler evolves a PDMP with linear deterministic
dynamics
$\phi((\bm{\theta}_{t},\bm{v}_{t}),s)=(\bm{\theta}_{t}+s\bm{v}_{t},\bm{v}_{t})$
between events. The canonical event rate is given by

$\lambda^{c}(\bm{z}_{t})=\max\left\\{0,\langle\bm{v}_{t},\nabla_{\bm{\theta}}U(\bm{\theta}_{t})\rangle\right\\}.$

At an event time $t$ generated according to this rate the velocity is updated
as $\bm{v}_{t}=R_{\bm{\theta}_{t-}}(\bm{v}_{t-})$ where

$R_{\bm{\theta}}(\bm{v})=\bm{v}-2\frac{\langle\bm{v},\nabla_{\bm{\theta}}U(\bm{\theta})\rangle}{\|\nabla_{\bm{\theta}}U(\bm{\theta})\|^{2}}\nabla_{\bm{\theta}}U(\bm{\theta}).$

This update can be interpreted as reflecting the velocity off the hyperplane
tangential to the gradient of $U(\bm{\theta})$. Additionally, at random times
$t$ an additional event occurs according to a homogeneous Poisson process with
rate $\lambda^{\text{ref}}>0$ in which the velocity is refreshed, drawing
$\bm{v}_{t}\sim p(\bm{v})$. Refreshment ensures that the Bouncy Particle
Sampler samples the invariant distribution and does not get “stuck” on
contours of the potential (Bouchard-Côté et al., 2018). Formally the Bouncy
Particle Sampler has event rate
$\lambda(\bm{z}_{t})=\lambda^{c}(\bm{z}_{t})+\lambda^{\text{ref}}$ with Markov
transition kernel

$q(\bm{z}|\bm{z}_{t-})=\frac{\lambda^{c}(\bm{z}_{t-})}{\lambda(\bm{z}_{t-})}\delta_{\bm{\theta}_{t-}}(\bm{\theta})\delta_{R_{\bm{\theta}_{t-}}(\bm{v}_{t-})}(\bm{v})+\frac{\lambda^{\text{ref}}}{\lambda(\bm{z}_{t-})}\delta_{\bm{\theta}_{t-}}(\bm{\theta})p(\bm{v}),$

where $\delta_{a}(A)$ denotes the Dirac delta function with point mass at $a$.
The Bouncy Particle Sampler is an example of a global PDMP as all components
of the velocity $\bm{v}$ are changed according to the Markov transition
kernel.

###  2.4 Zig-Zag sampler

The Zig-Zag sampler (Bierkens et al., 2019, 2021a) takes the space of
velocities, $\mathcal{V}$, to be $\\{-1,1\\}^{p}$ with $p(\bm{v})=2^{-p}$
uniform on this space. The Zig-Zag sampler also uses linear dynamics
$\phi((\bm{\theta}_{t},\bm{v}_{t}),s)=(\bm{\theta}_{t}+s\bm{v}_{t},\bm{v}_{t})$
between events, but unlike the BPS only a single element of the velocity
vector is updated at event times. The overall Zig-Zag process can be viewed as
a composition of $p$ local event rates where each event rate has a
corresponding Markov transition kernel which operates on a single component of
the velocity. Specifically, the $i$-th component of the velocity is updated
with local event rate

$\lambda_{i}(\bm{z}_{t})=\max\left\\{0,v_{i}\partial_{\theta_{i}}U(\bm{\theta}_{t})\right\\}$

for $i=1,...,p$. Simulating an event time in this local framework consists of
simulating an event time for each local event rate and then taking the minimum
time as the event time $t$ for the overall process. At an event time $t$
triggered by event rate $\lambda_{i}$ the velocity is updated as
$\bm{v}_{t}=F_{i}(\bm{v}_{t-})$ where

$F_{i}(\bm{v})=(v_{1},\dots,v_{i-1},-v_{i},v_{i+1},\dots,v_{p})^{T}.$

Once the velocity is updated for component $i$ local event times must be re-
simulated for all rates $\lambda_{j}(\bm{z}_{t})$ which are functions of
$\theta_{i}$. This fact can induce massive computational savings when there is
conditional independence between the components of $\bm{\theta}$. To see this,
consider the extreme case where the target is fully independent across
dimensions. In this case each event rate, $\lambda_{i}(\bm{z}_{t})$, will only
depend on the $i$-th components of $\bm{\theta}_{t}$ and $\bm{v}_{t}$. Thus at
each event, it is only the time of the next event for the component whose
velocity changes that needs to be recalculated. The computational cost of
simulation in this case will scale as $O(1)$ per event as opposed to $O(d)$
for the global Bouncy Particle Sampler.

###  2.5 Global and local PDMP methods

As described above, we can divide PDMP samplers into two main classes: global
methods and local methods. These PDMPs differ in how they operate on the
velocity vector when an event is triggered. Global methods are PDMP-based
methods where the transition kernel acts on the entire velocity vector. By
contrast, local methods are defined so that the transition kernel only affects
a subset of the velocity. In this section we will give a general algorithm for
constructing a local PDMP-based sampler (though see also Bouchard-Côté et al.,
2018). Let $\bm{v}^{[S]}$ denote the sub-vector of $\bm{v}$ with elements
indexed by $S\subseteq\\{1,\dots,p\\}$ and $\bm{v}^{[-S]}$ denote the sub-
vector of $\bm{v}$ without the elements indexed in $S$. For some set $S$
define the local kernel as one that satisfies

$q^{[S]}_{\bm{\theta}}(\bm{v}|\bm{v}_{t-})=\delta_{\bm{v}^{[-S]}_{t-}}\left(\bm{v}^{[-S]}\right)q_{\bm{\theta}}^{[S]}(\bm{v}|\bm{v}_{t-}),$

that is, $q^{[S]}_{\bm{\theta}}(\bm{v}|\bm{v}_{t-})$ only updates the sub-
vector $\bm{v}^{[S]}$. Suppose we have a partition of the components of
$\theta_{1},...,\theta_{p}$ given by
$\mathcal{S}=\\{S_{1},S_{2},\dots,S_{F}\\}$. Subset $S_{f}$ of $\bm{v}$ is
updated with local kernel $q^{[S_{f}]}_{\bm{\theta}}(\bm{v}|\bm{v}_{t-})$ and
associated event rate

$\lambda_{f}(\bm{z}_{t})=\max\left\\{0,\langle\bm{v}^{[S_{f}]},\nabla_{\bm{\theta}^{[S_{f}]}}U(\bm{\theta}_{t})\rangle\right\\},$

for $f=1,\dots,F$. Simulating a local PDMP with this structure involves
simulating an event time for each of the $F$ rates and applying the local
transition kernel corresponding to the smallest simulated event time. For the
partition $\mathcal{S}$ define $\mathcal{N}=\\{N_{1},\dots,N_{F}\\}$ with

$N_{f}=\\{m\in\\{1,\dots,F\\}\mid\bm{\theta}^{[S_{m}]}\not\\!\perp\\!\\!\\!\perp\bm{\theta}^{[S_{f}]}\\}$

for $f=1,\dots,F$ where
$\bm{\theta}^{[S_{f}]}\not\\!\perp\\!\\!\\!\perp\bm{\theta}^{[S_{m}]}$ denotes
that the rate of events associated with $S_{m}$ depends on some element of
$\bm{\theta}^{[S_{f}]}$. When an event is triggered for a local rate $f$ the
velocity $\bm{v}^{[S_{f}]}$ is updated and new event times are simulated for
rates which depend on the value of one or more $\theta_{i}$ for $i\in S_{f}$.
This simulation process is summarised in Algorithm 1.

In this framework the Zig-Zag sampler is a local PDMP with partition
$\mathcal{S}=\\{\\{1\\},\\{2\\},\dots,\\{p\\}\\}$ and local kernels
$q^{[S_{f}]}_{\bm{\theta}}(\bm{v}|\bm{v}_{t-})=\delta_{F_{f}(\bm{v}_{t-})}(\bm{v})$.
The local Bouncy Particle Sampler uses kernel
$q^{[S]}_{\bm{\theta}}(\bm{v}|\bm{v}_{t-})=\delta_{\bm{v}_{*}}(\bm{v})$ where
$\bm{v}_{*}^{[-S]}=\bm{v}^{[-S]}_{t-}$ and
$\bm{v}_{*}^{[S]}=R_{\bm{\theta}}^{[S]}(\bm{v}^{[S]}_{t-})$ with

$R_{\bm{\theta}}^{[S]}(\bm{v}^{[S]})=\bm{v}^{[S]}-2\frac{\langle\bm{v}^{[S]},\nabla_{\bm{\theta}^{[S]}}U(\bm{\theta})\rangle}{\|\nabla_{\bm{\theta}^{[S]}}U(\bm{\theta})\|^{2}}\nabla_{\bm{\theta}^{[S]}}U(\bm{\theta}).$

If there is a single factor $\mathcal{S}=\\{S\\}$ containing all indices
$S=\\{1,\dots,p\\}$ we recover the global Bouncy Particle Sampler.

Algorithm 1: Simulating a PDMP with local structure Inputs: $t_{1}=0$,
$\bm{\theta}_{t_{1}},\bm{v}_{t_{1}}$, factorisation
$(\mathcal{S},\mathcal{N})$ with $\mathcal{S}=\\{S_{1},S_{2},\dots,S_{F}\\}$
and $\mathcal{N}=\\{N_{1},N_{2},\dots,N_{F}\\}$. Sample $\tau_{f}$ for
$f=1,...,F$ where $\mathbb{P}(\tau_{f}\geq
t)=\exp\left(-\int_{0}^{t}\lambda_{f}(\bm{\theta}_{0}+\bm{v}_{0}s,\bm{v}_{0})\text{d}s\right)$
For $k=1,\dots,n$ (a) Let $f^{*}=\text{argmin}_{f}(\tau_{f})$ and update
$t_{k+1}=t_{k}+\tau_{f^{*}}$ (b) Update states:
$\bm{\theta}_{t_{k+1}}=\bm{\theta}_{t_{k}}+\tau_{f^{*}}\bm{v}_{t_{k}}$ (c)
Update velocity:
$\bm{v}_{t_{k+1}}=q_{\bm{\theta}_{t_{k+1}}}^{[S_{f^{*}}]}(\bm{v}_{t_{k}}|\bm{v})$
(d) For $f\in N_{f^{*}}$ update times, i.e. resample $\tau_{f}$ where
$\mathbb{P}(\tau_{f}\geq
t)=\exp\left(-\int_{0}^{t}\lambda_{f}(\bm{\theta}_{t_{k+1}}+\bm{v}_{t_{k+1}}s,\bm{v}_{t_{k+1}})\text{d}s\right)$
Outputs: PDMP skeleton
$\\{(t_{k},\bm{\theta}_{t_{k}},\bm{v}_{t_{k}})\\}_{k=1}^{n}$

##  3 Algorithms for event-time simulation

Simulating the first event time, $\tau$, from a Poisson process with event
rate $\lambda(t)$ is equivalent to simulating $u\sim\text{Uniform}[0,1]$ and
solving

$\tau=\text{argmin}_{t}\left\\{-\log(u)=\int_{0}^{t}\lambda(\bm{z}_{s})ds\right\\}.$

If the event rate is sufficiently simple this can be done exactly. For example
if the rate is constant or linear, then there are analytical solutions for
$\tau$ in terms of $u$. If the rate function is convex a solution may be found
using numerical methods (Bouchard-Côté et al., 2018). For more complex
functions, the rate can be simulated by a process known as thinning. This
process makes use of Theorem 1.

######  Theorem 1 (Lewis and Shedler (1979)).

Let $\lambda(t)$ and $\bar{\lambda}(t)$ be continuous functions where
$\lambda(t)\leq\bar{\lambda}(t)$ for all $0\leq t_{\max}$. Let
$\tau_{1},\tau_{2},...,\tau_{n}$ be event times of the Poisson process with
rate $\bar{\lambda}(t)$ over the interval $[0,t_{\max})$. For $i=1,...,n$
retain the event times $\tau_{i}$ with probability
$\lambda(\tau_{i})/\bar{\lambda}(\tau_{i})$ to obtain a set of event times
from a non-homogeneous Poisson process with rate $\lambda$ over the interval
$(0,t_{\max}]$.

The efficiency of simulating via thinning depends on how tightly the rate
$\bar{\lambda}$ upper-bounds $\lambda$ and how costly it is to simulate from
$\bar{\lambda}$. Another key tool for enabling the simulation of event times
is known as superposition and the general process is outlined in Theorem 2.

######  Theorem 2 (Kingman (1992)).

Suppose that $\Lambda_{1},...,\Lambda_{n}$ are a set of independent Poisson
processes with rates $\lambda_{1}(t),...,\lambda_{n}(t)$ with first arrival
times $\tau^{[1]},...,\tau^{[n]}$. The Poisson process with rate
$\lambda(t)=\sum_{i=1}^{n}\lambda_{i}(t)$ may be simulated by returning the
first arrival time as $\tau=\min_{i=1,...,n}\tau^{[i]}$.

Both superposition and thinning provide useful tools in the simulation of
event times from a non-homogeneous Poisson process. However, constructing the
upper-bounds required for thinning is largely an ad hoc and time consuming
process for the practitioner. A common desire for a practitioner is to reuse
simulation knowledge. For example if an event rate can be written as the sum
of several sub-event rates, that can be exactly simulated, ideally this
knowledge should be used to simulate from the sum. We refer to this class of
event rates as additive rates. Specifically, we will refer to a Poisson
process with a rate

$\lambda(\bm{z}_{t})=\max\\{0,f_{1}(t)+f_{2}(t)\\}$

as process with an additive rate. The superposition method from Theorem 2
facilitates additive event-time simulation by thinning using an event rate

$\bar{\lambda}(t)=\max\\{0,f_{1}(t)\\}+\max\\{0,f_{2}(t)\\},$

which satisfies $\lambda(\bm{z}_{t})\leq\bar{\lambda}(t)$. We can simulate the
first event from a process with this rate as $\tau=\min(\tau_{1},\tau_{2})$,
where $\tau_{1}$ is simulated with rate
$\bar{\lambda}_{1}=\max\\{0,f_{1}(t)\\}$ and $\tau_{2}$ with rate
$\bar{\lambda}_{2}=\max\\{0,f_{2}(t)\\}$. This simulated time will be accepted
with probability

$\frac{\max\\{0,f_{1}(\tau)+f_{2}(\tau)\\}}{\max\\{0,f_{1}(\tau))+\max(0,f_{2}(\tau)\\}},$

which will be one when both $f_{1}(\tau)>0$ and $f_{2}(\tau)>0$. This
procedure is useful since it allows the re-use of thinning procedures for
$f_{1}$ and $f_{2}$. Consequently a practitioner can build up to simulating
from a complex rate by combining smaller simpler rates. Thinning via
superposition is often recommended in the literature; notable examples include
Bouchard-Côté et al. (2018) who illustrate this approach when simulating from
an exponential family and Wu and Robert (2020) who discuss the approach
generally. The approach was also recommended broadly by Sen et al. (2019)
where it was suggested for practical implementation of the Zig-Zag where
$f_{1}$ and $f_{2}$ correspond to the terms from the likelihood and prior
respectively. Once simulation for a choice of prior or likelihood are
individually established this knowledge can be reused in future modelling.

##  4 Concave-convex adaptive thinning

Our general proposal for simulating events in a PDMP is based on concave-
convex adaptive thinning. As the dynamics between events are deterministic,
conditional on the current state of the PDMP, we can re-write the rate of the
next event as a function of time. With slight abuse of notation we will denote
this rate as $\lambda(t)$ with the argument of the function making it clear
whether we are viewing it as a function of the state, or as here, as a
function of time. Suppose $\lambda(t)=\max\\{0,f(t)\\}$ where $f(t)$ can be
decomposed as:

$\displaystyle f(t)=f_{u}(t)+f_{n}(t)$ (1)

on a finite interval $t\in[0,\tau_{\max})$ where $f_{u}(t)$ is a convex
function and $f_{n}(t)$ is a concave function in time. (We discuss below
general conditions where such a decomposition is possible.) The problem of
upper-bounding $f(t)$ is recast to finding upper-bounding piecewise linear
functions $\ell_{u}(t)$ and $\ell_{n}(t)$ such that $f_{u}(t)\leq\ell_{u}(t)$
and $f_{n}(t)\leq\ell_{n}(t)$. We may then apply the bound $f(t)\leq\ell(t)$
where $\ell(t)=\ell_{n}(t)+\ell_{u}(t)$. Since $\ell(t)$ is the sum of
piecewise linear functions it will also be a piecewise linear function and
direct simulation from a Poisson process with rate $\ell(t)$ is possible (see
Appendix A).
Concave-convex adaptive thinning proceeds by constructing a piecewise linear
function $\ell(t)$ over a set of $m$ abscissae
$t_{0}=0<t_{1}<\dots<t_{m}=\tau_{\max}$ at which the evaluation of
$f_{u}(t_{i}),f_{n}(t_{i})$ and derivative $f_{n}^{\prime}(t_{i})$ are known.
An event time $\tau$ is simulated from the Poisson process with rate
$\bar{\lambda}(t)=\max\\{0,\ell(t)\\}$ and is accepted with probability

$\frac{\max\\{0,f_{u}(\tau)+f_{n}(\tau)\\}}{\bar{\lambda}(\tau)}.$

If the event is rejected, we can reuse the information from the evaluation of
$f_{u}(\tau)$ and $f_{n}(\tau)$ with an additional evaluation of
$f^{\prime}_{n}(\tau)$ to refine the simulation on the abscissae
$t_{0}=\tau<t_{k}<\dots<t_{m}=\tau_{\max}$ where $t_{k-1}=\max\\{t_{i}\mid
t_{i}<\tau\\}$. If an event does not occur on the range $t\in[0,\tau_{\max})$
the PDMP process is evolved by $\tau_{\max}$ and the thinning process is
repeated. We give the general construction of the piecewise linear function
$\ell(t)$ for $t\in[t_{1},t_{2})$ and note this construction may be applied
iteratively over the abscissae. This construction is also depicted visually in
Figure 1.

Figure 1: Upper bounds for a rate function based on concave (left) and convex
(right) information on three abscissae. Concave bounds are formed using linear
segments from the gradient of $f_{n}(t)$. Convex bounds are constructed using
linear segments connecting evaluations of $f_{u}(t)$.

A convex function $f_{u}(t)$, is by definition a function that can be upper-
bounded by the line segment connecting any two function evaluations. So, for
any $t\in[t_{1},t_{2})$ we have the upper-bound $f_{u}(t)\leq\ell_{u}(t)$
where

$\ell_{u}(t)=f_{u}(t_{1})+\frac{f_{u}(t_{2})-f_{u}(t_{1})}{t_{2}-t_{1}}t.$

For a concave function $f_{n}(t)$ with derivative $f_{n}^{\prime}(t)$ the line
segment corresponding to the tangent of $f_{n}(t)$ will upper-bound the
function. Thus $f_{n}(t)$ will be upper-bounded by

$\ell_{n}(t)=\min\left\\{f_{n}(t_{1})+f_{n}^{\prime}(t_{1})t,f_{n}(t_{2})+f_{n}^{\prime}(t_{2})(t_{2}-t)\right\\}.$

for $t\in[t_{1},t_{2})$. This minimum will switch at the point of intersection
between the lines $f_{n}(t_{1})+f_{n}^{\prime}(t_{1})t$ and
$f_{n}(t_{2})+f_{n}^{\prime}(t_{2})(t_{2}-t)$. It is simple to find this
intersection point (Görür and Teh, 2011),

$t^{*}=\frac{f_{n}(t_{2})-f_{n}^{\prime}(t_{2})t_{2}-f_{n}(t_{1})+f_{n}^{\prime}(t_{1})t_{1}}{f_{n}^{\prime}(t_{1})-f_{n}^{\prime}(t_{2})}$

which will be a point on the interval $[t_{1},t_{2}]$ provided the derivatives
$f_{n}^{\prime}(t_{1})$ and $f_{n}^{\prime}(t_{2})$ are not equal. If
$f_{n}^{\prime}(t_{1})=f_{n}^{\prime}(t_{2})$, the linear segment will not
change over the interval and we take $t^{*}=t_{2}$. So we take

$\ell_{n}(t)=\begin{cases}f_{n}(t_{1})+f_{n}^{\prime}(t_{1})t&t\in[t_{1},t^{*})\\\
f_{n}(t_{2})+f_{n}^{\prime}(t_{2})(t_{2}-t)&t\in[t^{*},t_{2})\end{cases},$

and combining these bounds we have $\ell(t)=\ell_{u}(t)+\ell_{n}(t)$ which is
a piecewise linear function upper-bounding $f(t)\leq\ell(t)$ for
$t\in[t_{1},t_{2})$.

###  4.1 Concave-convex decompositions

The proposition below gives simple conditions for the class of non-homogeneous
Poisson processes that admit thinning using the concave-convex approach.

######  Proposition 1.

If a Poisson process with rate function $\lambda(t)$ can be written as
${\lambda(t)=\max\\{0,f(t)\\}}$ where $f(t)$ is continuous except possibly at
a finite number of known points, $0<t^{\star}_{1}<t^{\star}_{2}<\cdots
t_{r}^{\star}<\infty$ , then the process admits thinning using the concave-
convex adaptive thinning approach.

###### Proof.

Consider $f(t)$ on the interval $[0,\tau_{\max})$ where if $t_{1}^{\star}$
exists, $\tau_{\max}=t_{1}^{\star}$ otherwise $\tau_{\max}$ is some arbitrary
value $\tau_{\max}>0$. On this interval $f(t)$ is continuous and by the Stone-
Weierstrass theorem (Stone, 1948), for any $\epsilon>0$ there exists a
polynomial $g(t)$ on the interval $[0,\tau_{\max}]$ with $\|f-g\|<\epsilon$.
The function $g$ is a polynomial so it admits the concave-convex decomposition

$g(t)=\sum_{\\{i:a_{i}>0\\}}a_{i}t^{i}+\sum_{\\{i:a_{i}<0\\}}a_{i}t^{i}$

with convex function $g_{u}(t)=\sum_{\\{i:a_{i}>0\\}}a_{i}t^{i}$ and concave
function $g_{n}(t)=\sum_{\\{i:a_{i}<0\\}}a_{i}t^{i}$. The Process with rate
$\lambda(t)$ admits thinning on the interval $[0,\tau_{\max})$ using the
Poisson process with rate $\hat{\lambda}(t)=\max\\{0,g(t)+\epsilon\\}$ which
has a convex-concave decomposition. If an event does not occur on the interval
$[0,\tau_{\max})$ then the process is considered on the next continuous
interval. ∎

The proof of Proposition 1 relies on finding a polynomial that can upper-bound
the process over a finite interval. If a bound can be given for higher order
derivatives of $f(t)$ then there are two main approaches for finding an upper-
bounding polynomial. Using Taylor’s expansion of $f(t)$ about zero

$f(t)=f(0)+f^{\prime}(0)t+\frac{f^{\prime\prime}(0)}{2!}t^{2}+\cdots+\frac{f^{(k)}(t)}{k!}t^{k}+\int_{0}^{t}f^{(k+1)}(s)\frac{t^{k}}{k!}ds.$

If there is a known constant $M$ such that $|f^{(k+1)}(t)|\leq M$ we can
employ thinning based on the polynomial bound
$g(t)=f(0)+f^{\prime}(0)t+\frac{f^{\prime\prime}(0)}{2!}t^{2}+\cdots+\frac{f^{(k)}(t)}{k!}t^{k}+\frac{M}{(k+1)!}t^{(k+1)}$.
Alternatively polynomial upper-bounds can be constructed based on
interpolation where the error is controlled. For example using Lagrange
polynomial interpolation on $[0,\tau_{\max}]$ with $k$ evaluations of the
function has error is bounded by $\tau_{\max}^{k+1}\frac{M}{(k+1)!}$. Adding
this bound as a constant yields an upper-bounding polynomial without requiring
evaluation of the derivatives of $f(t)$ required for the Taylor series
expansion.

###  4.2 Concave-convex adaptive thinning and PDMP-based samplers

We will refer to the process of using concave-convex adaptive thinning for
simulating a PDMP-based sampler as CC-PDMP sampling. In this section we
illustrate how CC-PDMP sampling is applied when the target distribution
corresponds to a Bayesian posterior:

$\pi(\bm{\theta})\propto\exp(-U^{(\ell)}(\bm{\theta})-U^{(p)}(\bm{\theta}))$

where $U^{(\ell)}(\bm{\theta})$ is the potential of the likelihood and
$U^{(p)}(\bm{\theta})=-\log p_{0}(\bm{\theta})$ is the potential for the prior
distribution $p_{0}$. Consider the $k$th rate for a fully local PDMP such as
the Zig-Zag sampler

$\lambda_{k}(t)=\max\left\\{0,v_{k}\partial_{\theta_{k}}U^{(\ell)}(\bm{\theta}+t\bm{v})+v_{k}\partial_{\theta_{k}}U^{(p)}(\bm{\theta}+t\bm{v})\right\\},$

here we could apply concave-convex thinning with
$f_{k}(t)=v_{k}\partial_{\theta_{k}}U^{(\ell)}(\bm{\theta}+t\bm{v})+v_{k}\partial_{\theta_{k}}U^{(p)}(\bm{\theta}+t\bm{v})$.
For Bayesian inference problems the function $f_{k}$ is the sum
$f_{k}(t)=f_{k}^{1}(t)+f_{k}^{2}(t)$ of likelihood and prior part where
$f_{k}^{1}(t)=v_{k}\partial_{\theta_{k}}U^{(\ell)}(\bm{\theta}+t\bm{v})$ and
$f_{k}^{2}(t)=v_{k}\partial_{\theta_{k}}U^{(p)}(\bm{\theta}+t\bm{v})$. If
$f_{k}^{1}$ and $f_{k}^{2}$ have concave-convex decompositions then $f_{k}$
has a concave-convex decomposition. This allows for reuse of bounding
knowledge, priors and likelihoods can be swapped provided their decompositions
are known. Moreover, to implement this fully local method the rates
$\lambda_{k}(t)=\max\\{0,f_{k}(t)\\}$ must have a concave-convex decomposition
for $k=1,...,p.$ Thus we can reuse this information to find a decomposition
for a global PDMP-based sampler with rate

$\lambda(t)=\max\\{0,\langle\bm{v},\nabla_{\bm{\theta}}U(\bm{\theta}+t\bm{v})\rangle\\}=\max\left\\{0,\sum_{k=1}^{p}f_{k}(t)\right\\}.$

Terms can be trivially regrouped to facilitate simulation from any choice of
local PDMP structure as described in Section 2.5.

##  5 Tuning parameters

###  5.1 Choice of abscissae

A key user choice in CC-PDMP sampling is that of the position and number of
abscissae on the interval $[0,\tau_{\max})$. Suppose that CC-PDMP simulation
is implemented with abscissae $t_{0}=0<t_{1}<\cdots<t_{n}=\tau_{\max}$ if an
event does not occur on the interval $[0,\tau_{\max})$ the PDMP will be
evolved by $\tau_{\max}$ and the thinning process will repeat where the
function evaluations at $\tau_{\max}$ can be reused at $t_{0}=0$. Thus one
simulation with $n$ abscissae would be computationally equivalent to $n$
evolutions of the CC-PDMP approach using abscissae with only two points
$t_{0}=0$ and $t_{1}=\tau_{\max}$. This encourages choosing a minimal number
of abscissae and more carefully choosing the length of the interval.

There are two main issues that can occur when tuning the parameter
$\tau_{\max}$. If $\tau_{\max}$ is too small, then the proposal frequently lie
outside the interval $[0,\tau_{\max})$ and simulating an event will requiring
many iterations of bounding the event rate. Whilst if $\tau_{\max}$ is too
large then the bound on the rate may be be poor, leading to many events that
are rejected. In this article, we refer to an iteration of the CC-PDMP that
does not generate an event time (i.e. a rejected proposal event or not
generating on the interval) as a shadow event. The total number of iterations
will be the sum of the number of events and shadow events. Efficiency is
measured as the proportion of iterations that are events not shadow events,
calculated as the ratio of the number of events to the number of iterations.

Adapting the value of $\tau_{\max}$ will not change the sampling dynamics.
Consequently, unlike in adaptive MCMC, this parameter may be adapted
throughout the course of simulating the PDMP. Ideally this parameter should be
a little larger than the average event time so that events are proposed on the
interval. A simple automatic approach for choosing this parameter is to set
$\tau_{\max}$ equal to the $q$-th percentile of previous simulated event
times, and update $\tau_{\max}$ every 100 iterations of the sampler. We found
that setting $q=80$th percentile worked well for automatically selecting
$\tau_{\max}$.

We investigate the dependency on this tuning parameter by implementing the
Zig-Zag sampler on the Banana distribution. The 2-dimensional Banana
distribution has potential:

$U(\bm{\theta})=(\theta_{1}-1)^{2}+\kappa(\theta_{2}-\theta_{1}^{2})^{2}$

where $\kappa$ controls how much mass concentrates around the region
$\theta_{2}\approx\theta_{1}^{2}$.

Since the potential is polynomial in both parameters it is trivial to see that
the Zig-Zag rate functions will also be polynomial. Specifically the rates
will be polynomial functions of time $t$ of degree 3 and 2 for parameters
$\theta_{1}$ and $\theta_{2}$. Further implementation details may be found in
the supplementary material and associated GitHub.

Figure 2 shows the efficiency of the Zig-Zag sampler on the Banana
distribution with $\kappa=1$ when changing the interval length $\tau_{\max}$
from 0.1 to 4 with two abscissae. The Zig-Zag algorithm was simulated for
10000 events at each value of $\tau_{\max}$. The effect of selecting
$\tau_{\max}$ too small or too large is clearly seen and the red line shows
the performance using the adaptive approach. While it is possible to get low
efficiency if the parameter $\tau_{\max}$ is chosen poorly, we have found
either using the adaptive approach or a short pilot run of the sampler to tune
$\tau_{\max}$ gives reliable performance.

Figure 2: Efficiency, defined as the proportion of events that are not shadow
events, of Zig-Zag using the CC-PDMP thinning with varying choices for
$\tau_{\max}$ on the Banana target. The red line shows the performance of an
automatic choice for this parameter.

###  5.2 Choice of decomposition

The choice of concave-convex decomposition is not unique (Görür and Teh,
2011). Arbitrary convex functions can be added to $f_{u}(t)$ and subtracted
from $f_{n}(t)$ to give a new valid decomposition. A concave-convex
decomposition is minimal if there is no non-affine convex function that can be
added to $f_{n}(t)$ and subtracted from $f_{u}(t)$ while preserving the
convexity of the decomposition (Hartman, 1959). Consider the function

$f(t)=-t^{3}+3t^{2}-3t+3,$

where $f^{\prime\prime}(t)=-6t+6$. It is simple to see that the function is
convex for $t<1$ and concave for $t>1$. A minimal decomposition is given by
the piecewise functions

$f_{u}^{1}(t)=\begin{cases}-t^{3}+3t^{2}-3t+3&t\leq 1\\\
0&t>1\end{cases}\text{ and }f_{n}^{1}(t)=\begin{cases}0&t\leq 1\\\
-t^{3}+3t^{2}-3t+3&t>1\end{cases}.$

Görür and Teh (2011) give a general construction for a minimal concave-convex
decomposition. However, this construction relies on finding all points of
inflection, which are points where the function changes convexity.
Alternatively, Proposition 1 gives a simple decomposition
$f_{u}^{2}(t)=3t^{2}+3$ and $f_{n}^{2}(t)=-t^{3}-3t$. While the decomposition
from Proposition 1 is not minimal, it does not require finding all points of
inflection. These two decompositions are shown in Figure 3 using abscissae at
0 and 1. The optimal decomposition gives a tighter bound, here reducing the
bounding rate by approximately 0.33 on average across the interval. However,
using this bound comes with additional computation cost of finding inflection
points. This has the potential to reduce the overall efficiency of the method.
Our experience is that the efficiency gains for using the optimal polynomial
are generally not sufficient for it to be beneficial.

Figure 3: The upper-bounds resulting from two different concave-convex
decompositions of the function $f(t)=-t^{3}+3t^{2}-3t+3$, based on abscissae
at 0 and 1. The top row corresponds to the decomposition from Proposition 1
and the bottom row corresponds to an optimal decomposition which recognises
$f(t)$ as a convex function on $[0,1)$. The columns show the functions
$f_{u}$, $f_{n}$ and $f$ as well as their piece-wise linear upper-bounds.

##  6 Experiments

We now present empirical evaluation of CC-PDMP and comparison with other
approaches to simulate PDMPs. Our experiments were implemented using the R
package ccpdmp available at https://github.com/matt-sutton/ccpdmp.

The package enables you to simulate a PDMP provided one specifies the concave
and convex decomposition of the rate function. If the rate function is
polynomial (or bounded by a polynomial) the practitioner may parse this
instead and the concave-convex decomposition will be handled internally. The
package contains some basic example use cases and code to reproduce the
experiments. Code snippets are given in the additional material.

###  6.1 Application in generalised linear models

Generalised linear models (GLMs) provide a rich class of models that are
frequently used in Bayesian inference. Let
$\\{(y_{i},\bm{x}_{i})\\}_{i=1}^{n}$ be the observed data, where $y_{i}$ is an
observed response and $\bm{x}_{i}\in\mathbb{R}^{p}$ is a vector of associated
covariates for $i=1,\dots,n$. The expected value of $y_{i}$ is modelled by
$g^{-1}(\bm{x}_{i}^{T}\bm{\theta})$ where
$g^{-1}:\mathbb{R}\rightarrow\mathbb{R}$ is the inverse link function. The
potential of the likelihood has the form

$U^{(\ell)}(\bm{\theta})=\sum_{i=1}^{n}\phi(\bm{x}_{i}^{T}\bm{\theta},y_{i})$

where $\phi:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is the GLM mapping function
$\phi(a,y)=\log p(y|g^{-1}(a))$ returning the log likelihood for observation
$y$ given $a$. The partial derivative with respect to $\theta_{k}$ has the
form

$\partial_{\theta_{k}}U^{(\ell)}(\bm{\theta})=\sum_{i=1}^{n}\phi^{\prime}(\bm{x}_{i}^{T}\bm{\theta},y_{i})x_{ik}$

where $\phi^{\prime}(a,y)=\frac{\partial}{\partial a}\phi(a,y)$ and higher
order derivatives are defined similarly. Over the time interval $t\in[0,\tau]$
let $f_{k}(t)=v_{k}\partial_{\theta_{k}}U^{(\ell)}(\bm{\theta}+t\bm{v})$ which
is

$f_{k}(t)=v_{k}\sum_{i=1}^{n}\phi^{\prime}(a_{i}(t),y_{i})x_{ik}$

where $a_{i}(t)=\bm{x}_{i}^{T}(\bm{\theta}+t\bm{v})$. We can use this to
define local rates as $\lambda_{k}(t)=\max\\{0,f_{k}(t)\\}$ for $k=1,...,p$.
For GLMs, repeated application of the chain rule yields:

$f_{k}^{(j)}(t)=v_{k}\sum_{i=1}^{n}\phi^{(j+1)}(a_{i}(t),y_{i})\left(\frac{\partial
a_{i}}{\partial t}\right)^{j}x_{ik}$

where $\frac{\partial a_{i}}{\partial t}=\bm{x}_{i}^{T}\bm{v}$ and
$f_{k}^{(j)}$ denotes the $j$th derivative of $f_{k}(t)$. When there are
bounds on $\phi^{(j)}$ we can use the upper-bounding Taylor polynomial. In
Appendix B we provide bounds for several modelling choices.

Here we consider a logistic regression GLM which has the mapping function
$\phi(a,y)=\log(1+\exp(a))-ya$ where $y\in\\{0,1\\}$. We look at the
efficiency of CC-PMDP thinning for a three dimensional logistic regression
problem with $n=200$ observations. The covariates were generated from a
multivariate normal, $\bm{x}_{i}\sim\mathcal{N}(\bm{0},V^{-1})$ for
$i=1,..,200$, with mean zero, precision matrix

$V=\left(\begin{matrix}1&\rho&0&0&0\\\ \rho&1&0&0&0\\\ 0&0&1&0&0\\\
0&0&0&1&0\\\ 0&0&0&0&1\\\ \end{matrix}\right)$

and data generated using $\theta=(-1.25,0.5,-0.4,-0.4,-0.4)^{T}$. We take a
Gaussian prior $\theta_{j}\sim\mathcal{N}(0,1)$ for $j=1,\dots,5.$ As the
correlation is increased the thinning becomes more challenging. We investigate
the performance of CC-PDMP for increasing polynomial bounds with increasing
$\rho$.

| Correlation $(\rho)$
---|---
Polynomial Order | 0.00 | 0.25 | 0.50 | 0.65 | 0.75 | 0.85 | 0.95
1 | 0.53 | 0.50 | 0.45 | 0.39 | 0.34 | 0.27 | 0.15
2 | 0.80 | 0.80 | 0.79 | 0.78 | 0.76 | 0.71 | 0.46
3 | 0.82 | 0.82 | 0.82 | 0.82 | 0.81 | 0.79 | 0.62
Table 1: Efficiency of thinning in Zig-Zag sampling of logistic regression for
increasing order of Taylor series polynomial thinning bound. Efficiency is
measured as the average proportion of iterations that are events not shadow
events over 20 repetitions of the sampler.

To date only linear bounds for logistic regression have been used for thinning
(Bierkens et al., 2018; Bouchard-Côté et al., 2018). These are based on the
bound $|\phi^{\prime\prime}(a,y)|\leq 1/4$. Table 1 shows the efficiency for
thinning using polynomials of order 1-3 (using bounds
$|\phi^{\prime\prime}(a,y)|\leq 1/4$, $|\phi^{\prime\prime\prime}(a,y)|\leq
1/(6\sqrt{3})$ and $|\phi^{(4)}(a,y)|\leq 1/8$). The linear 1st order bound
matches the bound used in Bierkens et al. (2018) for logistic regression using
the Zig-Zag sampler. Table 1 shows that higher order polynomial thinning
facilitated through the CC-PDMP allows more efficient thinning. For the linear
bounds (Polynomial order 1) we used a fixed $\tau_{\max}=1$ as linear
proposals are exactly simulated using the concave-convex approach. For other
polynomial orders we used the adaptive procedure described in the tuning
section to select $\tau_{\max}$.

###  6.2 Comparison to thinning via superposition

In this example we demonstrate the advantages of CC-PDMP when the event rate
is additive. In particular, we compare thinning using the CC-PDMP approach
with thinning using superposition as outlined in Section 3. Let
$\\{y_{i}\\}_{i=1}^{n}$ be the observed data with the following model

$\displaystyle y_{i}\mid\theta_{i}$
$\displaystyle\sim\text{Poisson}(\exp(\theta_{i})),$ $\displaystyle\theta_{i}$
$\displaystyle\sim\mathcal{N}(0,1)$

independently for $i=1...,n$. The derivative of the potential is
$\partial_{\theta_{k}}U(\bm{\theta})=\theta_{k}-y_{k}+\exp(\theta_{k})$ and

$f_{k}(t)=v_{k}\partial_{\theta_{k}}U(\bm{\theta}+t\bm{v})=v_{k}(\theta_{k}+v_{k}t)-y_{k}v_{k}+v_{k}\exp(\theta_{k}+v_{k}t).$

Which has the concave-convex decomposition $f_{u}(t)=f_{k}(t)$ when $v_{k}>0$
and $f_{u}(t)=v_{k}(\theta_{k}+v_{k}t)-y_{k}v_{k}$,
$f_{n}(t)=v_{k}\exp(\theta_{k}+v_{k}t)$ when $v_{k}<0$. The global Bouncy
Particle Sampler rate can be defined as

$\lambda(t)=\max\left\\{0,\sum_{k=1}^{n}v_{k}\partial_{\theta_{k}}U(\bm{\theta}+t\bm{v})\right\\}=\max\left\\{0,\sum_{k=1}^{n}f_{k}(t)\right\\},$

which has concave-convex decomposition defined by the decompositions of the
individual $f_{k}$ functions. In contrast, thinning via superposition involves
simulating an event time for the linear and exponential terms of the rate
individually. The proposed event time is the minimum of all simulated times.
This approach to simulating event times for exponential likelihood and
Poisson-Gaussian Markov random fields has been previously used in implementing
PDMPs Bouchard-Côté et al. (2018); Vanetti et al. (2017).

Figure 4: Efficiency of thinning for the global BPS using thinning via super-
position and CC-PDMP thinning with increasing dimension. Efficiency of the
thinning (left) and total computation time (right) is averaged over 20
repeated runs. The Bouncy Particle Sampler samplers were simulated for 1000
event times on each run.

Figure 4 compares the thinning and overall computational efficiency for the
Bouncy Particle Sampler applied to the Poisson likelihood problem. In the CC-
PDMP approach we used the adaptive update for $\tau_{\max}$ as described in
the tuning section. As the dimension increases we see the proportion of
iterations that are events using superposition drops quickly, consequently the
overall computation time for this approach scales poorly. The poor performance
of the superposition approach with increasing dimension is the result of the
large number of exponential terms, $v_{k}\exp(\theta_{k}+v_{k}t)$, in the
event rate. The thinning acceptance rate for this proposal is

$\frac{\max\left\\{0,\sum_{k=1}^{n}v_{k}(\theta_{k}+v_{k}t-y_{k}+\exp(\theta_{k}+v_{k}t))\right\\}}{\max\left\\{0,\sum_{k=1}^{n}(v_{k}(\theta_{k}+v_{k}t)-y_{k}v_{k})\right\\}+\sum_{k=1}^{n}\max\left\\{0,v_{k}\exp(\theta_{k}+v_{k}t)\right\\}}$

When $v_{k}<0$ these exponential terms contribute zero to the denominator of
the superposition approach. In comparison the CC-PDMP approach uses a linear
upper-bound on these exponential terms which can be negative and reduce the
denominator term, leading to more efficient thinning proposals. This is seen
in Figure 4 where the thinning efficiency remains roughly constant with
increasing dimension.

###  6.3 Local Methods

Local PDMP methods take advantage of the conditional independence between
parameters. Consider the following extension to the previous Poisson example

$\displaystyle y_{i}\mid\theta_{i}$
$\displaystyle\sim\text{Poisson}(\exp(\theta_{i}))\qquad\text{independently
for }i=1,\dots,n$ $\displaystyle\theta_{1}$
$\displaystyle\sim\mathcal{N}(0,1/(1-\rho^{2}))$
$\displaystyle\theta_{i}\mid\theta_{i-1}$
$\displaystyle\sim\mathcal{N}(\rho\theta_{i-1},1),\qquad\text{for
}i=2,\dots,n$

where we fix $\rho=0.5$. The prior on $\bm{\theta}$ corresponds to an AR(1)
process. The partial derivative is

$\partial_{\theta_{k}}U(\bm{\theta})=\begin{cases}(1+\rho^{2})\theta_{k}-\rho\theta_{k+1}-y_{k}+\exp(\theta_{k})&k=1\\\
(1+\rho^{2})\theta_{k}-\rho\theta_{k-1}-\rho\theta_{k-1}-y_{k}+\exp(\theta_{k})&1<k<n\\\
\theta_{k}-\rho\theta_{k-1}-y_{k}+\exp(\theta_{k})&k=n\end{cases}$

which has the same linear and exponential form as the rate in Section 5.2 so
we can use an analogous concave-convex decomposition for this rate. In this
section we consider local PDMP implementations. The CC-PDMP implementation
facilitates simple construction of thinning bounds for local PDMP
factorisations. For the partition $\mathcal{S}=\\{S_{1},\dots,S_{F}\\}$ the
local rate for factor $f$ is

$\lambda_{f}(\bm{z}_{t})=\max\left\\{0,\sum_{k\in
S_{f}}v_{k}\partial_{\theta_{k}}U(\bm{\theta}_{t})\right\\}=\max\left\\{0,\sum_{k\in
S_{f}}f_{k}(t)\right\\},$

for $f=1,\dots,F$. In this section we consider a number of local
decompositions of the form $\mathcal{S}^{(j)}=\\{S_{1},\dots,S_{F}\\}$ where
$S_{1}=\\{1,2,...,j\\}$, $S_{2}=\\{j+1,j+2,..,2j\\}$,
$\dots,S_{F}=\\{p-j,\dots,p-1,p\\}$. If $p\neq Fj$ then the final factor will
have fewer elements than the previous factors. For these decompositions the
conditional independence between parameters gives the following neighbours
$\mathcal{N}=\\{N_{1},\dots,N_{F}\\}$ where

$N_{f}=\begin{cases}\\{f,f+1\\}&f=1\\\ \\{f-1,f,f+1\\}&1<f<F\\\
\\{f-1,f\\}&f=F\end{cases}.$

This local decomposition offers computational advantages since updating a
single rate requires resimulation of three additional rates regardless of the
dimension of the problem. However, this computational efficiency comes at a
loss of statistical efficiency. In particular, a global PDMP will be able to
move further in stochastic time than local methods before requiring a new
event simulation as the global rate function will lower bound the local rate.
The overall efficiency of the PDMP method should therefore be considered in
terms of both its computational and statistical efficiency.

In Figure 5 we investigate the scaling performance of local methods in
comparison with alternative state-of-the-art MCMC methods. In particular, we
compare local PDMP methods with well-tuned implementations of Metropolis
adjusted Langevin algorithm (MALA) and Hamiltonian Monte Carlo (HMC). MALA is
tuned to obtain acceptance rate approximately equal to 0.5 and HMC tuned to
scale with acceptance probability approximately equal to 0.6. For MALA this
involves scaling the variance in the proposal and for HMC this involves
adjusting the number of leapfrog steps per iteration. These implementations
are in line with theoretical scaling results when the parameters of the
distribution are all independent.

Figure 5: Log-log breakdown of computational, statistical and overall
empirical efficiency scaling with dimension. The ESS values are calculated
with respect to the first coordinate $\theta_{1}$ using the coda R package on
a discretised trajectory of the PDMPs. Plotted are the average rates and error
bars of all methods calculated from 50 repeated runs of the methods. The
legend shows the slope and intercept fitted for each method giving empirical
evidence for scaling rates.

As expected, the computational efficiency for the local methods remains
(approximately) constant with increasing dimension. Based on the fitted models
shown alongside the legends, it appears that computational cost for both MALA
and Bouncy Particle Sampler scales approximately as $O(d^{-1})$. The
computational cost for HMC scales at a rate greater than $O(d^{-1})$ due to
the increase in the number of leapfrog iterations. The statistical efficiency
for the local methods improves with larger local factors and is highest for
the Global Bouncy Particle Sampler. The statistical efficiency for MALA drops
at a rate roughly equal to $O(d^{-1/3})$. Overall efficiency can be seen as
the sum of the log computational and log statistical efficiencies. The Zig-Zag
and local Bouncy Particle Sampler methods attain the best overall efficiency
scaling rates around $O(d)$ or better. It is clear that there is a trade-off
between the statistical efficiency of larger local factors and the increased
computational cost. This decomposition can be thought of as an additional
choice for PDMP implementation that can easily be tuned using the CC-PDMP
approach for thinning. Despite a poorer scaling with dimension, HMC remains
competitive with local methods to a very high dimension.

##  7 Conclusion

PDMP-based samplers have shown advantages over traditional MCMC samplers, but
their use has been limited by the perceived difficulty in simulating the PDMP
dynamics. We have introduced CC-PDMP as a general approach that can simulate
the PDMP provided we can specify a concave-convex decomposition of the event
rate. This method has broad applicability, enables simple implementation of
local PDMP methods, and empirically outperforms alternative simulation
approaches. Additional generalisations of the CC-PDMP approach are also
possible by making use of other ideas for adaptive rejection sampling. For
example, to bound the concave term, $f_{n}(t)$, of the rate function we
require that the derivative $f_{n}^{\prime}(t)$ is known on the interval
$t\in[0,\tau_{\max})$. However, we could instead use a looser bound that does
not require this derivative information based on the adaptive rejection bounds
proposed in Gilks (1992). If the rate function is particularly computationally
intensive, a lower bound based on the abscissae where $f_{u},f_{n}$ and
$f^{\prime}_{n}$ and $f_{u}^{\prime}$ are evaluated may be used to perform
early rejection in the thinning. Additional generalisations may also be found
in the difference of convex functions programming literature (Le Thi and Pham
Dinh, 2018). A useful approach to finding a decomposition was to construct a
polynomial approximation that upper-bounds the rate. An approximate version of
our algorithm could also be implemented where this polynomial approximation is
estimated via interpolation. This would allow thinning without having to find
the concave-convex decomposition. However, the overall algorithm would be
biased unless error in the interpolation is added to the polynomial rate.

Finally while our CC-PDMP approach has been illustrated only on PDMPs with
linear dynamics the approach could be used more generally on PDMPs with
nonlinear dynamics, such as the Boomerang sampler of Bierkens et al. (2020).
The only requirement is that the rate function can be bounded by a function
with a concave-convex decomposition.

## References

  * Bierkens and Roberts (2017) Bierkens, J. and Roberts, G. (2017). A piecewise deterministic scaling limit of lifted Metropolis–Hastings in the Curie–Weiss model. The Annals of Applied Probability 27(2), 846–882.
  * Bierkens et al. (2018) Bierkens, J., Bouchard-Côté, A., Doucet, A., Duncan, A. B., Fearnhead, P., Lienart, T., Roberts, G. and Vollmer, S. J. (2018). Piecewise deterministic Markov processes for scalable Monte Carlo on restricted domains. Statistics & Probability Letters 136, 148–154.
  * Bierkens et al. (2019) Bierkens, J., Fearnhead, P. and Roberts, G. (2019). The zig-zag process and super-efficient sampling for Bayesian analysis of big data. The Annals of Statistics 47(3), 1288–1320.
  * Bierkens et al. (2020) Bierkens, J., Grazzi, S., Kamatani, K. and Roberts, G. (2020). The boomerang sampler. In: International Conference on Machine Learning, PMLR, 908–918.
  * Bierkens et al. (2021a) Bierkens, J., Grazzi, S., van der Meulen, F. and Schauer, M. (2021a). A piecewise deterministic Monte Carlo method for diffusion bridges. Statistics and Computing 31(3).
  * Bierkens et al. (2021b) Bierkens, J., Grazzi, S., van der Meulen, F. and Schauer, M. (2021b). Sticky PDMP samplers for sparse and local inference problems. ArXiv.2103.08478.
  * Bouchard-Côté et al. (2018) Bouchard-Côté, A., Vollmer, S. J. and Doucet, A. (2018). The bouncy particle sampler: A nonreversible rejection-free Markov chain Monte Carlo method. Journal of the American Statistical Association 113(522), 855–867.
  * Chevallier et al. (2020) Chevallier, A., Fearnhead, P. and Sutton, M. (2020). Reversible jump PDMP samplers for variable selection. ArXiv.2010.11771.
  * Chevallier et al. (2021) Chevallier, A., Power, S., Wang, A. Q. and Fearnhead, P. (2021). Pdmp monte carlo methods for piecewise-smooth densities. ArXiv.2111.05859.
  * Cotter et al. (2020) Cotter, S., House, T. and Pagani, F. (2020). The NuZZ: Numerical zigzag sampling for general models. ArXiv.2003.03636.
  * Davis (1993) Davis, M. (1993). Markov Models and Optimization. Chapman Hall.
  * Diaconis et al. (2000) Diaconis, P., Holmes, S. and Neal, R. M. (2000). Analysis of a nonreversible Markov chain sampler. Annals of Applied Probability 10, 726–752.
  * Fearnhead et al. (2018) Fearnhead, P., Bierkens, J., Pollock, M. and Roberts, G. O. (2018). Piecewise deterministic Markov processes for continuous-time Monte Carlo. Statistical Science 33(3), 386–412.
  * Gilks (1992) Gilks, W. R. (1992). Derivative-free adaptive rejection sampling for Gibbs sampling. Bayesian Statistics 4(4), 641–649.
  * Goan and Fookes (2021) Goan, E. and Fookes, C. (2021). Stochastic bouncy particle sampler for bayesian neural networks. In: Stochastic Bouncy Particle Sampler for Bayesian Neural Networks.
  * Goldman and Singh (2021) Goldman, J. V. and Singh, S. S. (2021). Spatiotemporal blocking of the bouncy particle sampler for efficient inference in state-space models. Statistics and Computing 31(5), 68.
  * Görür and Teh (2011) Görür, D. and Teh, Y. W. (2011). Concave-convex adaptive rejection sampling. Journal of Computational and Graphical Statistics 20(3), 670–691.
  * Hartman (1959) Hartman, P. (1959). On functions representable as a difference of convex functions. Pacific Journal of Mathematics 9, 707–713.
  * Kingman (1992) Kingman, J. F. C. (1992). Poisson Processes. Clarendon Press.
  * Le Thi and Pham Dinh (2018) Le Thi, H. A. and Pham Dinh, T. (2018). DC programming and DCA: thirty years of developments. Mathematical Programming. A Publication of the Mathematical Programming Society 169(1), 5–68.
  * Lewis and Shedler (1979) Lewis, P. A. W. and Shedler, G. S. (1979). Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly 26(3), 403–413.
  * Michel et al. (2014) Michel, M., Kapfer, S. C. and Krauth, W. (2014). Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps. The Journal of Chemical Physics 140(5), 054116\.
  * Michel et al. (2020) Michel, M., Durmus, A. and Sénécal, S. (2020). Forward event-chain Monte Carlo: Fast sampling by randomness control in irreversible Markov chains. Journal of Computational and Graphical Statistics 29, 689–702.
  * Pakman et al. (2017) Pakman, A., Gilboa, D., Carlson, D. and Paninski, L. (2017). Stochastic bouncy particle sampler. In: Proceedings of the 34th International Conference on Machine Learning (eds. D. Precup and Y. W. Teh), volume 70 of Proceedings of Machine Learning Research, PMLR, 2741–2750.
  * Peters and de With (2012) Peters, E. A. and de With, G. (2012). Rejection-free Monte Carlo sampling for general potentials. Physical Review E 85(2), 026703.
  * Sen et al. (2019) Sen, D., Sachs, M., Lu, J. and Dunson, D. (2019). Efficient posterior sampling for high-dimensional imbalanced logistic regression. arXiv:1905.11232 .
  * Stone (1948) Stone, M. H. (1948). The generalized Weierstrass approximation theorem. Mathematics Magazine 21(5), 237–254.
  * Vanetti et al. (2017) Vanetti, P., Bouchard-Côté, A., Deligiannidis, G. and Doucet, A. (2017). Piecewise-Deterministic Markov Chain Monte Carlo. arXiv:1707.05296 .
  * Wu and Robert (2020) Wu, C. and Robert, C. P. (2020). Coordinate sampler: a non-reversible Gibbs-like MCMC sampler. Statistics and Computing 30(3), 721–730.

##  Appendix A Supplementary Material

###  A.1 Simulating from a piecewise linear rate

Suppose the rate function $\lambda(t)=\max\\{0,\ell(t)\\}$ is piecewise linear
on times $t_{0},t_{1},\dots,t_{n}$ with

$\displaystyle\ell(t)=\begin{cases}a_{0}+b_{0}t&0<t<t_{1}\\\
a_{1}+b_{1}t&t_{1}<t<t_{2}\\\ \dots&\dots\\\
a_{n-1}+b_{n-1}t&t_{n-1}<t<t_{n}\end{cases}$ (2)

The event time can be simulated using the following process:

Simulate $u\sim$Uniform[0,1] Set $E=-\log(u)$ For $k=0,...,n-1$: 1\. If
$\int_{t_{k}}^{t_{k+1}}\max\\{0,a_{k}+b_{k}s\\}ds>E$, then return $\tau$
solving:
$\tau=\text{argmin}_{t}\left\\{E=\int_{t_{k}}^{t}\max(0,a_{k}+b_{k}s)ds\right\\}.$
2\. If $\int_{t_{k}}^{t_{k+1}}\max\\{0,a_{k}+b_{k}s\\}ds\leq E$, then update
$E$ as $E=E-\int_{t_{k}}^{t_{k+1}}\max\\{0,a_{k}+b_{k}s\\}ds$ Return

##  Appendix B Alternative GLM likelihood and prior specifications

The likelihood potential for a GLM posterior has the form

$U^{(\ell)}(\bm{\theta})=\sum_{i=1}^{n}\phi(\bm{x}_{i}^{T}\bm{\theta},y_{i})$

where $\phi:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is the GLM mapping function
$\phi(a,y)=\log(y|g^{-1}(a))$. As described in Section 5.1 polynomial rates
can be constructed if the following functions can be evaluated and bounded

$f_{k}^{(j)}(t)=v_{k}\sum_{i=1}^{n}\phi^{(j+1)}(a_{i}(t),y_{i})\left(\frac{\partial
a_{i}}{\partial t}\right)^{j}x_{ik}.$

Functions and derivatives can be found using mathematical software such as
Mathematica.

###  B.1 Robust regression

For robust regression of a response $y\in\mathbb{R}$ we could model the error
with wide tails using the Cauchy likelihood:

$\phi(a,y)=-\log\left(1+\frac{(a-y)^{2}}{b}\right)$

with scale $b>0$. The derivatives are
$\phi^{\prime}(a,y)=-2\frac{a-y}{b+(a-y)^{2}}$,
$\phi^{\prime\prime}(a,y)=2\frac{(a-y)^{2}-b}{(b+(a-y)^{2})^{2}}$ and
$\phi^{\prime\prime\prime}(a,y)=-4\frac{(a-y)((a-y)^{3}-3b)}{(b+(a-y)^{2})^{3}}$
with upper-bounds $|\phi^{\prime}(a,y)|\leq\frac{1}{\sqrt{b}}$,
$|\phi^{\prime\prime}(a,y)|\leq\frac{2}{b}$. For $b\geq 1$
$|\phi^{\prime\prime\prime}(a,y)|\leq 3$ (though better bounds are available).
For $b\leq 1$, the gradient can be bounded but closed form solutions for this
bound are not available.
Mixture of normals:

$\phi(a,y)=-\log\left(0.5\mathcal{N}(y-a;0,1)+0.5\mathcal{N}(y-a;0,10^{2})\right),$

upper-bounds on the derivatives for this function are
$|\phi^{\prime}(a,y)|\leq 2$, $|\phi^{\prime\prime}(a,y)|\leq 0.91$ and
$|\phi^{\prime\prime}(a,y)|\leq 1.8$.

###  B.2 Poisson regression

For $y\sim$Poisson($\exp(a)$) the mapping function is,

$\phi(a,y)=-\log\left(\exp(ya-\exp(a))\right).$

Let $a^{*}=\max\\{a(0),a(\tau_{\max})\\}$ where
$a(t)=\bm{x}^{T}(\bm{\theta}+t\bm{v})$. The derivatives for the mapping
function are $\phi^{(j)}(a,y)\leq\exp(a)$ and upper-bounds on the interval
$\tau\in[0,\tau_{\max}]$ are $\phi^{(j)}(a,y)\leq\exp(a^{*})$.

###  B.3 Gaussian Spike and Slab Prior

Following Chevallier et al. (2020) and Bierkens et al. (2021b), the spike and
slab prior can be simulated using PDMP dynamics. The Gaussian Spike and Slab
prior has the form $\theta_{j}\sim
w_{j}\mathcal{N}(0,\sigma)+(1-w_{j})\delta_{0}$ for $j=1,...,p$ where $w_{j}$
is the prior probability of inclusion and $\delta$ is the Dirac spike at zero.

For RJPDMPs there are three types of events that occur; a reversible jump (RJ)
move if a variable hits the axis ($\theta_{j}$ =0), events according to the
likelihood and prior $\mathcal{N}(0,\sigma)$ for $\theta_{j}$ non-zero, or the
RJ move to reintroduce a variable. The rate to reintroduce a variable is
calculated in Chevallier et al. (2020) where it is found to be

$\frac{p_{rj}}{\sqrt{2\pi\sigma^{2}}}\frac{w}{(1-w)}$

where $p_{rj}$ is a tuning parameter for the probability of jumping to the
reduced model when hitting the axis in Chevallier et al. (2020). The time to
hit an axis for parameter $\theta_{j}$ will be $-\theta_{j}/v_{j}$ when
$\theta_{j}v_{j}>0$ and $\infty$ otherwise. Finally simulating non RJ events
is equivalent to simulating from the Gaussian slab and thus the part of the
event rate contributed by the prior on parameter $\theta_{j}$ is
$v_{j}(\theta_{j}+v_{j}t)/\sigma^{2}$ for $t\in[0,\tau_{\max}]$. This is
linear in time so will be exactly simulated by the CC-PDMP approach.

###  B.4 Cauchy Prior

A Cauchy prior on $\theta$ with scale $b$ and mean $m$ has the form
$p(\theta)\propto\exp(-U(\theta))$ where

$U(\theta)=\log\left(1+\frac{(\theta-m)^{2}}{b}\right)$

the term contributed to the rate depends on

$v\partial_{\theta}U(\theta)=v\frac{2(m-\theta)}{b+(m-\theta)^{2}}$

where $v$ is the velocity corresponding to $\theta$. Following the bounds from
the Cauchy likelihood (B.1) we have

$v\partial_{\theta}U(\theta+tv)\leq f(0)+f^{\prime}(0)t+\frac{M}{2!}t^{2}$

where $f(0)=v\frac{2(m-\theta)}{b+(m-\theta)^{2}}$,
$f^{\prime}(0)=2v\frac{(\theta-m)^{2}-b}{(b+(\theta-m)^{2})^{2}}$ and $M=3$.

###  B.5 PDMP-based samplers on restricted domains

PDMP-based samplers can be implemented on a restricted domain following
Bierkens et al. (2018) or when the density is only piece-wise smooth following
Chevallier et al. (2021). We consider some PDMP samplers on the restricted
domain $\theta>0$. Simulation involves tracking the time to hitting the
boundary $\tau_{b}=\text{inf}\\{t>0:\theta+tv=0\\}$. If an event time occurs
before the boundary is hit the sampler evolves with the usual rate. If the
boundary is hit before the switching event, the sampler is reflected off of
the boundary by sampling a new velocity from a probability measure
concentrated on directions that are normal to the boundary. We refer the
reader to Bierkens et al. (2018) for details.

###  B.6 Generalised inverse Gaussian prior

The generalised inverse Gaussian prior has the form
$p(\theta)\propto\exp(-U(\theta))$ where

$U(\theta)=\theta+\theta^{-1}+2\log(\theta)$

and $\theta>0$. The term contributed to the rate depends on

$v\partial_{\theta}U(\theta)=v-v\theta^{-2}+2\frac{v}{\theta}.$

When $v>0$ we have decomposition $f_{u}(t)=v(1+\frac{2}{\theta+vt})$,
$f_{n}(t)=-\frac{v}{(\theta+tv)^{2}}$ and
$f_{n}^{\prime}(t)=v^{2}\frac{2}{(\theta+vt)^{3}}$.
When $v<0$ we have decomposition $f_{u}(t)=v(1+\frac{1}{(\theta+vt)^{2}})$,
$f_{n}(t)=\frac{v}{(\theta+tv)}$ and
$f_{n}^{\prime}(t)=-v^{2}\frac{2}{(\theta+vt)^{2}}$.

###  B.7 Gamma prior

The Gamma prior has the form $p(\theta)\propto\exp(-U(\theta))$ where

$U(\theta)=\beta\theta-(\alpha-1)\log(\theta)$

and $\theta>0$ for hyper-parameters $\alpha,\beta>0$. The term contributed to
the rate depends on

$v\partial_{\theta}U(\theta)=v\beta-v\frac{(\alpha-1)}{\theta}.$

When $v(\alpha-1)>0$ we have decomposition
$f_{u}(t)=v(\beta+\frac{(\alpha-1)}{\theta+vt})$, $f_{n}(t)=0$.
When $v(\alpha-1)<0$ we have decomposition $f_{u}(t)=v\beta$,
$f_{n}(t)=v\frac{(\alpha-1)}{\theta+vt}$ and
$f_{n}^{\prime}(t)=-v^{2}\frac{(\alpha-1)}{(\theta+vt)^{2}}$.

##  Appendix C Comparison of CC-PDMP and CC-ARS

We note a quick comparison between CC-PDMP and sampling directly with CC-ARS
where computational cost is measured by the number of proposals used in the
adaptive thinning (CC-PDMP) or adaptive sampling (CC-ARS). Both sampling
methods are implemented on a generalised inverse Gaussian distribution (GIG).
The unnormalised GIG density function is

$\pi(x)\propto\exp\left(x+x^{-1}+2\log(x)\right),$

which is log concave for $2\log(x)$ and log convex for $x+x^{-1}$, see section
B.6 for details on the CC-PDMP. We favour the performance of CC-ARS in the
1-dimensional setting since the method yields independent samples while CC-
PDMP has correlated samples. Event times in the PDMP-based method give a
continuous trajectory where samples can be taken at a consistent times along
the trajectory. Rows 1 and 2 of Figure 6 show the methods have similar
performance in exploring the density. Jumping from row 2 to 3 we can see that
the PDMP sampler only needs to go into the tail once to give a good
approximation to the tail of the density. In higher dimensional targets CC-ARS
must be used within Gibbs sampling where the efficiency of the sampler will be
reduced.

Figure 6: Comparison of typical sampling dynamics for PDMP and ARS sampling
using concave-convex upper-bounding in a 1-dimensional GIG distribution. Each
row shows the results for an increasing number of proposals using in the
sampling procedure. Continuous lines are shown for the PDMP dynamics and
points for the ARS method.


# Attribute noise robust binary classification

Aditya Petety1, Sandhya Tripathi2, N Hemachandra2
1National Institute of Science Education and Research, Bhubaneshwar
<EMAIL_ADDRESS>
2Indian Institute of Technology Bombay
2{sandhya.tripathi<EMAIL_ADDRESS>
###### Abstract

We consider the problem of learning linear classifiers when both features and
labels are binary. In addition, the features are noisy, i.e., they could be
flipped with an unknown probability. In Sy-De attribute noise model, where all
features could be noisy together with same probability, we show that $0$-$1$
loss ($l_{0-1}$) need not be robust but a popular surrogate, squared loss
($l_{sq}$) is. In Asy-In attribute noise model, we prove that $l_{0-1}$ is
robust for any distribution over 2 dimensional feature space. However, due to
computational intractability of $l_{0-1}$, we resort to $l_{sq}$ and observe
that it need not be Asy-In noise robust. Our empirical results support Sy-De
robustness of squared loss for low to moderate noise rates.

## Introduction

Quality of data is being compromised as its quantity is getting larger. In
classification setup, bad quality data could be due to noise in the labels or
noise in the features. Label noise research has gained a lot of attention in
last decade (?). In contrast, feature or attribute noise is still unexplored.
As opposed to continuous valued attributes, noise in categorical features,
particularly binary, can drastically change the relative location of a data
point and significantly impact the classifier’s performance.

(?) studied the effect of noise when the algorithms are decision trees. (?; ?)
study attribute noise from the perspective of detecting noisy data points and
correcting them.

Our major contributions lie in identifying loss functions that are robust (or
not) to attribute (binary valued) noise in Empirical Risk Minimization (ERM)
framework. This has an advantage that there is no need of either knowing the
true value or cross-validating over or estimating the noise rates.

## Problem description

Let $D$ be the joint distribution over $\mathbf{X}\times Y$, where
$\mathbf{X}\in\mathcal{X}\subseteq\\{-1,1\\}^{n}$ and
$Y\in\mathcal{Y}=\\{-1,1\\}$. Let the decision function be
$f:\mathbf{X}\mapsto\mathbb{R}$, hypothesis class of all measurable functions
be $\mathcal{H}$ and class of linear hypothesis be
$\mathcal{H}_{lin}=\\{(\bm{\beta},c),\bm{\beta}\in\mathbb{R}^{n},c\in\mathbb{R}:\|\bm{\beta}\|_{2}\leq
B\\}.$ We restrict our set of hypothesis to be in $\mathcal{H}_{lin}$. Let
$\tilde{D}$ denote the distribution on $\tilde{\mathbf{X}}\times{Y}$ obtained
by inducing noise to $D$ with
$\tilde{\mathbf{X}}\in\mathcal{X}\subseteq\\{-1,1\\}^{n}$. The corrupted
sample is
$\tilde{S}:=\\{(\tilde{\mathbf{x}}_{1},{y}_{1}),\ldots,(\tilde{\mathbf{x}}_{m},{y}_{m})\\}\sim\tilde{D}^{m}.$
The probability that the value of $i^{th}$ attribute is flipped is given by
$p_{i}=(\tilde{\mathbf{X}}=-x|\mathbf{X}=x,Y=y),i\in[n]$. We assume that the
class/label does not change with noise in the attributes.

Based on the flipping probability and the dependence between events of
flipping for different attributes, we identify two attribute noise models. If
all the attribute values are flipped together with same probability $p$, then
it is referred to as the symmetric dependent attribute noise model (Sy-De). If
each attribute $j$ flips with probability $p_{j}$ independently of any other
attribute $k\in[n]\backslash\\{j\\}$, then it is referred to as the asymmetric
independent attribute noise model (Asy-In). Even though Sy-De attribute noise
model is simple, it cannot be obtained by taking
$p_{i}=p_{j},\leavevmode\nobreak\ \forall i,j\in[n]$ in Asy-In attribute noise
model. Real world example of Sy-De (or Asy-In) noisy attributes: Consider a
room with many sensors connected in series (or with individual battery)
measuring temperature, humidity, etc., as binary value, i.e., high or low. A
power failure (or battery failures) will lead to all (or individual)
sensors/attributes providing noisy observations with same (or different)
probability.

We consider ERM framework for classification. A natural choice for loss
function is $0$-$1$ loss, i.e.,
$l_{0-1}(f(\mathbf{x},y))=\mathds{1}_{[f(\mathbf{x})y<0]}$. Bayes classifier
$f^{*}=\arg\min_{f\in\mathcal{H}}R_{D}(f)$ and Bayes risk is
$R_{D}(f^{*})=\min_{f\in\mathcal{H}}R_{D}(f)$ where
$R_{D}(f)=E_{D}[\mathds{1}_{[f(\mathbf{x})y<0]}]$. Corresponding quantities
for noisy distribution $\tilde{D}$ are
$R_{\tilde{D}}(\tilde{f}^{*})=\min_{\tilde{f}\in\mathcal{H}}E_{\tilde{D}}[\mathds{1}_{[\tilde{f}(\tilde{\mathbf{x}})y<0]}]$
and
$\tilde{f}^{*}=\arg\min_{\tilde{f}\in\mathcal{H}}R_{\tilde{D}}(\tilde{f})$.

Non-convex nature of $0$-$1$ loss makes it difficult to optimize and hence
convex upper bounds (surrogate losses) are used in practice. In this work, we
consider the squared loss $l_{sq}(f(\mathbf{x},y))=(y-f(\mathbf{x}))^{2}$, a
differentiable and convex surrogate loss function. Our restriction of
hypothesis to linear class $\mathcal{H}_{lin}$ can be interpreted as a form of
regularization. Expected squared clean and corrupted risks are
$R_{D,l_{sq}}(f)=E_{D}[(y-f(\mathbf{x}))^{2}]$ and
$R_{\tilde{D},l_{sq}}(\tilde{f})=E_{\tilde{D}}[(y-\tilde{f}(\tilde{\mathbf{x}}))^{2}]$.
Hypothesis in $\mathcal{H}_{lin}$ minimizing these clean and corrupted risks
are denoted by $f^{*}_{sq,lin}$ and $\tilde{f}^{*}_{sq,lin}$. Next, we define
attribute noise robustness of risk minimization scheme $A_{l}$.

######  Definition 1.

Let $f^{*}_{A_{l}}$ and $\tilde{f}^{*}_{A_{l}}$ be obtained from clean and
corrupted distribution $D$ and $\tilde{D}$ using any arbitrary scheme $A_{l}$.
Then, scheme $A_{l}$ is said to be attribute noise robust if

$R_{D}(f^{*}_{A_{l}})=R_{D}(\tilde{f}^{*}_{A_{l}}).$

Also, $l$ is said to be an attribute noise robust loss function.

## Attribute noise robust loss functions

We, first, consider Sy-De attribute noise model and present a counter example
(Example 1) to show that $0$-$1$ loss need not be robust to Sy-De attribute
noise. To circumvent this problem, we provide a positive result by showing
that squared loss is Sy-De attribute noise robust with origin passing linear
classifiers (Theorem 1). Our hypothesis set belongs to $\mathcal{H}_{lin}$
which could be further categorized into origin passing and non-origin passing
($c=0$ or not). Details of examples and proofs are available in Supplementary
Material (SM).

######  Example 1.

Consider a population of two data points (in 1-D) $(x,y)$ as $(-1,1)$ and
$(1,-1)$ with probability $0.25$ and $0.75$ with a classifier
$f_{lin}(x)=bx+c$. Then, the $l_{0-1}$ optimal clean classifier is
$f_{lin}^{*}=(b^{*},c^{*})=(-1,-0.1)$ with $R_{D}(f_{lin}^{*})=0$. Also, the
$l_{0-1}$ optimal Sy-De attribute noise ($p=0.4$) corrupted classifier is
$\tilde{f}^{*}_{lin}=(\tilde{b}^{*},\tilde{c}^{*})=(1,-2)$ with
$R_{D}(\tilde{f}_{lin}^{*})=0.25$. Since, $R_{D}(\tilde{f}_{lin}^{*})\neq
R_{D}(f_{lin}^{*})$, $0$-$1$ loss function need not be Sy-De attribute noise
robust.

######  Theorem 1.

Consider a clean distribution $D$ on $\mathbf{X}\times Y$ and Sy-De attribute
noise corrupted distribution $\tilde{D}$ on $\tilde{\mathbf{X}}\times Y$ with
noise rate $p<0.5$. Then, squared loss $l_{sq}$ with origin passing linear
classifiers is Sy-De attribute noise robust,i.e.,

$R_{D}(\tilde{f}_{lin,l_{sq}}^{*})=R_{D}({f}_{lin,l_{sq}}^{*})$ (1)

where $f^{*}_{lin,l_{sq}}=(\beta_{1}^{*},\ldots,\beta_{n}^{*})$ and
$\tilde{f}^{*}_{lin,l_{sq}}=(\tilde{\beta}_{1}^{*},\ldots,\tilde{\beta}_{n}*)$
correspond to optimal linear classifiers learnt using squared loss on clean
($D$) and corrupted ($\tilde{D})$ distribution.

######  Remark 1.

Sy-De robustness of squared loss is an interesting result because given an
attribute noise corrupted dataset, obtaining a linear classifier entails
solving only a linear system of equations. (Demonstrated on UCI datasets.)

Now, we consider Asy-In attribute noise model and show that $0$-$1$ loss is
robust to this noise with non-origin passing classifiers when $n=2$ (Theorem
2). As $l_{0-1}$ based ERM is computationally intractable, we consider
$l_{sq}$ and present a counter example to show that $l_{sq}$ need not be Asy-
In noise robust (Example 2).

######  Theorem 2.

Consider a clean distribution $D$ with probabilities
$\\{d_{1},d_{2},d_{3},d_{4}\\}$ on $\mathbf{X}\times Y$ with $n=2$ (population
of $2^{n}$ data points) and Asy-In attribute noise corrupted distribution
$\tilde{D}$ on $\tilde{\mathbf{X}}\times Y$ with noise rates $p_{1}<0.5$ and
$p_{2}<0.5$. Then, $0$-$1$ loss with non-origin passing linear classifiers is
Asy-In attribute noise robust, i.e.,

$R_{D}(\tilde{f}_{lin,l_{0-1}}^{*})=R_{D}({f}_{lin,l_{0-1}}^{*})$ (2)

where $f^{*}_{lin,l_{0-1}}=(\beta_{1}^{*},1,c^{*})$ and
$\tilde{f}^{*}_{lin,l_{0-1}}=(\tilde{\beta}_{1}^{*},1,\tilde{c}^{*})$
correspond to optimal linear classifiers learnt using $0$-$1$ loss on clean
($D$) and corrupted ($\tilde{D})$ distribution respectively.

######  Example 2.

Consider a population of 3 data points (in 2-D) $(x_{1},x_{2},y)$ as
$(-1,1,1)$, $(-1,-1,-1)$, and $(1,-1,1)$ with probabilities as
$(\frac{1}{4},\frac{1}{2},\frac{1}{4})$ with a classifier
$f_{lin}(\mathbf{x})=b_{1}x_{1}+b_{2}x_{2}$. Then, $l_{sq}$ optimal clean
classifier is $f_{lin,l_{sq}}^{*}=(b_{1}^{*},b_{2}^{*})=(0.5,0.5)$ with
$R_{D}(f_{lin,l_{sq}}^{*})=0.5$. Also, $l_{sq}$ optimal Asy-In attribute noise
($p_{1}=0.1,\leavevmode\nobreak\ p_{2}=0.2$) corrupted classifier is
$\tilde{f}_{lin,l_{sq}}^{*}=(\tilde{b}_{1}^{*},\tilde{b}_{2}^{*})=(0.4,0.3)$
with $R_{D}(\tilde{f}_{lin,l_{sq}}^{*})=0.25$. Since,
$R_{D}(f_{lin,l_{sq}}^{*})\neq R_{D}(\tilde{f}_{lin,l_{sq}}^{*})$, squared
loss need not be Asy-In attribute noise robust.

## Experiments

Figure 1 demonstrates Sy-De attribute noise robustness of squared loss on 3
UCI datasets (?); details in SM. As SPECT dataset is imbalanced, in addition
to accuracy, we also report arithmetic mean (AM). To account for randomness in
noise, results are averaged over 15 trials of train-test partitioning (80-20).
The low accuracy in comparison to clean classifier can be attributed to the
finite samples available for learning the classifiers.

Figure 1:  Test data performance of $l_{sq}$ with Sy-De attribute noise.

## Looking forward

Our work is an initial attempt in binary valued attribute noise; an extension
to general discrete valued attributes would be interesting. Asy-In attribute
noise model raises some non-trivial questions w.r.t. choice of loss functions
like robustness of $0$-$1$ for $n>2$, explanation for the surprising non-
robustness of squared loss as compared to robustness of a difficult to deal
$0$-$1$ loss, search for other surrogate loss functions that are robust.
Finally, we believe that attribute dimension $n$ could have a role to play in
noise robustness.

## References

  * [Dheeru and Karra Taniskidou, 2017] Dheeru, D. and Karra Taniskidou, E. (2017).  UCI machine learning repository.
  * [Khoshgoftaar and Van Hulse, 2009] Khoshgoftaar, T. M. and Van Hulse, J. (2009).  Empirical case studies in attribute noise detection.  IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews), 39(4):379–388.
  * [Quinlan, 1986] Quinlan, J. R. (1986).  The effect of noise on concept learning.  Machine learning: An artificial intelligence approach, 2:149–166.
  * [Sastry and Manwani, 2016] Sastry, P. S. and Manwani, N. (2016).  Robust Learning of Classifiers in the Presence of Label Noise, chapter 6, pages 167–197.
  * [Zhu and Wu, 2004] Zhu, X. and Wu, X. (2004).  Class noise vs. attribute noise: A quantitative study.  Artificial Intelligence Review, 22(3):177–210.

Supplementary material

##  Appendix A Proofs

### Proof of Theorem 1

###### Proof.

Consider squared loss based clean risk given as follows:

$R_{D,l_{sq}}(f)=\mathbb{E}_{\mathbf{X}\times
Y}[(Y-\bm{\beta}^{T}\mathbf{X})^{2}]$

We minimize this by differentiating the expectation term and equating it to 0.

$\mathbb{E}_{\mathbf{X}\times
Y}[2(Y-\bm{\beta}^{T}\mathbf{X})\mathbf{X}^{T}]=0$
$\implies\mathbb{E}_{\mathbf{X}\times
Y}[\mathbf{X}^{T}Y-\bm{\beta}^{T}\mathbf{X}\mathbf{X}^{T}]=0$
$\implies\bm{\beta}^{*}=(\mathbb{E}[\mathbf{XX}^{T}]^{-1})^{T}\mathbb{E}[Y\mathbf{X}]$

Now, to obtain the optimal noisy classifier $\tilde{f}^{*}_{lin,l_{sq}}$ we
minimize the following squared loss based expected risk.

$\displaystyle R_{\tilde{D},l_{sq}}(\tilde{f})$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\tilde{\mathbf{X}}\times
Y}[(Y-\tilde{\bm{\beta}}^{T}\tilde{\mathbf{X}})^{2}]$ (3) $\displaystyle=$
$\displaystyle(1-p)\mathbb{E}_{\mathbf{X}\times
Y}[(Y-\tilde{\bm{\beta}}^{T}\mathbf{X})^{2}]$
$\displaystyle+p\mathbb{E}_{\mathbf{X}\times
Y}[(Y+\tilde{\bm{\beta}}^{T}\mathbf{X})^{2}]$

We minimize the corrupted risk given in equation (3) by differentiating and
equating to 0 as described below:

$(1-p)\mathbb{E}_{\mathbf{X}\times
Y}[(Y-\tilde{\bm{\beta}}^{T}\mathbf{X})(-2\mathbf{X}^{T})]$
$+p\mathbb{E}_{\mathbf{X}\times
Y}[(Y+\tilde{\bm{\beta}}^{T}\mathbf{X})(2\mathbf{X}^{T})]=0$
$\implies(p-1)\mathbb{E}[\mathbf{X}^{T}Y]-(p-1)\tilde{\bm{\beta}}^{T}\mathbb{E}[\mathbf{X}\mathbf{X}^{T}]$
$+p\tilde{\bm{\beta}}^{T}\mathbb{E}[\mathbf{X}\mathbf{X}^{T}]+p\mathbb{E}[\mathbf{X}^{T}Y]=0$
$\implies(2p-1)\mathbb{E}[\mathbf{X}^{T}Y]+\tilde{\bm{\beta}}^{T}\mathbb{E}[\mathbf{X}\mathbf{X}^{T}]=0$
$\implies\tilde{\bm{\beta}}^{*}=(1-2p)(\mathbb{E}[\mathbf{XX}^{T}]^{-1})^{T}\mathbb{E}[Y\mathbf{X}]$
$\implies\tilde{\bm{\beta}}^{*}=(1-2p)\bm{\beta}^{*}$

We can see that the noisy classifier $\tilde{f}^{*}_{lin,l_{sq}}$ is just a
scaled version of the classifier, $f^{*}_{lin,l_{sq}}$, obtained from the
clean risk. Since, for attribute noise robustness it is sufficient that
$sign(\tilde{\bm{\beta}}^{*}\mathbf{x})=sign(\bm{\beta}^{*}\mathbf{x})$, the
aforementioned observation proves that the squared loss function is Sy-De
attribute noise robust when $p<0.5$. ∎

### Proof of Theorem 2

###### Proof.

We prove the robustness of $0$-$1$ loss in 2 dimension by taking an exhaustive
search approach. Even though we consider a particular non-uniform distribution
$\\{d_{1},d_{2},d_{3},d_{4}\\}$ over 4 data points in the population and
certain values of $p_{1}$ and $p_{2}$, the following claims hold for any
arbitrary value of $d_{i}\in[0,1],\leavevmode\nobreak\ \leavevmode\nobreak\
\sum\limits_{i=1}^{4}d_{i}=0,\leavevmode\nobreak\ i=1,\ldots,4$ and any pair
of noise rate $p_{1},p_{2}<0.5$. In the population, there are 16 ways in which
four data points can be assigned to two classes. By symmetry, we can restrict
ourselves to just these 4 cases.

The probabilities of $(\mathbf{x}^{(1)},y^{(1)})$,
$(\mathbf{x}^{(2)},y^{(2)})$, $(\mathbf{x}^{(3)},y^{(3)})$,
$(\mathbf{x}^{(4)},y^{(4)})$ in the distribution $D$ are taken to be
$0.25,0.33,0.39$ and $0.03$ respectively. Here,
$\mathbf{x}^{(i)}=[x_{i1},x_{i2}],\leavevmode\nobreak\ i=1,\ldots,4$. The
noise rates are $p_{1}=0.12$ and $p_{2}=0.23$. We consider a classifier of the
form $f(\mathbf{x})=\beta\mathbf{x}+c$ where $\beta=[b1,1]$.

The clean $0$-$1$ risk denoted by $R$ is given as follows:

$\displaystyle
R=R_{D}(f_{lin})=\mathbb{E}_{D}[\mathds{1}_{[(\beta\mathbf{x}+c)y<0]}]=0.25\mathds{1}_{[(b1x_{11}+x_{12}+c)y_{<}0]}$
$\displaystyle+0.33\mathds{1}_{[(b1x_{21}+x_{22}+c)y_{2}<0]}+0.39\mathds{1}_{[(b1x_{31}+x_{32}+c)y_{3}<0]}$
$\displaystyle+0.03\mathds{1}_{[(b1x_{41}+x_{42}+c)y_{4}<0]}$

The noisy $0$-$1$ risk denoted by $R^{\prime}$ is given as follows:

$\displaystyle
R^{\prime}=R_{\tilde{D}}(\tilde{f}_{lin})=\mathbb{E}_{\tilde{D}}[\mathds{1}_{[(\tilde{\mathbf{\beta}}\tilde{\mathbf{x}}+c^{\prime})y<0]}]$
$\displaystyle=(1-p_{1})(1-p_{2})\big{[}0.25\mathds{1}_{[(b1^{\prime}x_{11}+x_{12}+c^{\prime})y_{<}0]}+0.33\mathds{1}_{[(b1^{\prime}x_{21}+x_{22}+c^{\prime})y_{2}<0]}$
$\displaystyle+0.39\mathds{1}_{[(b1^{\prime}x_{31}+x_{32}+c^{\prime})y_{3}<0]}+0.03\mathds{1}_{[(b1^{\prime}x_{41}+x_{42}+c^{\prime})y_{4}<0]}\big{]}$
$\displaystyle+p_{1}p_{2}[0.25\mathds{1}_{[(-b1^{\prime}x_{11}-x_{12}+c^{\prime})y_{<}0]}+0.33\mathds{1}_{[(-b1^{\prime}x_{21}-x_{22}+c^{\prime})y_{2}<0]}$
$\displaystyle+0.39\mathds{1}_{[(-b1^{\prime}x_{31}-x_{32}+c^{\prime})y_{3}<0]}+0.03\mathds{1}_{[(-b1^{\prime}x_{41}-x_{42}+c^{\prime})y_{4}<0]}]$
$\displaystyle+(1-p_{1})p_{2}[0.25\mathds{1}_{[(b1^{\prime}x_{11}-x_{12}+c^{\prime})y_{<}0]}+0.33\mathds{1}_{[(b1^{\prime}x_{21}-x_{22}+c^{\prime})y_{2}<0]}$
$\displaystyle+0.39\mathds{1}_{[(b1^{\prime}x_{31}-x_{32}+c^{\prime})y_{3}<0]}+0.03\mathds{1}_{[(b1^{\prime}x_{41}-x_{42}+c^{\prime})y_{4}<0]}]$
$\displaystyle+p_{1}(1-p_{2})[0.25\mathds{1}_{[(-b1^{\prime}x_{11}+x_{12}+c^{\prime})y_{<}0]}+0.33\mathds{1}_{[(-b1^{\prime}x_{21}+x_{22}+c^{\prime})y_{2}<0]}$
$\displaystyle+0.39\mathds{1}_{[(-b1^{\prime}x_{31}+x_{32}+c^{\prime})y_{3}<0]}+0.03\mathds{1}_{[(-b1^{\prime}x_{41}+x_{42}+c^{\prime})y_{4}<0]}]$

To find the minimizer of the clean and corrupted risks $R$ and $R^{\prime}$,
we plot them (in MATLAB) as a function of $(b1,c)$ and
$(b1^{\prime},c^{\prime})$. In the 4 cases considered below, we observed that
even though the minimum value for $R$ and $R^{\prime}$ are different, they
have same minimizers. This implies that the clean $0$-$1$ risk corresponding
to the optimal clean and noisy classifier would be same. And hence, condition
for attribute noise robustness in Definition 1 is satisfied.

Next, we provide the details for each case.

$\text{{Case 1} }(\mathbf{x}^{(1)},y^{(1)})=(\begin{bmatrix}+1\\\
+1\end{bmatrix},+1)\quad(\mathbf{x}^{(2)},y^{(2)})=(\begin{bmatrix}+1\\\
-1\end{bmatrix},-1)$ $(\mathbf{x}^{(3)},y^{(3)})=(\begin{bmatrix}-1\\\
-1\end{bmatrix},+1)\quad(\mathbf{x}^{(4)},y^{(4)})=(\begin{bmatrix}-1\\\
+1\end{bmatrix},-1)$

The set of classifiers obtained by minimizing the 0-1 risk is given in Fig 2.

Figure 2: Clean $0$-$1$ risk $R$ for case 1.

Now, the set of classifiers obtained by minimizing the 0-1 risk is given in
Figure 3.

Figure 3: Noisy $0$-$1$ risk $R^{\prime}$ for case 1.

We can see that $f^{*}=-4x_{1}+x_{2}+4$ is a classifier which minimizes the
clean as well as noisy risk.

$\text{{Case 2} }(\mathbf{x}^{(1)},y^{(1)})=(\begin{bmatrix}+1\\\
+1\end{bmatrix},+1)\quad(\mathbf{x}^{(2)},y^{(2)})=(\begin{bmatrix}+1\\\
-1\end{bmatrix},-1)$ $(\mathbf{x}^{(3)},y^{(3)})=(\begin{bmatrix}-1\\\
-1\end{bmatrix},-1)\quad(\mathbf{x}^{(4)},y^{(4)})=(\begin{bmatrix}-1\\\
+1\end{bmatrix},+1)$

The set of classifiers obtained by minimizing the 0-1 risk is given in Fig 4.

Figure 4: Clean $0$-$1$ risk $R$ for case 2.

Now, the set of classifiers obtained by minimizing the 0-1 risk is given in
Figure 5.

Figure 5: Noisy $0$-$1$ risk $R^{\prime}$ for case 2.

We can see that $f^{*}=4x_{1}+x_{2}-4$ is a classifier which minimizes the
clean as well as noisy risk.

$\text{{Case 3} }(\mathbf{x}^{(1)},y^{(1)})=(\begin{bmatrix}+1\\\
+1\end{bmatrix},+1)\quad(\mathbf{x}^{(2)},y^{(2)})=(\begin{bmatrix}+1\\\
-1\end{bmatrix},-1)$ $(\mathbf{x}^{(3)},y^{(3)})=(\begin{bmatrix}-1\\\
-1\end{bmatrix},-1)\quad(\mathbf{x}^{(4)},y^{(4)})=(\begin{bmatrix}-1\\\
+1\end{bmatrix},-1)$

The set of classifiers obtained by minimizing the 0-1 risk is given in Fig 6.

Figure 6: Clean $0$-$1$ risk $R$ for case 3.

Now, the set of classifiers obtained by minimizing the 0-1 risk is given in
Figure 7.

Figure 7: Noisy $0$-$1$ risk $R^{\prime}$ for case 3.

We can see that $f^{*}=4x_{1}+x_{2}-4$ is a classifier which minimizes the
clean as well as noisy risk.

$\text{{Case 4} }(\mathbf{x}^{(1)},y^{(1)})=(\begin{bmatrix}+1\\\
+1\end{bmatrix},-1)\quad(\mathbf{x}^{(2)},y^{(2)})=(\begin{bmatrix}+1\\\
-1\end{bmatrix},-1)$ $(\mathbf{x}^{(3)},y^{(3)})=(\begin{bmatrix}-1\\\
-1\end{bmatrix},-1)\quad(\mathbf{x}^{(4)},y^{(4)})=(\begin{bmatrix}-1\\\
+1\end{bmatrix},-1)$

The set of classifiers obtained by minimizing the 0-1 risk is given in Fig 8.

Figure 8: Clean $0$-$1$ risk $R$ for case 4.

Now, the set of classifiers obtained by minimizing the 0-1 risk is given in
Figure 9.

Figure 9: Noisy $0$-$1$ risk $R^{\prime}$ for case 4.

We can see that $f^{*}=x_{2}-3.9$ is a classifier which minimizes the clean as
well as noisy risk.

As seen in the above four cases, we conclude that the set of minimizers of
clean $0$-$1$ risk and corrupted $0$-$1$ are same. And hence, the classifiers
learnt by minimizing the 0-1 risk noisy distributions are Asy-In attribute
noise robust with $n=2$.

∎

##  Appendix B Details of counter examples

### Counter example 1

Consider a population of two data points of the form $(x,y)$ in 1-dimension as
$(-1,1)$ and $(1,-1)$ with probability $0.25$ and $0.75$ and consider the
classifier of the form $f_{lin}(x)=bx+c$. Then, the clean $0$-$1$ risk is
given as follows:

$R_{D}(f_{lin})=0.25\mathds{1}_{[-b+c<0]}+0.75\mathds{1}_{[b+c>0]}$

Minimizing $R_{D}(f_{lin})$ w.r.t. $b$ and $c$ gives us the optimal clean
classifier and optimal clean risk as

$f_{lin}^{*}=(b^{*},c^{*})=(-1,-0.1),$ $R_{D}(f_{lin}^{*})=0.$

Next, we consider the Sy-De attribute noise corrupted risk with $p=0.4$ as
follows:

$\displaystyle R_{\tilde{D}}(\tilde{f}_{lin})$ $\displaystyle=$
$\displaystyle(1-p)[0.25\mathds{1}_{[-\tilde{b}+\tilde{c}<0]}+0.75\mathds{1}_{[\tilde{b}+\tilde{c}>0]}]$
(4)
$\displaystyle+p[0.25\mathds{1}_{[\tilde{b}+\tilde{c}<0]}+0.75\mathds{1}_{[-\tilde{b}+\tilde{c}>0]}]$

Minimizing $R_{\tilde{D}}(\tilde{f}_{lin})$ w.r.t. $\tilde{b}$ and $\tilde{c}$
gives us the following optimal noisy classifier and $0$-$1$ clean risk as

$\tilde{f}^{*}_{lin}=(\tilde{b}^{*},\tilde{c}^{*})=(1,-2)$
$R_{D}(\tilde{f}_{lin}^{*})=0.25$

Clearly, $R_{D}(f_{lin}^{*})\neq R_{D}(\tilde{f}_{lin}^{*})$ implying that
$0$-$1$ loss need not be Sy-De attribute noise robust.

### Counter example 2

Consider a population of 3 data points of the form $(x_{1},x_{2},y)$ in 2
dimensions as $(-1,1,1),(-1,-1,-1)$ and $(1,-1,+1)$ with probabilities as
$(\frac{1}{4},\frac{1}{2},\frac{1}{4})$. We consider the classifier of the
form $f_{lin}(\mathbf{x})=b_{1}x_{1}+b_{2}x_{2}$. Then, the clean squared loss
based risk is given as follows:

$\displaystyle R_{D,l_{sq}}(f_{lin,sq})$ $\displaystyle=$
$\displaystyle\frac{1}{4}(-b_{1}+b_{2}-1)^{2}+\frac{1}{2}(-b_{1}-b_{2}+1)^{2}$
$\displaystyle+\frac{1}{4}(b_{1}-b_{2}-1)^{2}$

Minimizing $R_{D,l_{sq}}(f_{lin,l_{sq}})$ w.r.t. $b_{1}$ and $b_{2}$ leads to
following system of equations:

$\displaystyle 2b_{1}-1$ $\displaystyle=$ $\displaystyle 0$ $\displaystyle
2b_{2}-1$ $\displaystyle=$ $\displaystyle 0$

Solving the above system of equations gives us the optimal squared loss based
clean classifier and clean $0$-$1$ risk as

$f^{*}_{lin,l_{sq}}=(b_{1}^{*},b_{2}^{*})=(0.5,0.5)$
$R_{D}(f^{*}_{lin,l_{sq}})=\frac{1}{2}$

Next, we consider the Asy-In attribute noise corrupted risk in terms of
$p_{1}$ and $p_{2}$ as follows:

$\displaystyle R_{\tilde{D},l_{sq}}(\tilde{f}_{lin,l_{sq}})$
$\displaystyle=(1-p_{1})(1-p_{2})\bigg{[}\frac{1}{4}(-\tilde{b}_{1}+\tilde{b}_{2}-1)^{2}+\frac{1}{2}(-\tilde{b}_{1}-\tilde{b}_{2}+1)^{2}$
$\displaystyle+\frac{1}{4}(\tilde{b}_{1}-\tilde{b}_{2}-1)^{2}\bigg{]}+p_{1}p_{2}\Bigg{[}\frac{1}{4}(\tilde{b}_{1}-\tilde{b}_{2}-1)^{2}+\frac{1}{2}(\tilde{b}_{1}+\tilde{b}_{2}+1)^{2}$
$\displaystyle+\frac{1}{4}(-\tilde{b}_{1}+\tilde{b}_{2}-1)^{2}\bigg{]}+(1-p_{1})p_{2}\Bigg{[}\frac{1}{4}(-\tilde{b}_{1}-\tilde{b}_{2}-1)^{2}$
$\displaystyle+\frac{1}{2}(-\tilde{b}_{1}+\tilde{b}_{2}+1)^{2}+\frac{1}{4}(\tilde{b}_{1}+\tilde{b}_{2}-1)^{2}\bigg{]}$
$\displaystyle+p_{1}(1-p_{2})\Bigg{[}\frac{1}{4}(\tilde{b}_{1}+\tilde{b}_{2}-1)^{2}+\frac{1}{2}(\tilde{b}_{1}-\tilde{b}_{2}+1)^{2}$
$\displaystyle+\frac{1}{4}(-\tilde{b}_{1}-\tilde{b}_{2}-1)^{2}\Bigg{]}$

Minimizing $R_{\tilde{D},l_{sq}}(\tilde{f}_{lin,l_{sq}})$ w.r.t.
$\tilde{b}_{1},\tilde{b}_{2}$ leads to the following system of equations:

$\displaystyle 2\tilde{b}_{1}$ $\displaystyle=$ $\displaystyle(1-2p_{1})$
$\displaystyle 2\tilde{b}_{2}$ $\displaystyle=$ $\displaystyle(1-2p_{2})$

Now, if we take the noise rate values as $p_{1}=0.1$ and $p_{2}=0.2$, then
solving for above system of equations gives us the optimal squared loss based
noisy linear classifier and clean $0$-$1$ risk as

$\tilde{f}^{*}_{lin,l_{sq}}=(\tilde{b}_{1}^{*},\tilde{b}_{2}^{*})=(0.4,0.3)$
$R_{D}(\tilde{f}^{*}_{lin,l_{sq}})=\frac{1}{4}$

Clearly, $R_{D}(f_{lin,l_{sq}}^{*})\neq R_{D}(\tilde{f}_{lin,l_{sq}}^{*})$
implying that squared loss need not be Asy-In attribute noise robust.

##  Appendix C UCI dataset details

In this section, we provide some details on the way we processed the datasets
to obtain the experimental results. Number of features and the number of data
points (number of negative and positive labelled separately) are provided in
Table 1. We provide individual pre-processing details for each dataset as
follows:

  * •

Vote: This is a voting dataset with original size of 435 data points. Since
missing entry in a cell meant that the person takes a neutral stand, we
removed such instances to fit in the framework of binary valued attributes and
finally used 232 data points. Label $+1$ corresponds to “Democrat” and label
$-1$ corresponds to ”Republican”.

  * •

SPECT: This dataset has information extracted from cardiac SPECT images with
values as $0$ and $1$. To be consistent with the binary format, we replaced
all $0$’s by $-1$ without loss of generality in both attributes and labels.

  * •

KR-vs-KP: This dataset originally has 36 features but we removed feature
number 15 as it had three categorical values and finally used the dataset with
35 attributes. Also, the attributes are processed as follows: “f” replaced by
“+1”, “t” replaced by “-1”, “n” replaced by “+1”, “g” replaced by “+1”, “l”
replaced by “-1”. Finally, label $+1$ corresponds to “Won” and label $-1$
corresponds to “Nowin”.

S.no | Dataset name | n | m ($m_{+}$,$m_{-}$)
---|---|---|---
1 | Vote | 16 | 232 (124,108)
2 | SPECT | 22 | 267 (212,55)
3 | KR-vs-KP | 35 | 3196 (1569,1527)
Table 1: Details about the number of features (n), number of total data points
(m), number of positively labelled data points ($m_{+}$) and number of
negatively labelled data points ($m_{-}$).

In Tables 2 and 3, we present the values which we used to generate the plots
in Figure 1. The results are averaged over 15 trials. We observe that at high
noise rates, the theoretically proven robustness of squared loss doesn’t work
because of the finite samples.

Dataset | Vote | SPECT | KR-vs-KP
---|---|---|---
(m,n) | (232 (124,108), 16) | (267(212,55),22) | (3196(1569,1527),35)
| Mean | SD | Mean | SD | Mean | SD
p | 0.969 | 0.0218 | 0.732 | 0.0589 | 0.940 | 0.0109
0 | 0.966 | 0.0203 | 0.693 | 0.0746 | 0.938 | 0.0096
0.1 | 0.955 | 0.0363 | 0.680 | 0.0681 | 0.924 | 0.0106
0.2 | 0.865 | 0.0839 | 0.638 | 0.0797 | 0.900 | 0.0150
0.3 | 0.821 | 0.0941 | 0.625 | 0.0664 | 0.866 | 0.0193
0.35 | 0.694 | 0.0919 | 0.579 | 0.0955 | 0.795 | 0.0263
0.4 | 0.573 | 0.1041 | 0.505 | 0.0793 | 0.706 | 0.0393
Table 2: Average (Mean) test accuracy along with standard deviation (SD) over 15 trials obtained by using squared loss based linear classifier learnt on Sy-De (noise rate p) attribute noise corrupted data. Even though squared loss is theoretically shown to be Sy-De noise robust, it doesn’t show good performance at high noise rates. This could be because the result in Theorem 1 is in expectation. In particular, finite sample size starts showing its effect at high noise rates and the performance deteriorates. Dataset | Vote | SPECT | KR-vs-KP
---|---|---|---
(m,n) | (232 (124,108), 16) | (267(212,55),22) | (3196(1569,1527),35)
| Mean | SD | Mean | SD | Mean | SD
p | 0.970 | 0.0207 | 0.665 | 0.0639 | 0.940 | 0.0111
0 | 0.967 | 0.0193 | 0.613 | 0.1214 | 0.937 | 0.0097
0.1 | 0.956 | 0.0373 | 0.605 | 0.1338 | 0.923 | 0.0109
0.2 | 0.868 | 0.0839 | 0.568 | 0.1353 | 0.899 | 0.0152
0.3 | 0.823 | 0.0950 | 0.586 | 0.1188 | 0.866 | 0.0186
0.35 | 0.698 | 0.0903 | 0.544 | 0.1542 | 0.794 | 0.0261
0.4 | 0.575 | 0.1045 | 0.540 | 0.1472 | 0.706 | 0.0399
Table 3: Average (Mean) test AM value along with standard deviation (SD) over
15 trials obtained by using squared loss based linear classifier learnt on Sy-
De (noise rate p) attribute noise corrupted data. Due to the imbalanced nature
of SPECT dataset, AM is a more suitable evaluation metric.

##  Appendix D Additional examples

######  Example 3.

This is another example which demonstrates that $0$-$1$ is robust to Asy-In
attribute noise with $n=2$. Let the input data be two dimensional and the
clean training set be

$(\mathbf{x}^{(1)},y^{(1)})=(\begin{bmatrix}+1\\\
+1\end{bmatrix},+1)\quad(\mathbf{x}^{(2)},y^{(2)})=(\begin{bmatrix}+1\\\
-1\end{bmatrix},-1)$ $(\mathbf{x}^{(3)},y^{(3)})=(\begin{bmatrix}-1\\\
-1\end{bmatrix},-1)$

and uniformly distributed with
$\mathbf{x}^{(i)}=[x_{i1},x_{i2}],\leavevmode\nobreak\ i=1,2,3$. Let the
flipping probabilities of the first and second component be $p_{1}=0.12$ and
$p_{2}=0.23$ respectively. Let us consider the loss function to be the $0$-$1$
loss and the classifier be of the form
$f_{lin}(\mathbf{x}^{(i)})=b_{1}x_{i1}+x_{i2}+c,i=1,2,3$. We calculate
$f^{*}_{lin}$ as

$f_{lin}^{*}=\arg\min\limits_{f_{lin}}\mathbb{E}_{D}[l_{0-1}(y,f_{lin})]$

We get a range of values for $b_{1}^{*}$ and $c^{*}$, we can choose one of the
value which is $b_{1}^{*}=1/3=c^{*}$ that gives the risk to be 0.

Now, for the corrupted case, we minimize the noisy $0$-$1$ classifier to
obtain
$\tilde{f^{*}_{lin}}(\tilde{x}^{(i)})=\tilde{b}_{1}\tilde{x}_{i1}+\tilde{x}_{i2}+\tilde{c},\leavevmode\nobreak\
i=1,2,3$ given as below:

$\displaystyle\tilde{f^{*}_{lin}}=\arg\min\limits_{\tilde{f^{*}}_{lin}}\mathbb{E}_{\tilde{D}}[l_{0-1}(y,\tilde{f}_{lin})]$
$\displaystyle=\arg\min\limits_{\tilde{b}_{1},\tilde{c}}\frac{1}{3}(1-p_{1})(1-p_{2})\Big{[}\mathds{1}_{[\tilde{b}_{1}+1+\tilde{c}<0,y=1]}$
$\displaystyle+\mathds{1}_{[\tilde{b}_{1}-1+\tilde{c}>0,y=-1]}+\mathds{1}_{[-\tilde{b}_{1}-1+\tilde{c}>0,y=-1]}\Big{]}$
$\displaystyle+\frac{1}{3}\Big{[}p_{1})(p_{2})(\mathds{1}_{[-\tilde{b}_{1}-1+\tilde{c}<0,y=1]}+\mathds{1}_{[-\tilde{b}_{1}+1+\tilde{c}>0,y=-1]}$
$\displaystyle+\mathds{1}_{[\tilde{b}_{1}+1+\tilde{c}>0,y=-1]}\Big{]}+\frac{1}{3}(1-p_{1})(p2)\Big{[}\mathds{1}_{[\tilde{b}_{1}-1+\tilde{c}<0,y=1]}$
$\displaystyle+\mathds{1}_{[\tilde{b}_{1}+1+\tilde{c}>0,y=-1]}+\mathds{1}_{[-\tilde{b}_{1}+1+\tilde{c}>0,y=-1]}\Big{]}$
$\displaystyle+\frac{1}{3}(p_{1})(1-p_{2})\Big{[}\mathds{1}_{[-\tilde{b}_{1}+1+\tilde{c}<0,y=1]}+\mathds{1}_{[-\tilde{b}_{1}-1+\tilde{c}>0,y=-1]}$
$\displaystyle+\mathds{1}_{[\tilde{b}_{1}-1+\tilde{c}>0,y=-1]}\Big{]}$

The minimizer for the above equation is calculated by plotting a graph (in
MATLAB) given in Figure 10. Here, the values of $p_{1}$ and $p_{2}$ are taken
to be 0.12 and 0.23 respectively. The same pattern is observed for all values
of $p_{1},p_{2}<0.5$

Figure 10: Noisy $0$-$1$ risk is plotted on $z$ axis. The $x$ and $y$ axis
labels are to be read as $\tilde{b}_{1}$ and $\tilde{c}$ as we are looking for
the minimizers of noisy risk.

The minimum value is obtained at $\tilde{b}_{1}^{*}=3$ and $\tilde{c}^{*}=-3$.

Comparing the clean $0$-$1$ risk of the classifiers $f^{*}_{lin}$ and
$\tilde{f^{*}}_{lin}$, we observe that they are equal(0) and hence the $0$-$1$
loss function is attribute noise robust in this case.