Pychastic: Precise Brownian Dynamics using Taylor-Itō integrators in Python

Radost Waszkiewicz1$\star$, Maciej Bartczak1, Kamil Kolasa1, and Maciej
Lisicki1

1 Institute of Theoretical Physics, Faculty of Physics, University of Warsaw

L. Pasteura 5, 02-093 Warsaw, Poland

⋆<EMAIL_ADDRESS>
## Abstract

In the last decade, Python-powered physics simulations ecosystem has been
growing steadily, allowing greater interoperability, and becoming an important
tool in numerical exploration of physical phenomena, particularly in soft
matter systems. Driven by the need for fast and precise numerical integration
in colloidal dynamics, here we formulate the problem of Brownian Dynamics (BD)
in a mathematically consistent formalism of the Itō calculus, and develop a
Python package to assist numerical computations. We show that, thanks to the
automatic differentiation packages, the classical truncated Taylor-Itō
integrators can be implemented without the burden of computing the derivatives
of the coefficient functions beforehand. Furthermore, we show how to
circumvent the difficulties of BD simulations such as calculations of the
divergence of the mobility tensor in the diffusion equation and discontinuous
trajectories encountered when working with dynamics on $S^{2}$ and $SO(3)$.
The resulting Python package, Pychastic, is capable of performing BD
simulations including hydrodynamic interactions at speeds comparable to
dedicated implementations in lower-level programming languages, but with a
much simpler end-user interface.

###### Contents

  1. 1 Introduction
  2. 2 Three vantage points: Langevin, Fokker-Planck, and Itō
  3. 3 Pychastic: description of the package
    1. 3.1 Available numerical integration packages
    2. 3.2 Implementation details
  4. 4 Examples of usage
    1. 4.1 First passage problems, polar random walk
    2. 4.2 Rotational Brownian motion, step_post_processing function
    3. 4.3 Bead models with hydrodynamic interactions, pygrpy package integration
  5. 5 Conclusion
  6. A Appendix: Typos in integration schemes in Ref. [1]
  7. B Appendix: Rotational Brownian motion

##  1 Introduction

The dynamics of soft and colloidal matter systems is of importance for
numerous technological and industrial processes, such as food products,
pastes, creams, and gels [2]. Another important example are biological systems
which involve aqueous suspensions of colloids, macromolecules, polymers, and
cells. The diversity of the constituent elements, together with the tunability
of their direct interactions, and the presence of hydrodynamic interactions
(HI) mediated by the suspending fluid [3], gives rise to a multitude of
complex dynamic phenomena, which can only be explored using appropriate
numerical techniques [4, 5].

The choice of a suitable method derives from the characteristic time scales of
the investigated dynamics [6]. In the context of macromolecules, when short-
time effects are of interest, Molecular Dynamics (MD) simulations are a
popular choice [7]. Short times may be comparable to the solvent relaxation
time scale, which for a fluid of kinematic viscosity $\nu$ and speed of sound
$c$ scales as $\nu/c^{2}$. For water, this time is of the order of $10^{-13}$
s. The idea behind MD simulations is to solve Newtonian equations of motion
for atoms or molecules, which are a set of second-order ODEs with effective
interaction potentials. In this case, HI can be resolved using either explicit
solvent methods that resolve the molecular structure of the solvent or more
approximate implicit solvent models. However, the typical time scales of
colloidal motion are much longer; the velocity relaxation for colloidal
particles of size $a$ and density comparable to that of the fluid is
approximately $a^{2}/\nu$, amounting to ca. $10^{-8}$ s for a 100 nm colloid.
On coarser time scales, the colloidal velocity relaxes multiple times, and the
Rayleigh description in terms of velocity becomes irrelevant [6, 8]. Instead,
the description in terms of the position of a Brownian particle becomes
possible on diffusive time scales, during which a particle diffuses over a
distance of its own radius. This time scale is $a^{2}/D\sim\eta a^{3}/k_{B}T$,
where $\eta$ is the dynamic viscosity of the fluid, $k_{B}$ is the Boltzmann
constant, and $T$ is the temperature. For the aforementioned particle, this
time scale is ca. $10^{-3}$ s, and corresponds to the minimal time resolution
of typical light scattering and microscopy experiments. Thus, the description
of the dynamics on such coarse scales is made solely using the position of the
particle. The clear separation of discussed time scales, combined with the
numerical stiffness of Newton’s equations, renders MD simulations not
applicable to explore the diffusive dynamics.

Instead, BD simulations present a convenient alternative, building on models
suitable for diffusive time scales, and involving accordingly long time steps
for numerical computations. Since the velocity of particles varies very
rapidly on such time scales, the position becomes a random variable whose
properties are dictated by the fluctuation-dissipation theorem [8].

The development of a suitable theoretical description to explain the Brownian
motion and its subsequent experimental verification at the beginning of the
20th century were a groundbreaking step that confirmed the atomic structure of
matter [9]. Early works by pioneers such as Einstein [10], Smoluchowski [11],
Fokker [12], and Planck [13], led to the description of diffusive processes
through the underlying probability density functions (PDFs) for the positions
of the particles and their temporal evolution, rather than representations of
individual trajectories. The realisation that the path of a Brownian particle
resembles a nowhere differentiable function paved the route for a different
approach, now called stochastic calculus, with the classical works of
Kolmogorov [14], Wiener, and finally Itō [15]. With these tools, the problem
of Brownian motion can be recast in the form of stochastic differential
equations (SDEs), which offer a practical route for numerical simulations.

In this contribution, we present our approach to implementation of such
Brownian dynamics. Our goal was to show how modern programming techniques
available in Python: reflective and functional programming paradigms allow one
to develop concise implementations of physically relevant problems. These
constitute relatively recent developments in computer modelling, and we follow
this trend to facilitate faster modelling of diffusive phenomena.

First, in Sec. 2 we review three popular ways of describing the Brownian
motion using the Langevin equation, the Fokker-Planck equation, and Itō
equation. We then provide an overview of available numerical integration
packages for SDEs in Sec. 3.1. Next, we present our new package, Pychastic, in
Sec. 3, along with three examples of problems in which the packages greatly
facilitate numerical computations. We conclude the work in Sec. 5.

##  2 Three vantage points: Langevin, Fokker-Planck, and Itō

Perhaps the most popular approach to modelling the trajectories of Brownian
particles dates back to Langevin’s work of 1908 [16, 17]. The Langevin
equation can be rationalised as an extension of Newtonian mechanics to include
the effects of fluctuations by adding a stochastic force, $F_{n}$, acting on a
particle. We write it below for a particle with a constant friction
coefficient $\zeta$. The average value of the force is zero, while its
covariance at temperature $T$ is given by the fluctuation-dissipation theorem
[8], $\langle
F_{n}(t)F_{n}(t^{\prime})\rangle=2k_{B}T\zeta\delta(t-t^{\prime})$, with
$k_{B}$ being the Boltzmann constant. The equation of motion of the particle
in the presence of a deterministic force $F$ reads then

$m\ddot{x}=-\zeta\dot{x}+F(x,t)+F_{n}(t),$ (1)

For the case of a spherical particle of radius $a$ in Stokes flow the friction
coefficient is $6\pi\eta a$, $\eta$ being the dynamic viscosity. We note that
even though $x(t)$ is not differentiable once, and certainly not
differentiable twice at any point even in the usual distributional sense
(because the Lebesgue integral requires finite variation and realisations of
the Wiener process have infinite variation), Eq. (1) can be given a proper
interpretation by transforming it into an integral form. However, the above
equation proved practical for numerical calculations, e.g. using an
integration scheme analogous to the Euler-Maruyama method.

An alternative description via the Fokker-Planck equation circumvents the
difficulties of interpreting the dynamics on a single trajectory in the
Langevin equation by focussing on a probabilistic description. On the Brownian
time scale, the position of a particle is now treated as a random variable,
and the probability density function (PDF) $P(x,t)$ of finding the particle at
a location $x$ at a time $t$ evolves according to the partial differential
equation

$\frac{\partial P(x,t)}{\partial t}=\zeta^{-1}\frac{\partial}{\partial
x}\left(F(x,t)P(x,t)\right)+D\frac{\partial^{2}P(x,t)}{\partial x^{2}},$ (2)

with the diffusion coefficient defined by $D=k_{B}T/\zeta$. The problem is now
well-posed mathematically, and the tools of mathematical analysis of PDEs can
now be employed to study the evolution of the underlying probability
distribution, e.g. by examining its moments. However, questions involving
individual trajectories become less straightforward: problems involving first
passage time have to be dealt with with a careful treatment of the boundary
conditions (e.g. substances vanishing in chemical reactions). Moreover, from
the point of view of numerical simulations, this approach becomes impractical
for large systems. Using finite-difference methods, for $N$ particles in a 3D
simulation box of size $L$ and mesh size $\Delta x$, we typically require
$(L/\Delta x)^{3N}$ points to track the PDF, which may become prohibitively
large.

The mathematical difficulties of the Langevin equation are absent in the
proper treatment of SDEs in the Itō formalism. If now $dX$ denotes the
position increment of a particle in a time interval $dt$, we can write it as

$dX=Fdt+\sqrt{2D}dW,$ (3)

where $Fdt$ denotes the systematic drift of the particle and $\sqrt{2D}dW$
denotes the stochastic (diffusive) component of motion (provided that $D$ is
constant and the metric tensor is constant as well). This equation is meant in
the distributional sense with respect to the Itō integral, that is

$X(T)=\int_{0}^{T}dX=\int_{0}^{T}F(X,t)dt+\int_{0}^{T}\sqrt{2D}dW.$ (4)

Here, all integrals are taken in the Itō sense. Such a trajectory-focused
formulation effectively deals with the ill-defined derivatives in the Langevin
equation. Using Itō’s lemma imposes formal rules of transformation of the
coefficients, freeing us from a canonical coordinate description. Finally, the
estimation of observable quantities such as expected values, correlation
functions or equilibrium distributions of low-dimensional projections of
evolving variables can be recovered using a Monte Carlo approach, which for
$M$ simulations converges as $\sqrt{M}$, regardless of dimensionality. The
constant of proportionality of this convergence is controlled only by the
variance of the variable of interest, which is often independent of the
dimensionality of the equation. For example, in the problem of diffusion of a
single particle in a semidilute suspension, the variance of the mean squared
displacement of the tracer is weakly dependent on the number of particles in
the simulation volume.

##  3 Pychastic: description of the package

###  3.1 Available numerical integration packages

Only a few SDE integration packages have been made available in recent years.
Two notable examples are DifferentialEquations.jl for Julia and ItoProcess
being a part of Mathematica. Now, we present Pychastic for Python, which takes
advantage of the popularity of this language. Our package source code is
available on Github, up-to-date documentation on ReadTheDocs and ready to
install via pip via Python Package Index. As a preliminary comparison, we note
that DifferentialEquations.jl has the largest variety of integrators (for
example, many options for stiff equations), while Mathematica’s ItoProcess and
Pychastic contain essentially the same algorithms. However, the biggest
drawback of ItoProcess is the lack of step post-processing. When working with
SDEs defined on manifolds whose universal cover is not $\mathbb{R}^{n}$, such
as a sphere $S^{2}$ or the space of rigid rotations $SO(3)$, any
parameterisation of $\mathbb{R}^{n}$ will contain singularities, thus an
integrator which cannot handle discontinuous paths cannot reproduce, e.g. 3-
dimensional rotational dynamics [18]. Table 1 contains a synthetic comparison
of the three mentioned packages.

Package | DifferentialEquations.jl | ItoProcess | Pychastic
---|---|---|---
Language | Julia | Mathematica | Python
License | MIT | proprietary | MIT
Codebase | open | closed | open
scalar SDEs | yes | yes | yes
vector SDEs | yes | yes | yes
strong convergence | up to order 1.5 | up to order 1.5 | up to order 1.5
weak convergence | up to order 2.0 | up to order 2.0 | up to order 2.0
supports events | yes | no | yes
Table 1: Available SDE integration packages.

###  3.2 Implementation details

Pychastic contains implementations of three numerical integration schemes
based on the Taylor-Itō expansion. They are the schemes of strong order $1/2$,
$1$, and $3/2$ which in the package are referred to as Euler, Milstein, and
Wagner-Platen schemes, respectively.

The basis for these integration schemes is the Taylor-Itō expansion [1] which
generalises the deterministic Taylor expansion to SDEs. We write it here for a
one-dimensional problem driven by a Wiener process $W$

$dX=a(X)dt+b(X)dW.$ (5)

Following the notation of Kloeden & Platen [1], the Taylor-Itō expansion up to
strong order 3/2 for a scalar function $X$ has the form

$\begin{split}X(T)=&X(0)+a\Delta+b\Delta W\\\
&+\frac{1}{2}bb^{\prime}\left((\Delta W)^{2}-\Delta\right)\\\
&+a^{\prime}b\Delta
Z+\frac{1}{2}\left(aa^{\prime}+\frac{1}{2}b^{2}a^{\prime\prime}\right)\Delta^{2}\\\
&+(ab^{\prime}+\frac{1}{2}b^{2}b^{\prime\prime})\left(\Delta W\Delta-\Delta
Z\right)+\frac{1}{2}b\left(bb^{\prime\prime}+(b^{\prime})^{2}\right)\left(\frac{1}{3}(\Delta
W)^{2}-\Delta\right)\Delta W,\end{split}$ (6)

where

$\displaystyle\Delta$ $\displaystyle=\int_{0}^{T}dt=T,$ (7)
$\displaystyle\Delta W$ $\displaystyle=\int_{0}^{T}dW=W(T)-W(0),$ (8)
$\displaystyle\Delta Z$ $\displaystyle=\int_{0}^{T}\int_{0}^{s}dW(s)\,ds.$ (9)

The stochastic Euler scheme is based on the first line of this expansion, with
three terms only. The Milstein scheme includes the next term, namely the
second line of Eq. (6). Finally, the Wagner-Platen scheme includes all the
terms mentioned above. We note from this expansion that even the simplest
schemes require correct computation of coefficients multiplying principal
Wiener integrals ($\Delta$, $\Delta W$, $\Delta Z$, and others in the
multidimensional case) and sampling from correct distributions corresponding
to these integrals. Importantly, an analogous expansion for vector quantities
is considerably more complex, and other Wiener integrals arise, as detailed in
[1]. For brevity, we do not write this expansion here explicitly, but we have
implemented the vector Taylor-Itō expansion in Pychastic to enable simulations
of both scalar and vector processes.

Although expressions for these coefficient functions can be written explicitly
in principle, in Pychastic we take advantage of functional programming tools,
which results in a greatly simplified implementation. First, by introducing
the $\mathcal{L}^{0}$ and $\mathcal{L}^{j}$ operators (again using the Kloeden
& Platen notation), we can express all coefficient functions by repeated
application of the $\mathcal{L}$ operators to the $a$ and $b$ functions. Using
the jax.grad functionality, this is implemented directly as

⬇

1def tensordot1(a, b):

2 return jax.numpy.tensordot(a, b, axes=1)

3

4def tensordot2(a, b):

5 return jax.numpy.tensordot(a, b, axes=2)

6

7# Taylor-Ito expansion operators

8def L_t_operator(f,problem):

9 @wraps(f)

10 def wrapped(x):

11 b_val = problem.b(x)

12 val = tensordot1(jax.jacobian(f)(x), problem.a(x)) + 0 .5 * tensordot2(

13 jax.hessian(f)(x), tensordot1(b_val, b_val.T)

14 )

15 return val[:,jnp.newaxis,...] #indexing convention [spatial, time, ... =
noiseterms/time]

16

17 return wrapped

18

19def L_w_operator(f,problem):

20 @wraps(f)

21 def wrapped(x):

22 val = tensordot1(jax.jacobian(f)(x), problem.b(x))[:,jnp.newaxis,...]

23 return jnp.swapaxes(val,1,-1)[:,..., 0 ] # indexing convention [spatial,
noiseterms, ... = noiseterms/time]

24

25 return wrapped

Second, different integration methods – Euler, Milstein, and Wagner-Platen –
were implemented following the book by Kloeden & Platen [1]. These methods
make use of samples of the principal Wiener integrals listed above.
Unfortunately, the text contains typographic errors, which were found by
examining the mean, variance and covariance properties of fundamental multiple
Itō and Stratonovich integrals using a testing suite that is part of the
Pychastic package. The errors we found were corrected appropriately and are
listed in the Appendix A.

Finally, using the jax.lax.scan functionality of the package jax, we can
generate many trajectories in parallel by programmatically vectorising the
code representing $a$ and $b$ provided by the user.

##  4 Examples of usage

In the following, we show a few examples of usage that take advantage of
various functionalities of Pychastic.

###  4.1 First passage problems, polar random walk

In nearly all physical applications of Brownian Dynamics, the simulated
properties of the systems are observed via their moments rather than
individual realisations of the underlying stochastic process. In fact, in many
situations with random forcing, the paths of the process are only a model of
reality and cannot be directly observed, as opposed to their statistical
effect. For a smooth quantity $g$ observed through its expected value
$\mathbb{E}$ in thermodynamic equilibrium, it is the weak convergence rate
that controls the error of $\mathbb{E}(g(X_{t}))$ at fixed $t$. With this in
mind, it is tempting to dismiss methods with a high strong order of
convergence as impractical.

A natural candidate for a counterexample are first passage times, where the
answer requires a more subtle reasoning. Since the first passage problem is
_path-dependent_ , it would seem that strong convergence is important.
However, the theorem regarding the first passage times established by Whitt
[19] shows that the weak convergence of the process approximations $X_{n}\to
X$ implies a weak convergence of the first passage times $T(X_{n})\to T(X)$.
In conclusion, for all popular physical applications, the weak order of
convergence is the important one. In consequence, the Milstein scheme is never
a good choice, since it is equivalent to the Euler scheme in terms of weak
order of convergence, but is more computationally intensive and harder to
implement. We illustrate the surprising result of [19] by numerical simulation
in a familiar setting.

A simple analytically solvable case of the first passage problem is a two-
dimensional diffusion process in the $(X,Y)$ plane with constant drift
velocity $\bm{v}$ and constant diffusion coefficient $\sigma^{2}$ in a system
with an absorbing barrier at $Y=Y_{\text{b}}=2$. The problem admits an exact
solution, in which the probability density function for the hitting time
$t_{\text{hit}}=\mathrm{argmin}_{t}(Y_{t}>Y_{\text{b}})$ is given by

$p_{t_{\text{hit}}}(\tau)=\frac{Y_{\text{b}}}{\sqrt{2\pi\sigma^{2}\tau^{3}}}\exp\left(-\frac{(Y_{\text{b}}-v_{y}\tau)^{2}}{2\sigma^{2}\tau}\right),$
(10)

and so $\mathbb{E}[t_{\text{hit}}]=1/2$.

Numerically, we calculate the sample of the first passage times and locations
for $\bm{v}=(v_{x},v_{y})=(0,4)$ with stopping condition $Y=Y_{\text{b}}=2$
and initial condition $(X,Y)=(2,0)$. The Itō equations for the two-dimensional
diffusion driven by Wiener processes $W$ and $W^{\prime}$ with drift
velocities $[v_{x},v_{y}]$ take the form of

$\displaystyle dX=$ $\displaystyle v_{x}dt$ $\displaystyle+\sigma dW_{t},$
(11) $\displaystyle dY=$ $\displaystyle v_{y}dt$ $\displaystyle+\sigma
dW^{\prime}_{t},$ (12)

where $\sigma$ is a parameter that describes the strength of the noise
compared to the drift. We first transform this problem into a different set of
coordinates, for which we choose polar coordinates. While we retain the
analytical solution, we acquire non-linear terms in the governing equations,
which will help us study the convergence issues which arise when noise and
drift terms are generally dependent on the position. Eq. (12) transformed into
polar coordinates $(r,\phi)$ is expressed as

$\displaystyle dr=$
$\displaystyle\left(\frac{1}{2r}+v_{y}\sin\phi+v_{x}\cos\phi\right)dt$
$\displaystyle+\sigma\cos\phi dW_{t}+\sigma\sin\phi dW^{\prime}_{t},$ (13)
$\displaystyle d\phi=$
$\displaystyle\frac{1}{r}\left({v_{y}\cos\phi+v_{x}\sin\phi}\right)dt$
$\displaystyle-\sigma\frac{\sin\phi}{r}dW_{t}+\sigma\frac{\cos\phi}{r}dW^{\prime}_{t}.$
(14)

Figure 1: a) Four sample trajectories generated from Eq. (14) using the Euler
algorithm together with the barrier $r\sin\phi=y=2$. b) Distribution of
hitting times for two different step sizes in the Euler algorithm together
with the analytic solution of the problem. Large step sizes typically
overestimate the hitting time. c) Strong (path-wise) error of the hitting time
for different step sizes and solvers together with power-law eye guides:
$dt^{1/4}$, $dt^{1/2}$, $dt$. d) Error of the expected value of the first
passage time together with eye-guides: $dt^{1/2}$, $dt$. Contrary to the
strong error case, here Euler and Milstein algorithms coincide by Whitt’s
theorem.

The transformed equations are nonlinear, coupled, with nondiagonal,
noncommutative noise terms. Vector equations in the package Pychastic are
defined by providing both the drift vector a and the noise matrix b as
callable to the SDEProblem constructor. Denoting the configuration by
$\bm{q}=(r,\phi)$, we write

⬇

1import jaq.numpy as jnp

2

3y_drift = 4. 0

4sigma = 1

5y_barrier = 2

6

7def drift(q):

8 return jnp.array(

9 [1 / (2 * q[ 0 ]) + y_drift * jnp.sin(q[1]), y_drift * jnp.cos(q[1]) / q[ 0
]]

10 )

11

12def noise(q):

13 return sigma * jnp.array(

14 [[jnp.cos(q[1]), jnp.sin(q[1])], [-jnp.sin(q[1]) / q[ 0 ], jnp.cos(q[1]) /
q[ 0 ]]]

15 )

16

17problem = pychastic.sde_problem.SDEProblem(

18 a=drift,

19 b=noise,

20 x0=jnp.array([2. 0 , 0 . 0 ]),

21 tmax=2. 0 ,

22)

The jax.numpy package is a functional cousin of numpy [20], which is the
fundamental package for scientific computing with Python [21]. The functional
focus of jax.numpy enables automatic differentiation, and thus facilitates the
use of the high-order method without the need for numerical differentiation,
which can be imprecise and computationally costly. We obtain the trajectories
using the solve method of SDESolver.

⬇

1solver = pychastic.sde_solver.SDESolver(dt=2**(-5),
scheme=method=hex,ActualText=27’milsteinmethod=hex,ActualText=27’)

2solution = solver.solve_many(problem,n_samples,seed= 0 )

Pychastic supports simulating many trajectories simultaneously with the
solve_many method. Our strategy allows for concurrency in computation of
random variables and equation coefficients thanks to the jax package we use as
back-end (which in turn relies on vectorised mathematical operations provided
by the XLA supporting architecture of most modern processors). The jax package
is a fusion of three capabilities: numpy-like API for array-based computing,
functional transformations (such as vectorisation, parallelisation and
automatic differentiation), and modular back-end, allowing developers to test
their work with just CPU and deploy the same code later on a GPU (or TPU)
capable hardware. Furthermore, we avoid notoriously slow Python loops by using
jax.lax.scan routines, taking advantage of just-in-time compilation and
asynchronous dispatch, avoiding the problem of global interpreter lock, which
cripples imperatively coded Python programmes. A more detailed guide on the
advantages (and common pifalls) of jax package can be found in its
documentation on ReadTheDocs.

Importantly, even though the Milstein and Wagner-Platen schemes require the
values of spatial derivatives of noise and drift terms, we did not have to
provide them explicitly. They were calculated using an automatic
differentiation procedure from the coefficient functions provided in the
SDEProblem constructor.

We compare the results obtained with different solvers and time steps in Fig.
1. Sample trajectories starting from the point $(r,\phi)=(2,0)$ are shown in
Fig. 1a. The barrier at $y=2$ becomes a curve in polar parameterisation. In
Fig. 1 we present the resulting distribution of hitting times measured from an
ensemble of $N=1000$ trajectories with two different time-step sizes and
compared to the exact solutions. Choosing a too large step size can lead to an
overestimation of the typical hitting time. Having implemented three
stochastic integration algorithms, we compare the strong error of estimation
of the hitting time in Fig. 1c. It is a measure of the numerical error per
trajectory, defined as $\langle|\tilde{\tau}-\tau|\rangle$, where
$\tilde{\tau}$ is the numerical estimate of hitting time and $\tau$ is actual
hitting time that can be computed using the exact solution and a particular
realisation of the Wiener process. The average is taken over an ensemble of
realisations of the Wiener process. For different schemes, the strong error
scales differently with the step size – from the linear dependence on
$dt^{1/4}$ for the Euler scheme, through $dt^{1/2}$ for the Milstein
algorithm, to $dt$ for Wagner-Platen. We note that these exponents are
different (smaller) than for strong convergence at a fixed time [1].
Similarly, in Fig. 1d we present the weak error of estimation of the average
value of the hitting time. The weak error is defined as
$|\langle\tilde{\tau}\rangle-\langle\tau\rangle|$, and reflects the difference
between the ensemble average value of hitting time estimate and the true
ensemble average hitting time (computed directly from the theoretical
distribution). This quantity highlights the lack of difference between Euler
and Milstein schemes, in agreement with Whitt’s theorem [19] stating the
dependence of the expected first passage time on the weak convergence rate of
the scheme only. The power laws for convergence – $dt^{1/4}$ for Euler and
Milstein algorithms and $dt$ for the Wagner-Platen scheme – again have smaller
exponents than for those for weak convergence at a fixed time.

###  4.2 Rotational Brownian motion, step_post_processing function

To highlight the ease-of-use features of the Pychastic package, we present the
problem of rotational Brownian diffusion simulations. General BD simulations
that involve rigid bodies require the application of finite rotations to
diffusing objects. The resulting algorithms are widely applied to study
problems in physics and biology. While resolving translational motion is
straightforward, e.g. with the standard Ermak-McCammon (Euler) algorithm [22],
the rotational part is more involved to simulate, because the domain of
rotational motion is $SO(3)$ and this has to be taken into account when
solving the equations of motion. One difficulty lies in the commonly used
rotational coordinate systems, such as Euler angles [23], which contain strong
singularities around the polar orientations. When curvilinear coordinates are
used to describe the motion, the metric tensor gives rise to nontrivial
additional terms in the equations of motion, the so-called metric or drift
terms, which are frequently overlooked in rotational BD algorithms. For a
summary, see Ref. [24].

To overcome these limitations, we reimplement the problem following the
rotation-vector-based formulation of Evensen et al. [25]. We describe the
angular position of the particle by the angle of rotation $\Phi$ around a unit
vector $\bm{\delta}$ collinear with the axis of rotation. Generalised
coordinates are encoded in the vector $\bm{q}=\Phi\bm{\delta}$. The rotational
mobility matrix of the particle in the body-fixed frame is
$\bm{\mu}_{\text{body}}$. We transform it into the laboratory frame by

$\bm{\mu}=\bm{\Omega}^{T}\cdot\bm{\mu}_{\text{body}}\cdot\bm{\Omega},$ (15)

where $\bm{\Omega}$ is the relevant rotation matrix. We also define a velocity
transformation matrix $\bm{\Xi}$ to further introduce the transformed mobility
$\widehat{\bm{\mu}}$ by

$\widehat{\bm{\mu}}=\bm{\Xi}\cdot\bm{\mu}\cdot\bm{\Xi}^{T}.$ (16)

Now we can write the Itō SDE corresponding to the evolution of the generalised
coordinates as

$d\bm{q}=\widehat{\bm{\mu}}\cdot\left(\frac{\partial}{\partial\bm{q}}\log
V\right)dt+k_{B}T\left(\frac{\partial}{\partial\bm{q}}\cdot\widehat{\bm{\mu}}\right)dt+\sqrt{2k_{B}T}\widehat{\bm{\mu}}^{1/2}\cdot\,\,d\bm{W}.$
(17)

Here, $V$ is the density of the volume element ($SO(3)$ Haar measure) with
respect to the Lebesgue measure on $\mathbb{R}^{3}$. In the literature
${\partial\log V/\partial\bm{q}}$ is often called the metric force
$\bm{F}^{(m)}$. Explicit expressions for $\bm{\Xi}$, $\bm{\Omega}$ and $\log
V$ are given in Appendix B.

a) b)

Figure 2:  a) Distribution of rotation angles $\Phi$ after a long time in
$N=10^{4}$ simulations compared with the $SO(3)$ symmetric distribution
$P^{eq}(\Phi)$ given by Eq. (18). b) Comparison of the time-dependent
correlation structure arising from equation (17) for a spherical particle.
Dashed lines show exact predictions from Ref. [26]. Solid, colored lines show
average values of $10^{3}$ simulations (blue line corresponds to equation (17)
and orange points correspond to Ref. [25] containing a typo in their Eq. 12).
Note that even though equilibrium distributions coincide, time dependent
correlation structure is different.

Since the evolution equation depends on the divergence of a product of two
orientation-dependent matrices, writing Eq. (17) explicitly is quite
cumbersome. Earlier works of [25] and [27] avoid this difficulty by
approximating the gradients by sampling the transformed mobility matrix at
nearby locations in the phase space. We can address this problematic term
directly thanks to the automatic differentiation capabilities of the jax
package.

⬇

1def metric_force(q):

2 # Metric force, providing the Boltzmann distribution in equilibrium.

3 phi = jnp.sqrt(jnp.sum(q ** 2))

4 scale = jax.lax.cond( # Taylor expansion for the polar orientations.

5 phi < 0 . 0 1,

6 lambda t: -t / 6. 0 ,

7 lambda t: jnp.sin(t) / (1. 0 - jnp.cos(t)) - 2. 0 / t,

8 phi,

9 )

10 return jax.lax.cond(

11 phi > 0 . 0 , lambda: (q / phi) * scale, lambda: jnp.array([ 0 . 0 , 0 . 0
, 0 . 0 ])

12 )

13

14def t_mobility(q):

15 # Mobility matrix transformed to rotation vector coordinates.

16 Return transformation_matrix(q) @ mobility @ (transformation_matrix(q).T)

17

18def drift(q):

19 # jax.jacobian has the differentiation index last (like mu_ij d_k) so
divergence is contraction of the first and last axis.

20 return (

21 t_mobility(q) @ metric_force(q)

22 + jnp.einsum("iji->j", jax.jacobian(t_mobility)(q))

23 )

First, including the divergence of the mobility matrix in the new coordinate
system is as simple as adding a correct contraction of jax.jacobian of the
mobility tensor. To obtain this, we use the convenient jnp.einsum tool. This
function, similarly to its analogue np.einsum, takes in a tensor, represented
in the memory by an array of arrays of arrays etc. with the number of levels
corresponding to the tensor rank, and performs contraction described
symbolically by the index notation. In our case "iji->j" means ‘contract the
first index with the last index of a rank-3 tensor to produce a vector.’

Second, the domain of the equation is $SO(3)$, which is not covered by
$\mathbb{R}^{3}$ (since $\pi_{1}(SO(3))=\mathbb{Z}_{2}$ [18]). As a result,
some trajectories, continuous in $SO(3)$, will necessarily be discontinuous in
the $\mathbb{R}^{3}$ parameterisation. We adopt the method proposed by Ref.
[25], where after each step of integration we project on the principal value
of the rotation angle. In pychastic, we achieve this using the
step_post_processing capability.

⬇

1def canonicalize_coordinates(q):

2 phi = jnp.sqrt(jnp.sum(q ** 2))

3 max_phi = jnp.pi

4 canonical_phi = jnp.fmod(phi + max_phi, 2. 0 * max_phi) - max_phi

5 return jax.lax.cond(

6 phi > max_phi,

7 lambda canonical_phi, phi, q: (canonical_phi / phi) * q,

8 lambda canonical_phi, phi, q: q,

9 canonical_phi,

10 phi,

11 q,

12 )

13solver = pychastic.sde_solver.SDESolver(dt= 0 . 0 1)

14trajectories = solver.solve_many(

15 problem,

16 step_post_processing=canonicalize_coordinates,

17 n_trajectories=1 0 0 0

18)

The simulations were validated by checking the equilibrium properties and the
time-dependent correlation structure. The equilibrium distribution
$P^{eq}(\Phi)$ for the rotation angle $\Phi$ is given by

$\displaystyle P^{eq}(\Phi)=\frac{1-\cos\Phi}{\pi}.$ (18)

To characterise the time evolution, we use the components of the modified
rotation vector components $\Delta
u_{k}(t)=-\frac{1}{2}\epsilon_{ijk}\Omega_{ij}(t)$, for which exact
theoretical predictions were derived by Cichocki et al. [26] for an
arbitrarily shaped molecule. For a spherical particle, the correlation
functions are given by

$\displaystyle\langle\Delta u_{k}(t)\Delta
u_{l}(t)\rangle_{0}=\left[\frac{1}{6}-\frac{5}{12}e^{-6D_{r}t}+\frac{1}{4}e^{-2D_{r}t}\right]\delta^{K}_{kl},$
(19)

where $\delta^{K}$ is the Kronecker delta. The diffusion coefficient for a
sphere is $D_{r}=k_{B}T/\pi\eta d^{3}$, where $d$ is the diameter of the
particle.

In Fig. 2, we present a comparison between the results obtained with an
algorithm based on Eq. (17) (blue lines), and a similar method based on
equations (5-12) from Ref. [25] (orange). Fig. 2a shows the equilibrium
distribution of the rotation angle $\Phi$, together with the analytical
prediction of Eq. (18). We note that although the equilibrium distribution
predicted by both algorithms shows exact agreement with the theoretical
predictions, the transformation matrix $\bm{\Xi}$ in Eq. (12) of Ref. [25]
contains an error. We corrected this typographical error and present the
proper formulation of the equations of motion matrices in the Appendix B. The
discrepancy is visible in the time-dependent correlation function in Fig. 2b.
Although at long times the system tends to the same equilibrium solution for
both approaches, we see the agreement of our algorithm with the theoretical
predictions of Eq. (19) (dashed line) and the deviation of the orange lines.
This discrepancy highlights the need for better test cases that reliably test
all properties of the simulated equation, and not solely the equilibrium
distribution.

As mentioned above, the algorithm based on Ref. [25] is singularity-free,
contrary to approaches based on Euler angles, which contain strong
singularities around the polar orientations. However, it should be noted that
the implementation of the code may still face some limitations, such as the
computational inability to calculate $\sin\Phi/\Phi$ for $\Phi\rightarrow 0$.
For these numerically restricted cases, we implemented the Taylor expansion of
$\bm{\Xi}$, $\bm{\Omega}$ and $\log V$, based on Ref. [27]. However, we also
report some typos in this publication. Therefore, in Appendix B we provide
correct formulations of the Taylor-expanded terms.

###  4.3 Bead models with hydrodynamic interactions, pygrpy package
integration

Hydrodynamically interacting beads (with and without springs) have been used
successfully in modelling the properties of elastic macromolecules [28, 22,
29, 30, 31, 32]. Some questions about elastic macromolecules can be answered
by computing the equilibrium ensemble of conformations (with methods such as
Markov Chain Monte Carlo [33]). Simulations involving hydrodynamic
interactions can, on the other hand, provide access to the dynamics and answer
questions, e.g., about timescales of conformational change [34], mechanisms of
protein association [35], pore translocation [36, 37, 38], near-wall hindered
diffusion [39], and many other dynamical processes.

The starting point is often the Hamiltonian, which describes intramolecular
interactions between the constitutive subunits of the molecule, modelled by a
collection of beads. The potential energy landscape can then be used to
compute the interaction forces. Thanks to the jax autograd capabilities, the
forces arising from many-body mechanical interactions can be automatically
calculated from the potential energy of the system.

The computation of hydrodynamic interactions, encoded in the mobility tensors
of the respective macromolecules is a more involved task. Numerous methods are
available, varying in scope, degree of precision, and complexity. For a review
of popular methods applied to macromolecules, see Ref. [40]. For completeness,
we supplemented pychastic with the package pygrpy for the calculation of grand
mobility tensors in the Rotne-Prager-Yamakawa (RPY) approximation [41, 42],
generalised to beads of different sizes [43]. The procedure is a Python port
of the GRPY package [40]. The Python package pygrpy simplifies the
implementation of similar bead-spring simulations, both stochastic and
deterministic. Furthermore, pygrpy is compatible with the jax functional
paradigm and allows automatic differentiation and vectorisation.

The RPY approximation is by far the most popular method of accounting for
hydrodynamic interaction in numerical models of soft matter systems [44]. The
mobility tensors calculated in this way preserve positive definiteness and are
divergence-free, significantly simplifying the BD algorithm [22]. In essence,
RPY is a far-field approximation that includes all terms that decay slower
than the inverse third power of the interparticle distances, but is less
accurate at smaller distances. When the particles come close together, it is
necessary to include higher-order terms of the multipole expansion [45] and
lubrication corrections [46]. However, as shown by Żuk et al. [43], the RPY
approximation can be generalised to overlapping particles, and thus it can
also be used to model complex-shaped particles as conglomerates of rigidly
glued overlapping spheres, in particular to calculate hydrodynamic properties
of biological macromolecules, as in the GRPY method [40].

To highlight the ease of use and interoperability of pychastic and pygrpy, we
implement a benchmark problem proposed by Cichocki et al. [47]. It concerns
the diffusion coefficient of an elastic "macromolecule" composed of 4 beads of
radii $r_{i}\in\\{3,1,1,1\\}$ joined into a string, as shown in Fig. 3a. Here,
the length scale is the radius of a small bead $a$. The neighbouring beads
interact directly via harmonic potentials, with the equilibrium distance
$d_{i}=4$ and the spring constant $k=5.5$ $k_{B}T/a$, and indirectly through
hydrodynamic interactions.

a) b)

Figure 3:  a) Representative configurations of four co-diffusing beads
connected with harmonic springs at time (bottom to top)
$t/\tau_{d}=2,\,2000,\,4000,\,6000$. b) Mean square displacement of three
different tracked locations: the small bead at the end of the chain (blue),
the large bead at the other end (orange) and the weighted average of all four
beads with weights corresponding to the effective centre of diffusion, defined
according to Ref. [47]. The dashed line corresponds to numerical results
reported in [47] based on an extremely long, single-trajectory simulation.

We determine the diffusion coefficient $D$ by tracking the mean square
displacement $\langle(\bm{x}(t)-\bm{x}(0))^{2}\rangle$ (MSD) of a point
$\bm{x}$ of the molecule. If the observation time is long enough,
$\text{MSD}\approx 6Dt$. However, for shorter times, the coefficient of
proportionality in the $\text{MSD}(t)$ curve is different and depends on the
choice of the reference point on the molecule. Although biophysical
experiments typically measure the long-time diffusion coefficient, numerical
simulations have easier access to short-time diffusivity. To minimise the
difference, one should choose a particular point, called the centre of
diffusion, which can be constructed as a weighted average of the positions of
beads with weights determined from the hydrodynamic mobilities of an ensemble
of equilibrium configurations of the elastic molecule [47].

We use this system to test our algorithm. In Fig. 3b, we show the temporal
evolution of the MSD when different reference points on the elastic molecule
are chosen: the small terminal bead, the large bead, and the centre of
diffusion. The time is scaled by $\tau_{d}=\pi\eta d^{3}/k_{B}T$. We computed
hydrodynamic interactions using the package pygrpy, which allows computations
of mobility tensors for macromolecules composed of unequally sized and
possibly overlapping spherical beads in the Rotne-Prager-Yamakawa
approximation. The dashed line is the numerical result of [47] where the
diffusion coefficient was estimated using a single very long BD trajectory.
Indeed, we see that tracking the centre of diffusion provides a good estimate
of the long-time behaviour of the MSD. When one chooses to track one of the
smaller beads or one of the beads further away from the middle of the
molecule, precision decreases primarily because rotational diffusion plays a
bigger role in the motion of such tracers. The strength of our algorithm is
that we use many trajectories instead of a single long one, thus avoiding the
need for a more complex method of calculating the diffusion coefficient,
involving recursive subdivisions of the simulation interval described by
Frenkel and Smit [48].

⬇

1radii = jnp.array([3. 0 ,1. 0 ,1. 0 ,1. 0 ]) # sizes of spheres used

2n_beads = len(radii)

3equilibrium_dist = 4. 0

4spring_constant = 5.5

5

6def pot_ene(x): # potential energy

7 locations = jnp.reshape(x,(n_beads,3))

8 distance_ab = jnp.sqrt(jnp.sum((locations[ 0 ] - locations[1])**2))

9 distance_bc = jnp.sqrt(jnp.sum((locations[1] - locations[2])**2))

10 distance_cd = jnp.sqrt(jnp.sum((locations[2] - locations[3])**2))

11 ene = 0 .5*spring_constant*jnp.sum(

12 ([distance_ab,distance_bc,distance_ca] - equilibrium_dist)**2

13 )

14 return ene

15

16def drift(x):

17 locations = jnp.reshape(x,(n_beads,3))

18 mu = pygrpy.jax_grpy_tensors.muTT(locations,radii)

19 force = -jax.grad(pot_ene)(x)

20 return jnp.matmul(mu,force)

21

22def noise(x):

23 locations = jnp.reshape(x,(n_beads,3))

24 mu = pygrpy.jax_grpy_tensors.muTT(locations,radii)

25 return jnp.sqrt(2)*jnp.linalg.cholesky(mu)

26

27problem = pychastic.sde_problem.SDEProblem(

28 drift,

29 noise,

30 x0 = jnp.reshape(jnp.array([

31 [-2., 0 ., 0 .],

32 [2., 0 ., 0 .],

33 [6., 0 ., 0 .],

34 [1 0 ., 0 ., 0 .]

35 ]),(3*n_beads,)),

36 tmax = 1 0 0 0 0 . 0 )

##  5 Conclusion

We developed a novel Python package pychastic dedicated to efficient numerical
solutions of SDEs. The package implements the classical truncated Taylor-Itō
integrators up to strong order $\mathcal{O}(dt^{3/2})$ providing a precise
treatment of both weak (e.g. equilibrium distributions, diffusion
coefficients) and strong problems. We included a set of simple test cases that
provide exact reference points for testing future stochastic integration
algorithms. The analysis of three-dimensional rotational Brownian motion
benchmarks is particularly important because it encompasses many of the
difficult aspects of SDE approaches to Brownian dynamics: divergence terms in
the evolution equation, handling discontinuous trajectories (unavoidable in
the case of $SO(3)$), and spurious agreements when testing only the
equilibrium distribution. We hope that pychastic will ease future studies of
Brownian dynamics problems, especially problems that involve hydrodynamic
interactions. The project is open source, and we hope to encourage
collaboration and its further development.

## Acknowledgements

The authors thank Piotr Szymczak for his insightful feedback.

#### Funding information

The work of ML, RW, and MB was supported by the National Science Centre of
Poland (FundRef DOI: http://dx.doi.org/10.13039/501100004281) grant Sonata to
ML no. 2018/31/D/ST3/02408.

##  Appendix A Appendix: Typos in integration schemes in Ref. [1]

In the book: _Numerical solution of stochastic differential equations_ by
Peter E. Kloeden and Eckhard Platen in equation 5.8.11 we have found the
following typos:

  1. 1.

In the definition of $D^{p}_{j1,j2,j3}$: the fourth summand should read
$\zeta_{j1,r}\eta_{j3,l+r}$ instead of $j1$.

  2. 2.

In the definition of $b_{j}$: there is a $1/\pi$ factor missing before sum.

  3. 3.

In the definition of $C^{p}_{j1,j2}$: second term should be
$+1/r\eta_{j1,r}\eta_{j2,l}$ instead of $-l/r$

##  Appendix B Appendix: Rotational Brownian motion

In equation (17) we recall the equations of rotational Brownian motion,
formulated in Ref. [25]. Angular orientation is described by the vector
$\bm{a}=\left(a_{1},a_{2},a_{3}\right)=\Phi\bm{\delta}$, where $\Phi$ is the
angle of rotation around a unit vector $\bm{\delta}$, which corresponds to the
axis of rotation. Because the expressions in Ref. [25] contain typos, here we
present the proper expressions for coordinate transformation and rotation
matrices $\bm{\Xi}$ and $\bm{\Omega}$:

$\displaystyle\bm{\Xi}=$
$\displaystyle\left(\frac{1}{\Phi^{2}}-\frac{\sin\Phi}{2\Phi(1-\cos\Phi)}\right)\begin{pmatrix}a_{1}a_{1}&a_{1}a_{2}&a_{1}a_{3}\\\
a_{2}a_{1}&a_{2}a_{2}&a_{2}a_{3}\\\
a_{3}a_{1}&a_{3}a_{2}&a_{3}a_{3}\end{pmatrix}+\frac{1}{2}\begin{pmatrix}\frac{\Phi\sin\Phi}{1-\cos\Phi}&-a_{3}&a_{2}\\\
a_{3}&\frac{\Phi\sin\Phi}{1-\cos\Phi}&-a_{1}\\\
-a_{2}&a_{1}&\frac{\Phi\sin\Phi}{1-\cos\Phi}\end{pmatrix},$ (20)
$\displaystyle\bm{\Omega}=$
$\displaystyle\frac{1}{\Phi^{2}}\begin{pmatrix}\Phi^{2}\cos\Phi&-a_{3}\Phi\sin\Phi&a_{2}\Phi\sin\Phi\\\
a_{3}\Phi\sin\Phi&\Phi^{2}\cos\Phi&-a_{1}\Phi\sin\Phi\\\
-a_{2}\Phi\sin\Phi&a_{1}\Phi\sin\Phi&\Phi^{2}\cos\Phi\end{pmatrix}+\frac{1-\cos\Phi}{\Phi^{2}}\begin{pmatrix}a_{1}a_{1}&a_{1}a_{2}&a_{1}a_{3}\\\
a_{2}a_{1}&a_{2}a_{2}&a_{2}a_{3}\\\
a_{3}a_{1}&a_{3}a_{2}&a_{3}a_{3}\end{pmatrix}.$ (21)

The matrix $\bm{\Xi}$ transforms the velocities from the Cartesian coordinate
system to the one described by $\bm{a}$. The matrix $\bm{\Omega}$ is a simple
rotation matrix. These two matrices combined allow transformation of the body-
fixed mobility matrix $\bm{\mu}_{\text{body}}$ to the lab-fixed,
$\bm{a}$-described mobility matrix, given as
$\widehat{\bm{\mu}}=\bm{\Xi}\cdot\bm{\Omega}^{T}\cdot\bm{\mu}_{\text{body}}\cdot\bm{\Omega}\cdot\bm{\Xi}^{T}$.

The term $\left(\frac{\partial}{\partial\bm{q}}\log V\right)$ from Eq. (17)
can be associated with the metric force $\bm{F}^{(m)}$, which guarantees the
Boltzmann distribution in equilibrium and is given as

$\left(\frac{\partial}{\partial\bm{q}}\log
V\right)=\bm{F}^{(m)}=-k_{B}T\left(\frac{\sin(\Phi)}{1-\cos(\Phi)}-\frac{2}{\Phi}\right)\bm{\delta}.$
(22)

Due to the numerical limitations, it is sometimes inevitable to perform the
Taylor expansion of the above quantities, as presented in Ref. [27]. However,
we also found typos in this case. Therefore, we provide the correctly expanded
matrices $\bm{\Xi}$ and $\bm{\Omega}$ for $\Phi\rightarrow 0$,

$\displaystyle\bm{\Xi}=$
$\displaystyle\frac{1}{12}\begin{pmatrix}a_{1}a_{1}&a_{1}a_{2}&a_{1}a_{3}\\\
a_{2}a_{1}&a_{2}a_{2}&a_{2}a_{3}\\\
a_{3}a_{1}&a_{3}a_{2}&a_{3}a_{3}\end{pmatrix}+\frac{1}{2}\begin{pmatrix}2&-a_{3}&a_{2}\\\
a_{3}&2&-a_{1}\\\ -a_{2}&a_{1}&2\end{pmatrix}+O(\Phi^{2}),$ (23)
$\displaystyle\bm{\Omega}=$ $\displaystyle\begin{pmatrix}1&-a_{3}&a_{2}\\\
a_{3}&1&-a_{1}\\\
-a_{2}&a_{1}&1\end{pmatrix}+\frac{1}{2}\begin{pmatrix}a_{1}a_{1}&a_{1}a_{2}&a_{1}a_{3}\\\
a_{2}a_{1}&a_{2}a_{2}&a_{2}a_{3}\\\
a_{3}a_{1}&a_{3}a_{2}&a_{3}a_{3}\end{pmatrix}+O(\Phi^{2}).$ (24)

The Taylor-expanded metric force around $\Phi\rightarrow 0$ becomes

$\bm{F}^{(m)}=k_{B}T\frac{\Phi}{6}\bm{\delta}+O(\Phi^{3}).$ (25)

## References

  * [1] P. Kloeden and E. Platen,  _Numerical Solution of Stochastic Differential Equations_ ,  Stochastic Modelling and Applied Probability. Springer Berlin Heidelberg,  ISBN 9783540540625 (2011).
  * [2] P.-G. De Gennes,  _Soft matter_ ,  Science 256(5056), 495 (1992).
  * [3] S. Kim and S. J. Karrila,  _Microhydrodynamics: principles and selected applications_ ,  Courier Corporation (2013).
  * [4] M. Praprotnik, L. D. Site and K. Kremer,  _Multiscale simulation of soft matter: From scale bridging to adaptive resolution_ ,  Annu. Rev. Phys. Chem. 59, 545 (2008).
  * [5] U. D. Schiller, T. Krüger and O. Henrich,  _Mesoscopic modelling and simulation of soft matter_ ,  Soft Matter 14, 9 (2018),  10.1039/C7SM01711A.
  * [6] J. K. G. Dhont,  _An introduction to dynamics of colloids_ , vol. 2 of _Studies in Interface Science_ ,  Elsevier,  ISBN 9780444820099 (1996).
  * [7] M. Karplus and G. A. Petsko,  _Molecular dynamics simulations in biology_ ,  Nature 347(6294), 631 (1990),  10.1038/347631a0.
  * [8] N. Van Kampen,  _Stochastic Processes in Physics and Chemistry_ ,  North-Holland Personal Library. Elsevier Science,  ISBN 9780080571386 (1992).
  * [9] R. M. Mazo,  _Brownian Motion: Fluctuations, Dynamics, and Applications_ ,  International Series of Monographs on Physics. Oxford University Press, Oxford,  ISBN 978-0-19-955644-1,  10.1093/acprof:oso/9780199556441.001.0001 (2008).
  * [10] A. Einstein,  _Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen_ ,  Annalen der Physik 322(8), 549 (1905),  https://doi.org/10.1002/andp.19053220806,  https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19053220806.
  * [11] M. von Smoluchowski,  _Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen_ ,  Annalen der Physik 326(14), 756 (1906),  https://doi.org/10.1002/andp.19063261405,  https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19063261405.
  * [12] A. D. Fokker,  _Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld_ ,  Annalen der Physik 348(5), 810 (1914),  https://doi.org/10.1002/andp.19143480507,  https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19143480507.
  * [13] M. Planck,  vol. Jan-Dec 1917,  Berlin, Deutsche Akademie der Wissenschaften zu Berlin, 1882-1918,  Https://www.biodiversitylibrary.org/bibliography/42231 (1917).
  * [14] A. Kolmogoroff,  _Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung_ ,  Mathematische Annalen 104, 415 (1931),  10.1007/BF01457949.
  * [15] K. Itô,  _On stochastic differential equations_ ,  4\. American Mathematical Soc. (1951).
  * [16] P. Langevin,  _Sur la théorie du mouvement brownien_ ,  C. R. Acad. Sci. (Paris) 146, 530–533 (1908).
  * [17] _Paul langevin’s 1908 paper “on the theory of brownian motion” [“sur la théorie du mouvement brownien,” c. r. acad. sci. (paris) 146, 530–533 (1908)]_ .
  * [18] B. C. Hall,  _Lie Groups, Lie Algebras, and Representations: An Elementary Introduction_ ,  Graduate Texts in Mathematics 222. Springer-Verlag New York, 1 edn.,  ISBN 0387401229,9780387401225 (2003).
  * [19] W. Whitt,  _Weak convergence of first passage time processes_ ,  J. Appl. Prob. 8(2), 417–422 (1971),  10.2307/3211913.
  * [20] Numpy contributors,  _Numpy homepage_ (2022).
  * [21] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus _et al._ ,  _Array programming with NumPy_ ,  Nature 585(7825), 357 (2020),  10.1038/s41586-020-2649-2.
  * [22] D. L. Ermak and J. A. McCammon,  _Brownian Dynamics with hydrodynamic interactions_ ,  J. Chem. Phys. 69(4), 1352 (1978),  10.1063/1.436761.
  * [23] H. Goldstein,  _Classical Mechanics_ ,  Pearson Education,  ISBN 9788177582833 (2002).
  * [24] I. M. Ilie, W. J. Briels and W. K. den Otter,  _An elementary singularity-free rotational brownian dynamics algorithm for anisotropic particles_ ,  The Journal of Chemical Physics 142(11), 114103 (2015),  10.1063/1.4914322,  https://doi.org/10.1063/1.4914322.
  * [25] T. Evensen, S. Naess and A. Elgsaeter,  _Brownian dynamics simulations of rotational diffusion using the cartesian components of the rotation vector as generalized coordinates_ ,  Macromolecular Theory and Simulations 17, 403 (2008),  10.1002/mats.200800031.
  * [26] B. Cichocki, M. L. Ekiel-Jeżewska and E. Wajnryb,  _Brownian motion of a particle with arbitrary shape_ ,  J. Chem. Phys. 142(21), 214902 (2015),  10.1063/1.4921729.
  * [27] I. M. Ilie, W. K. den Otter and W. J. Briels,  _Rotational brownian dynamics simulations of clathrin cage formation_ ,  J. Chem. Phys. 141(6), 065101 (2014),  10.1063/1.4891306.
  * [28] R. M. Jendrejack, J. J. de Pablo and M. D. Graham,  _Stochastic simulations of dna in flow: Dynamics and the effects of hydrodynamic interactions_ ,  J. Chem. Phys. 116(17), 7752 (2002),  10.1063/1.1466831,  https://doi.org/10.1063/1.1466831.
  * [29] R. G. Larson, T. T. Perkins, D. E. Smith and S. Chu,  _Hydrodynamics of a dna molecule in a flow field_ ,  Phys. Rev. E 55, 1794 (1997),  10.1103/PhysRevE.55.1794.
  * [30] A. Cressman, Y. Togashi, A. S. Mikhailov and R. Kapral,  _Mesoscale modeling of molecular machines: Cyclic dynamics and hydrodynamical fluctuations_ ,  Phys. Rev. E 77, 050901 (2008),  10.1103/PhysRevE.77.050901.
  * [31] C.-C. Hsieh, S. Jain and R. G. Larson,  _Brownian dynamics simulations with stiff finitely extensible nonlinear elastic-fraenkel springs as approximations to rods in bead-rod models_ ,  J. Chem. Phys. 124(4), 044911 (2006),  10.1063/1.2161210,  https://doi.org/10.1063/1.2161210.
  * [32] J. Huang and T. Schlick,  _Macroscopic modeling and simulations of supercoiled dna with bound proteins_ ,  J. Chem. Phys. 117(18), 8573 (2002),  10.1063/1.1511506,  https://doi.org/10.1063/1.1511506.
  * [33] B. A. Berg,  _Markov Chain Monte Carlo Simulations and Their Statistical Analysis_ ,  World Scientific,  10.1142/5602 (2004), https://www.worldscientific.com/doi/pdf/10.1142/5602.
  * [34] S. He and H. A. Scheraga,  _Macromolecular conformational dynamics in torsional angle space_ ,  J. Chem. Phys. 108(1), 271 (1998),  10.1063/1.475378,  https://doi.org/10.1063/1.475378.
  * [35] S. H. Northrup, J. A. Luton, J. O. Boles and J. C. L. Reynolds,  _Brownian dynamics simulation of protein association_ ,  Journal of Computer-Aided Molecular Design 1, 291 (1988),  10.1007/BF01677278.
  * [36] T. Yamanouchi and A. Satoh,  _Improvement of trapping performance of magnetic particles by magnetic multi-poles via Brownian dynamics simulations of magnetic rod-like particles in a Hagen-Poiseuille flow_ ,  Molecular Physics 120(11), e2067503 (2022),  10.1080/00268976.2022.2067503,  https://doi.org/10.1080/00268976.2022.2067503.
  * [37] R. Waszkiewicz and M. Lisicki,  _Hydrodynamic effects in the capture of rod-like molecules by a nanopore_ ,  J. Phys.: Cond. Matt. 33(10), 104005 (2020),  10.1088/1361-648x/abd11b.
  * [38] L. Qiao and G. W. Slater,  _Capture and translocation of a rod-like molecule by a nanopore: orientation, charge distribution and hydrodynamics_ ,  Phys. Chem. Chem. Phys. 24, 6444 (2022),  10.1039/D2CP00313A.
  * [39] S. A. Rogers, M. Lisicki, B. Cichocki, J. K. G. Dhont and P. R. Lang,  _Rotational diffusion of spherical colloids close to a wall_ ,  Phys. Rev. Lett. 109, 098305 (2012),  10.1103/PhysRevLett.109.098305.
  * [40] P. J. Zuk, B. Cichocki and P. Szymczak,  _GRPY: An Accurate Bead Method for Calculation of Hydrodynamic Properties of Rigid Biomacromolecules_ ,  Biophysical Journal 115(5), 782 (2018),  https://doi.org/10.1016/j.bpj.2018.07.015.
  * [41] J. Rotne and S. Prager,  _Variational Treatment of Hydrodynamic Interaction in Polymers_ ,  J. Chem. Phys. 50(11), 4831 (1969),  10.1063/1.1670977.
  * [42] H. Yamakawa,  _Transport Properties of Polymer Chains in Dilute Solution: Hydrodynamic Interaction_ ,  J. Chem. Phys. 53(1), 436 (1970),  10.1063/1.1673799.
  * [43] P. J. Zuk, E. Wajnryb, K. A. Mizerski and P. Szymczak,  _Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead models_ ,  J. Fluid Mech. 741, R5 (2014),  10.1017/jfm.2013.668.
  * [44] G. Nägele,  _Brownian Dynamics simulations_ ,  In S. Bluegel, G. Gompper, E. Koch, H. Mueller-Krumbhaar, R. Spatschek and R. G. Winkler, eds., _Computational Condensed Matter Physics_ , vol. B4:1 (2006).
  * [45] B. Cichocki, R. B. Jones, R. Kutteh and E. Wajnryb,  _Friction and mobility for colloidal spheres in Stokes flow near a boundary: The multipole method and applications_ ,  J. Chem. Phys. 112(5), 2548 (2000).
  * [46] M. L. Ekiel-Jeżewska and E. Wajnryb,  _Precise multipole method for calculating hydrodynamic interactions between spherical particles in the stokes flow_ ,  In F. Feuillebois and A. Sellier, eds., _Theoretical Methods in Micro Scale Viscous Flows_ , pp. 127–172. Transworld Research Network (2009).
  * [47] B. Cichocki, M. Rubin, A. Niedzwiecka and P. Szymczak,  _Diffusion coefficients of elastic macromolecules_ ,  J. Fluid Mech. 878, R3 (2019),  10.1017/jfm.2019.652.
  * [48] D. Frenkel and S. Berend,  _Understanding Molecular Simulation. From Algorithms to Applications._ ,  Elsevier (2002).


# Towards On-Device Face Recognition in Body-worn Cameras

Ali Almadan and Ajita Rattani  Dept. of Elec. Eng. and Computer Science
Wichita State University, Wichita, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Face recognition technology related to recognizing identities is widely
adopted in intelligence gathering, law enforcement, surveillance, and consumer
applications. Recently, this technology has been ported to smartphones and
body-worn cameras (BWC). Face recognition technology in body-worn cameras is
used for surveillance, situational awareness, and keeping the officer safe.
Only a handful of academic studies exist in face recognition using the body-
worn camera. A recent study has assembled BWCFace facial image dataset
acquired using a body-worn camera and evaluated the ResNet-50 model for face
identification. However, for real-time inference in resource constraint body-
worn cameras and privacy concerns involving facial images, on-device face
recognition is required. To this end, this study evaluates lightweight
MobileNet-V2, EfficientNet-B0, LightCNN-9 and LightCNN-29 models for face
identification using body-worn camera. Experiments are performed on a publicly
available BWCface dataset. The real-time inference is evaluated on three
mobile devices. The comparative analysis is done with heavy-weight VGG-16 and
ResNet-50 models along with six hand-crafted features to evaluate the trade-
off between the performance and model size. Experimental results suggest the
difference in maximum rank-1 accuracy of lightweight LightCNN-29 over best-
performing ResNet-50 is 1.85% and the reduction in model parameters is 23.49M.
Most of the deep models obtained similar performances at rank-5 and rank-10.
The inference time of LightCNNs is 2.1x faster than other models on mobile
devices. The least performance difference of 14% is noted between LightCNN-29
and Local Phase Quantization (LPQ) descriptor at rank-1. In most of the
experimental settings, lightweight LightCNN models offered the best trade-off
between accuracy and the model size in comparison to most of the models.

###### Index Terms:

Face Recognition, Body-worn Camera, Deep Learning, Edge AI, Mobile Computing,
Real-time Inference.

##  I Introduction

Facial recognition technology is used for verifying or identifying individuals
from a digital image or a video frame from a video source [1, 9]. Face
verification determines whether a pair of faces belong to the same identity,
while face identification classifies a _probe_ face image to a _template
gallery_ of the specific identity. It is a widely adopted technology for the
recognition of identities in surveillance, border control, healthcare, and
banking services.

A typical face recognition pipeline consists of face image acquisition, face
detection, facial image representation, and matching [1, 9]. Over the course
of time, various facial image representation methods have been proposed for
person recognition ranging from holistic such as Principal Component Analysis
(PCA) and Linear Discriminant Analysis (LDA) [4], to local image descriptors
such as Local Binary Patterns (LBP) [9] and Histogram of Oriented Gradients
(HOG) [23]. These local image descriptors also called _hand-crafted
descriptors or features_ have received significant attention over holistic
methods due to their lower computational cost and are non-training based
(i.e., hand-crafted features can be extracted from a single image without the
need of training a feature extraction module).

Deep learning has obtained significant advances in various computer vision
applications. Face recognition accuracy has been boosted to an unprecedented
level using deep learning [15]. Convolutional Neural Networks (CNNs) such as
VGG [17], ResNet [6], and InceptionNet [21] have been trained and fine-tuned
for face recognition and obtained hallmark accuracy rate of more than $99\%$
[15]. However, these high-performing networks have enormous space and
computational complexity due to the millions of parameters and computations
involved [8]. These requirements make the deployment of deep learning models
to resource-constrained environments such as mobile phones or wearable
devices, challenging. In order to meet the computational requirements of deep
learning in resource-constrained environments, a common approach is to
leverage cloud computing where deep learning training and inference is
performed on the server. However, cloud computing introduces  _high latency
and privacy concerns_.

Recently, face recognition technology has been ported to mobile and wearable
devices. For instance, Apple introduced the Face-ID moniker in its iPhone X
series for mobile user authentication. Further, recent interest has been drawn
on deploying face recognition technology on _body-worn cameras_ (BWC) by law
enforcement practitioners in order to keep officers safe, enable situational
awareness and surveillance [11, 13, 12]. With $88\%$ of Americans and $95\%$
of Dutch people supporting wearable cameras on police officers, there is
strong public support for this technology [12]. Only a handful of work exists
on face recognition111The term face recognition and identification is used
interchangeably. using body-worn cameras [14, 3, 10]. Recent study by the
authors [20] assembled _BWCFace dataset_ consisting of facial images captured
using body-worn cameras in in-door and daylight conditions, and evaluated
ResNet-50 [6] based CNN model for face identification on this data.

The _real-time inference_ is critical for applications involving face
recognition using body-worn cameras. Latency and privacy issues are involved
in the transfer of facial images to the cloud for network training and
inference. Therefore, towards on-device face recognition in body-worn cameras,
we evaluate lightweight CNN models in this study in order to evaluate the
tradeoff between the model size and the face identification accuracy obtained.

The contributions of this work are as follows:

  * •

Evaluation of light weight convolutional neural network (CNN) models namely,
MobileNet-V2 [28], EfficientNet-B0 [18], LightCNN-9 and LightCNN-29 [29] for
face identification via body-worn camera using publicly available BWCFace
dataset [20].

  * •

Comparative performance analysis is done with VGG-16 [17] and ResNet-50 [6]
models to understand the trade-off between the model size and the performance
obtained. In addition, six hand-crafted features namely, Local Binary Pattern
(LBP) [19], Modified Local Binary Patterns (mLBP) [24], Histograms of Oriented
Gradients (HOG) [23], Pyramid Histogram of Oriented Gradients (PHOG) [22],
Local Phase Quantization (LPQ) [25], and Local Ternary Patterns (LTP) [26] are
also evaluated for performance comparison.

This paper is organized as follows: Section II discusses the prior work on
face recognition using body-worn cameras. Deep learning-based models evaluated
in this study are discussed in section III. Hand-crafted descriptors used for
comparative performance evaluation are discussed in section IV. Experimental
protocol and the dataset is discussed in section V. Results are discussed in
section VI. Conclusions are drawn in section VII.

##  II Prior Work

One of the earliest studies on face recognition using body-worn cameras can be
dated back to $2011$ by Al-Obaydy and Sellahewa [14]. Authors collected videos
from $20$ subjects using an iOPTEC-P300 body-worn digital video camera and
have evaluated traditional subspace methods, namely, PCA, LDA, and discrete
wavelet transforms, for facial image representation. The cosine similarity
metric was used to compute scores between feature vectors from the pair of
facial images for the final identity assignment. Rank-1 recognition accuracy
values in the range [$65.83\%$, $76.04\%$] was obtained.

Brown and Fan [10] evaluated three face detectors: Viola & Jones, Aggregate
Channel Feature, and Faster R-CNN on an in-house dataset containing $638$ face
images captured using the body-worn camera. R-CNN approach achieved the
highest accuracy of $94\%$ over other face detection methods. The challenges
included detecting blurry, dark, and occluded faces captured in an
uncontrolled environment.

Recently, Almadan et al. [20] assembled a publicly available BWCFace dataset
consisting of a total of $178$K facial images from $132$ subjects captured
using body-worn cameras in in-door and daylight conditions. BWCFace is the
latest publicly available dataset with _largest number of subjects_ captured
using Body-worn camera till date. ResNet-50 [6] model was used for open-set
face identification on the collected dataset. Rank-1 to Rank-10 identification
accuracy values were used for performance evaluation. The model obtained the
highest Rank-1 accuracy of $99\%$ when templates and probes were acquired in
daylight. The model obtained the highest Rank-1 accuracy of $95.19\%$ in
cross-lighting conditions.

Table I summarizes the existing literature on facial analysis using the body-
worn camera.

TABLE I: Summary of the existing facial analysis studies using the body-worn camera. Reference | Dataset Size | Method
---|---|---
Al-Obaydy | 20 subjects | PCA, LDA
and Sellahewa [14] | x 96 images | and Discrete Wavelet
| (UBHSD dataset) |
Brown and Fan [10] |  $638$ face images | Viola & Jones,
|  | Agg. Feature,
|  | Faster R-CNN
Bryan [3] |  $3,600$ face images | Neurotechnology SDKs
Almadan et al. [20] |  $132$ subjects | ResNet-50
| x $178$K | for face identification
| (BWCFace dataset) |

##  III CNN Architectures

A Convolutional Neural Network (CNN) is a type of feed-forward artificial
neural network in which the connectivity pattern of its neurons, with shared
learnable weights and biases, is inspired by the organization of the visual
cortex (http://deeplearning.net/tutorial/lenet.html). A CNN consists of an
input and an output layer, as well as multiple hidden layers: convolutional
layers, pooling layers, fully connected layers, and normalization layers. We
evaluated lightweight; LightCNN, MobileNet, and EfficientNet for face
identification using the body-worn camera. ResNet-50 and VGG-16 are used for
comparative performance assessment with light CNN models. Deep features are
extracted from one of the fully connected layers of these trained CNNs for
facial image representation.

Next, we discuss these models as follows:

###  III-A Light-Weight CNN Models

  * •

LightCNN: This model extensively uses the Max-Feature-Map (MFM) operation
instead of ReLu activation, which acts as a feature filter after each
convolutional layer [29]. The operation takes two feature maps ($x^{1}$ and
$x^{2}$), eliminates the element-wise minimums, and returns element-wise
maximums $(h(x)$ as shown in Figure 1 (a). By doing so across feature
channels, only $50\%$ of the information-bearing nodes from each layer reach
the next layer. Consequently, during training, each layer is forced to
preserve only compact feature maps. Therefore, model parameters and the
extracted features are significantly reduced. We used the lightCNN-$29$ model
consisting of 12.6M parameters and 5.5M parameters for the LightCNN-$9$ model.

  * •

MobileNet: MobileNet [27] and [28] is one of the most popular mobile-centric
deep learning architectures which is small in size as well as computationally
efficient. The main idea of MobileNet is that instead of using regular
$3\times 3$ convolution filters, the operation is split into depth-wise
separable $3\times 3$ convolution filters followed by $1\times 1$ convolutions
(Figure 1 (b)). While achieving the same filtering and combination process as
a regular convolution, the new architecture requires a fewer number of
operations and parameters. In this study, we used MobileNet-V$2$ [28] which
consist of 3.4M parameters.

  * •

EfficientNet: The network consists of MBConv [18] which is similar to the
inverted residual blocks in MobileNetv2. This block adds squeeze-and-
excitation optimization. The residual blocks are capable of forming a shortcut
connection between the beginning of a convolutional block to the end. The
total number of parameters of the EfficientNet-B0 is 5.3M.

###  III-B Heavy-weight CNNs Models Used for Comparison

  * •

ResNet-50: ResNet [6] is a short form of residual network based on the idea of
“identity shortcut connection” where input features may skip certain layers.
The residual or shortcut connections introduced in ResNet allow for identity
mappings to propagate around multiple nonlinear layers, preconditioning the
optimization and alleviating the vanishing gradient problem. In this study, we
used ResNet-50 model which has 23.5M parameters.

  * •

VGG-16: VGG-16 [17] architecture consists of sequentially stacked $3\times 3$
convolutional layers with intermediate max-pooling layers followed by a fully
connected layers for feature extraction. The VGG-16 model has a total of 138M
parameters.

##  IV Handcrafted Descriptors

These hand-crafted features reported in the literature have been manually
designed, or “handcrafted” for feature representation from an image using
pixel intensity comparison or gradient information. The handcrafted features
used in this study are discussed below.

  * •

LBP: Local Binary Pattern (LBP) [19] is a low computational, efficient
descriptor that thresholds the neighbors of each pixel in an image. LBP
identifies local regions by considering a N x N size window to calculate the
subtraction between the central pixels and their neighboring pixels. If the
resulting intensity is greater or equal to one, then the local region is set
to 1; otherwise, it is set to zero. The binary code is computed for each
region and converted into a decimal number and the histogram of the
frequencies is computed and used as a feature vector.

  * •

mLBP: Modified Local Binary Patterns (mLBP) [24] compares the neighboring
pixels to the average of values in the N x N window, instead of the value of
the center pixel by LBP.

  * •

HOG: Histograms of Oriented Gradients (HOG) [23] can re-represent the shape of
a local object by using its distribution of edge directions without regards to
the spatial distribution of the image. Its advantage over other descriptors is
that it is variant only to object orientation.

  * •

PHOG: Pyramid Histogram of Oriented Gradients (PHOG) [22] is an extension to
the HOG descriptors. To preserve the spatial layout of the image, the PHOG
divides a given image into different cells at various pyramid levels. In
contrast, HOG operates at a fixed size image.

  * •

LPQ: Local Phase Quantization (LPQ) [25] is a LBP variant that is insensitive
to image blur. LPQ utilizes the blur invariance attribute of the Fourier phase
spectrum to extract the local phase of an image.

  * •

LTP: Local Ternary Patterns (LTP) [26] is another variant of LBP which is
powerful against noise; however, LTP may be sensitive to illumination change.
The LTP computes three-level thresholding, unlike LBP with two-level
thresholding.

These handcrafted features are used for facial feature extraction. We used
block size of $32$ and $3\times 3$ window size for handcrafted feature
extraction. These features are used for comparative performance analysis with
lightweight CNN models.

Figure 1: Main component of (a) LightCNN and (b) MobileNet. [29]

##  V Experimental Dataset and Protocol

In this section, we describe BWCFace dataset and the experimental protocol
used in this study.

BWCFace Dataset [20]: This dataset is collected from students at Wichita State
University using the Body-worn Camera (BWC) [20]. For the data collection,
videos from $132$ subjects were collected from Boblov $1296$P WiFi Body
Mounted Camera222https://www.boblov.com/. The subjects were recorded using a
$30$s high definition video of resolution $1920\times 1080$ pixels at $30$fps
in in-door and daylight settings ($10$ to $15$ minutes apart). The body-worn
camera was mounted on the center chest of the Research Assistant following the
recommendation in [3]. During recording, subjects were asked to act naturally
in order to represent realistic conditions. The recording was done in an
uncontrolled environment with non-uniform background and varying distances
between the subject and the body-worn camera.

Ethnic groups of the captured subjects include Asian, White American, Black,
and Middle Eastern. After the data collection, frames were extracted from the
collected videos, and facial regions were detected and cropped using Dlib
library [5].

Most of the images in BWCface are blurry or have a very low resolution because
of the shooting environment and the specification of the camera device, such
as not having autofocus. Also, the camera records at a wide angle of 170 ∘ for
clear and wide footage to cover the whole scene for evidence, unlike mugshots
or selfie images for the task of face recognition. In fact, the camera holder
and the target subjects are not stable while shooting in an un-collaborative
environment.

The selection of images was random from a large set of images at each specific
lighting condition. For instance, the chosen subsets of ”Day vs. Office” may
not be overlapping with the subset for ”Office vs. Day”. Therefore, there
might be a drop in performance for different scenarios.

The CNN models used in this study are trained from scratch on MS1Mv2 [7]
dataset. The MS1Mv2 dataset [7] is a cleaned version of the MS1M dataset [16],
containing around $5.8$ million images from $85,742$ subjects were used for
training the models. For the purpose of experiments, all the images were
resized to $224\times 224$. The models were trained using an early stopping
mechanism using Adam optimizer and cross-entropy loss function.

For the open-set evaluation of the trained deep learning models for face
identification on BWCFace dataset, the template gallery set consists of $12$
face images per subject. The probe set consists of $100$ face images randomly
selected per subject. From the gallery and probe set of all the subjects, deep
features of size $512$ are extracted for each image from the fully connected
layer of the CNN architectures discussed in Section III. The matching score
between the deep features from the pair of gallery and probe image is computed
using cosine similarity given in eq 1:

$\operatorname{sim}(u,v)=\frac{u\cdot
v}{|u||v|}=\frac{\sum_{i=1}^{N}u_{i}v_{i}}{\sqrt{\left(\sum_{i=1}^{N}u_{i}^{2}\right)\left(\sum_{i=1}^{N}v_{i}^{2}\right)}}$
(1)

where $u$ and $v$ are the two deep feature vectors:
$u=\\{u_{1},u_{2},\dots,u_{N}\\}$ and $v=\\{v_{1},v_{2},\dots,v_{N}\\}$.

The evaluation of handcrafted features follows the same protocol. The matching
scores of handcrafted features from a pair of gallery-probe face images were
calculated using Euclidean distance given in eq 2.

$d\left(u,v\right)=\sqrt{\sum_{i=1}^{N}\left(v_{i}-u_{i}\right)^{2}}$ (2)

where $u$ and $v$ are the handcrafted feature vectors.

The PyTorch framework was used for CNN training and feature extraction and
MATLAB was used for handcrafted feature extraction. The average of the scores
from the multiple gallery images per probe is used for the final identity
assignment. In cross-lighting conditions, the template (gallery set) was
selected from a different lighting condition with that of probe set for all
the identities.

TABLE II: Real-time Inference: Deep feature extraction time per sample measured in milliseconds (ms) on three real mobile devices for each CNN model. CNN Model | Num. of Params (M) | Device | Extraction Time (ms)
---|---|---|---
ResNet-50 |  23.5 | iPhone 6 |  2386
|  | iPhone X |  941
|  | iPhone XR |  834
VGG-16 |  138 | iPhone 6 | memory constraints
|  | iPhone X |  4138
|  | iPhone XR |  3747
MobileNetV2 |  3.4 | iPhone 6 |  1022
|  | iPhone X |  357
|  | iPhone XR |  240
EfficientNet-B0 |  5.3 | iPhone 6 |  1321
|  | iPhone X |  450
|  | iPhone XR |  329
LightCNN-29 |  12.6 | iPhone 6 |  1426
|  | iPhone X |  604
|  | iPhone XR |  430
LightCNN-9 |  5.5 | iPhone 6 |  562
|  | iPhone X |  248
|  | iPhone XR |  171
TABLE III: Face identification accuracy of deep CNN models on BWCFace dataset [20] in same and cross-lighting conditions at rank-1, rank-5 and rank-10. Light Condition | Rank-1 [%] | Rank-5 [%] | Rank-10 [%]
---|---|---|---
LightCNN-9
Office vs. Office | 93.35 | 97.78 | 98.56
Day. vs. Day | 78.70 | 90.47 | 93.76
Day vs. Office | 22.24 | 35.63 | 50.65
Office vs. Day | 22.41 | 32.36 | 44.75
LightCNN-29
Office vs. Office | 94.48 | 98.27 | 99.17
Day. vs. Day | 78.07 | 89.41 | 93.52
Day vs. Office | 22.82 | 32.86 | 46.29
Office vs. Day | 22.93 | 37.32 | 47.18
MobileNet-V2
Office vs. Office | 85.50 | 97.38 | 98.99
Day. vs. Day | 72.43 | 93.98 | 97.72
Day vs. Office | 25.18 | 33.30 | 44.60
Office vs. Day | 22.35 | 32.30 | 44.60
EfficientNet-B0
Office vs. Office | 80.95 | 94.93 | 97.61
Day. vs. Day | 68.24 | 91.18 | 96.13
Day vs. Office | 27.23 | 42.23 | 56.08
Office vs. Day | 27.92 | 41.92 | 55.59
VGG-16
Office vs. Office | 81.65 | 95.32 | 97.71
Day. vs. Day | 69.49 | 88.73 | 94.53
Day vs. Office | 31.43 | 49.92 | 68.08
Office vs. Day | 27.22 | 44.48 | 58.34
ResNet-50 [20]
Office vs. Office | 96.33 | 98.47 | 99.00
Day. vs. Day | 99.36 | 99.81 | 99.85
Day vs. Office | 84.25 | 95.61 | 98.06
Office vs. Day | 94.45 | 98.99 | 99.63
TABLE IV: Face identification accuracy of the handcrafted descriptors on BWCFace dataset [20] for same and cross-lighting conditions at rank-1, rank-5 and rank-10. Light Condition | Rank-1 [%] | Rank-5 [%] | Rank-10 [%]
---|---|---|---
LBP
Office vs. Office | 60.05 | 83.20 | 90.61
Day. vs. Day | 55.13 | 82.18 | 91.01
Office vs. Day | 16.66 | 44.07 | 62.16
Day vs. Office | 14.55 | 40.22 | 56.65
mLBP
Office vs. Office | 56.33 | 81.05 | 88.53
Day. vs. Day | 51.82 | 80.37 | 89.87
Office vs. Day | 16.79 | 43.87 | 62.57
Day vs. Office | 16.18 | 41.86 | 56.75
HOG
Office vs. Office | 65.90 | 84.26 | 89.51
Day. vs. Day | 70.39 | 89.46 | 94.12
Office vs. Day | 20.64 | 47.31 | 64.01
Day vs. Office | 17.96 | 45.67 | 63.04
PHOG
Office vs. Office | 55.10 | 76.71 | 84.31
Day. vs. Day | 64.09 | 84.44 | 91.20
Office vs. Day | 11.38 | 33.10 | 45.47
Day vs. Office | 10.49 | 32.86 | 49.73
LPQ
Office vs. Office | 80.06 | 92.29 | 95.48
Day. vs. Day | 77.77 | 91.78 | 95.68
Office vs. Day | 28.67 | 56.71 | 69.71
Day vs. Office | 26.59 | 54.82 | 69.92
LTP
Office vs. Office | 75.41 | 90.26 | 93.99
Day. vs. Day | 10.44 | 32.06 | 45.69
Office vs. Day | 10.85 | 33.31 | 47.49
Day vs. Office | 10.92 | 38.29 | 51.94

##  VI Experimental Results

Table III shows the face identification accuracy using deep features extracted
from the CNN architectures detailed in section III. The deep features have
been evaluated in the same and across lighting conditions. For instance,
office vs. day represents cross lighting conditions in which templates are
acquired in the office light and probe images are acquired in daylight
condition.

At office vs. office lighting conditions, both LightCNNs scored an average of
93.92% at rank-1, followed by $85.50\%$ for MobileNet-V2 at the same rank.
VGG16 and EfficientNet-B0 obtained rank-1 accuracy of 81.65% and 80.95%,
respectively. The highest rank-1 accuracy of $96.33\%$ was obtained by
ResNet-50 over all the models. However, all the models performed similarly at
rank-5 and rank-10 ranging from $97.38\%$ to $99.17\%$ accuracy.

For the day vs. day, the identification accuracy of the ResNet-50 model at
rank-1 was higher than VGG16 and EfficientNet-B0 by about $14.68\%$,
MobileNet-V2 by $10.83\%$ points, and than both LightCNNs by $2.42\%$ points.
However, ResNet-50 and LightCNN-29 had almost similar results at rank-5 and
rank-10 with a minimal difference in the accuracy of $0.2\%$.

The difference in rank-1 accuracy of best-performing lightweight LightCNN over
best performing ResNet-50 is $4.88\%$ and the drop in model parameters is
$23.49M$. The overall difference in the performance of LightCNN-9 over
LightCNN-29 is $0.39\%$ at rank-1 and rank-10, and $0.41\%$ at rank-5. The
difference in the model parameters is $7.1K$.

However, for all the deep models, performance degraded significantly across
lighting conditions. Across lighting, a maximum accuracy of $94.45\%$ at
rank-1 was obtained by ResNet-50 in the day vs. office conditions followed by
EfficientNet-B0 and VGG-16 obtaining $27.57\%$ at the same rank and lighting
condition.

ResNet-50 model obtained a least performance drop between same-lighting and
cross-lighting of $8.5\%$, $1.84\%$, $0.58\%$ at rank-1, rank-5, and rank-10,
respectively. The accuracy dropped by about $5\%$ for lighter models in cross-
lighting condition.

Table IV shows the face identification accuracy obtained by handcrafted
descriptors (LPB, mLPB, HOG, PHOG, LPQ, and LTP) for same-lighting and cross-
lighting conditions. It can be observed that the best rank-1 identification
accuracy of $80.06\%$ is obtained using the LPQ descriptor, while the accuracy
of $75.41\%$ was obtained using the LTP descriptors. The rest of the
descriptors scored similarly across all ranks with an average of $59.34\%$,
$81.31\%$, and $88.24\%$ at rank-1, rank-5, and rank-10, respectively. LTP and
LPQ descriptors scored the highest with a mean accuracy of $91.28\%$ at rank-5
and $94.73\%$ at rank-10. For hand-crafted features, the average decrease in
accuracy in cross-lighting conditions is approximately $43.40\%$, $37.99\%$,
$29.21\%$ at rank-1, rank-5, and rank-10, respectively.

In summary, the difference in maximum rank-1 accuracy of highest performing
lightweight LightCNN-29 over ResNet-50 and VGG-16 is $-1.85\%$ and $+12.83\%$.
The reduction in the model parameters of lightCNN-29 over ResNet-50 and VGG-16
are $23.49M$ and $137.98M$, respectively. The difference further reduces to
$1.54\%$ and $8.5\%$, and to $0.68\%$ and $1.46\%$, at rank-5 and rank-10,
respectively. Further, the difference in maximum rank-1 accuracy of
lightweight LightCNN-29 over best performing LPQ is $14.42\%$ and $5.74\%$ in
the same and cross-lighting condition. This difference reduces to $5.98\%$ and
$3.69\%$ at rank-5 and rank-10, respectively. The _superior performance of LPQ
could be due to its blur invariance property and the face images captured
using the body-worn camera contain motion blur_. The performance of most of
the models and feature descriptors dropped significantly in cross-lighting
conditions. ResNet-50 obtained the best overall performance. Among,
LightCNN-29 and LightCNN-9, the difference in the model parameters is $7.1K$
and the accuracy difference is $0.39\%$ at rank-1. The _superior performance
of LightCNN over other lightweight models could be attributed to the use of
Max-Feature-Map operation, instead of ReLu, which results in compact and
optimized deep feature representation_.

Table LABEL:tab:inf_time shows the real-time inference in terms of feature
extraction time per sample on three iPhone devices and the testing environment
is Xcode version 12.4. LightCNN-9 obtained the least feature extraction time
being 2.1x faster across the three iPhones than EfficientNets. EfficientNets
obtained the second lowest execution time although being 1.5x larger than
MobileNets. This is due to the efficient, systematic scaling up of the depth,
width, and regulation of CNN models.

These results suggest LightCNN models offer the best trade-off between model
size and the recognition performance for on-device face recognition in body-
worn cameras.

##  VII Conclusion

This paper evaluates different convolutional neural networks (CNNs) and six
handcrafted descriptors for face identification using a body-worn camera using
a publicly available BWCFace dataset. Deep features are extracted using light
and heavy CNN models to measure the trade-off between the model size and
identification performance obtained in order to facilitate on-device face
recognition. Real-time inference is evaluated on three mobile devices. The
experimental results suggest the efficacy of LightCNN models in offering the
best tradeoff between the model size (being 2.1x faster) and the recognition
performance for the body-worn camera. This could be attributed to the use of
Max-Feature-Map operation in LightCNN models which acts as a feature filter in
obtaining compact and optimized deep feature representation. Mostly, a large
performance difference is noted between lightweight CNNs and the handcrafted
descriptors, especially at rank-1. LPQ obtained superior performance over
other descriptors which could be due to its robustness to motion blur present
in facial images captured using the body-worn camera. As a part of future
work, custom compact size CNN models will be developed and evaluated against
LightCNN models and ResNet-50, for on-device face recognition using the body-
worn camera.

## Acknowledgment

This work is supported in part by Award for Research/Creative projects by
Wichita State University. BWCFace dataset is available by contacting the
authors.

## References

  * [1] Anil K Jain and Stan Z Li.  Handbook of face recognition, volume 1.  Springer, 2011.
  * [2] Timo Ahonen, Abdenour Hadid, and Matti Pietikainen.  Face Description with Local Binary Patterns: Application to Face Recognition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(12):2037–2041, December 2006.
  * [3] Julia Bryan.  Effects of Movement on Biometric Facial Recognition in Body-Worn Cameras.  PhD thesis, Purdue University Graduate School, 2020.
  * [4] Kamran Etemad and Rama Chellappa.  Discriminant analysis for recognition of human face images.  Josa a, 14(8):1724–1733, 1997.
  * [5] Davis E King.  Dlib-ml: A machine learning toolkit.  Journal of Machine Learning Research, 10(Jul):1755–1758, 2009.
  * [6] Kaiming He, Xiangyu Zhang, Shaoqing Ren, Jian Sun.  Deep Residual Learning for Image Recognition.  In arXiv, 2015.
  * [7] Jiankang Deng, Jia Guo, Niannan Xue, and Stefanos Zafeiriou.  ArcFace: Additive Angular Margin Loss for Deep Face Recognition.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 4685–4694, Long Beach, CA, USA, June 2019. IEEE.
  * [8] J. Chen and X. Ran.  Deep Learning With Edge Computing: A Review.  In Proceedings of the IEEE, pages 1655-1674, 107, 8, 2019.
  * [9] Timo Ahonen, Abdenour Hadid, and Matti Pietikainen.  Face Description with Local Binary Patterns: Application to Face Recognition.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(12):2037–2041, December 2006.
  * [10] Lisa M Brown and Quanfu Fan.  Enhanced face detection using body part detections for wearable cameras.  In 2016 23rd International Conference on Pattern Recognition (ICPR), pages 715–720. IEEE, 2016.
  * [11] Jay Stanley.  Police body-mounted cameras: With right policies in place, a win for all.  American Civil Liberties Union, 2015.
  * [12] Daniel E Bromberg, Étienne Charbonneau, and Andrew Smith.  Public support for facial recognition via police body-worn cameras: Findings from a list experiment.  Government Information Quarterly, 37(1):101415, 2020.
  * [13] Kelly Blount.  Body worn cameras with facial recognition technology: when it constitutes a search.  Crim. L. Prac., 3:61, 2015.
  * [14] Wasseem Al-Obaydy and Harin Sellahewa.  In Vielhauer C., Dittmann J., Drygajlo A., Juul N.C., Fairhurst M.C. (eds) Biometrics and ID Management. BioID 2011.  Lecture Notes in Computer Science, vol 6583. Springer, Berlin, Heidelberg.
  * [15] Shervin Minaee, Amirali Abdolrashidi, Hang Su, Mohammed Bennamoun, and David Zhang.  Biometric recognition using deep learning: A survey.  arXiv preprint arXiv:1912.00271, 2019.
  * [16] Yandong Guo, Lei Zhang, Yuxiao Hu, Xiaodong He, and Jianfeng Gao.  Ms-celeb-1m: A dataset and benchmark for large-scale face recognition.  In European conference on computer vision, pages 87–102. Springer, 2016.
  * [17] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [18] Mingxing Tan and Quoc Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In International Conference on Machine Learning, pages 6105–6114. PMLR, 2019.
  * [19] Timo Ojala, Matti Pietikäinen, and David Harwood.  A comparative study of texture measures with classification based on featured distributions.  Pattern recognition, 29(1):51–59, 1996.
  * [20] Ali Almadan, Anoop Krishnan, and Ajita Rattani.  Bwcface: Open-set face recognition using body-worn camera.  19th IEEE International Conference on Machine Learning and Applications (ICMLA), 2020.
  * [21] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1–9, 2015.
  * [22] Anna Bosch, Andrew Zisserman, and Xavier Munoz.  Representing shape with a spatial pyramid kernel.  In Proceedings of the 6th ACM international conference on Image and video retrieval, pages 401–408, 2007.
  * [23] Navneet Dalal and Bill Triggs.  Histograms of oriented gradients for human detection.  In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 1, pages 886–893. Ieee, 2005.
  * [24] Timo Ojala, Matti Pietikainen, and Topi Maenpaa.  Multiresolution gray-scale and rotation invariant texture classification with local binary patterns.  IEEE Transactions on pattern analysis and machine intelligence, 24(7):971–987, 2002.
  * [25] Ville Ojansivu and Janne Heikkilä.  Blur insensitive texture classification using local phase quantization.  In International conference on image and signal processing, pages 236–243. Springer, 2008.
  * [26] Xiaoyang Tan and Bill Triggs.  Enhanced local texture feature sets for face recognition under difficult lighting conditions.  IEEE transactions on image processing, 19(6):1635–1650, 2010.
  * [27] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam.  Mobilenets: Efficient convolutional neural networks for mobile vision applications.  arXiv preprint arXiv:1704.04861, 2017.
  * [28] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov, and L.C. Chen.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In IEEE Conference on CVPR, pages 4510–4520, 2018.
  * [29] X. Wu, R. He, Z. Sun, and T. Tan.  A light cnn for deep face representation with noisy labels.  IEEE TIFS, 13(11):2884–2896, 2018.