# Coarse-grained curvature tensor on polygonal surfaces

Charlie Duclut<EMAIL_ADDRESS>Max-Planck-Institut für Physik Komplexer
Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany.  Aboutaleb Amiri
<EMAIL_ADDRESS>Max-Planck-Institut für Physik Komplexer Systeme, Nöthnitzer
Str. 38, 01187 Dresden, Germany.  Joris Paijmans  Max-Planck-Institut für
Physik Komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany.  Frank
Jülicher<EMAIL_ADDRESS>Max-Planck-Institut für Physik Komplexer
Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany.  Center for Systems
Biology Dresden, Pfotenhauerstr. 108, 01307 Dresden, Germany.  Cluster of
Excellence Physics of Life, TU Dresden, 01062 Dresden, Germany.

###### Abstract

Using concepts from integral geometry, we propose a definition for a local
coarse-grained curvature tensor that is well-defined on polygonal surfaces.
This coarse-grained curvature tensor shows fast convergence to the curvature
tensor of smooth surfaces, capturing with accuracy not only the principal
curvatures but also the principal directions of curvature. Thanks to the
additivity of the integrated curvature tensor, coarse-graining procedures can
be implemented to compute it over arbitrary patches of polygons. When computed
for a closed surface, the integrated curvature tensor is identical to a rank-2
Minkowski tensor. We also provide an algorithm to extend an existing C++
package, that can be used to compute efficiently local curvature tensors on
triangulated surfaces.

In the last decades, biophysics and soft matter physics have provided many
examples where the geometry of the system plays a major role to understand the
structural and mechanical properties of materials such as colloids, liquid
crystals, membranes, cells and tissues Seifert (1997); Bowick and Giomi
(2009); Shankar  _et al._ (2020); Paulose  _et al._ (2012); Mietke  _et al._
(2019); Messal  _et al._ (2019). In the case of a fluid membrane for instance,
the local mean and Gaussian curvatures of the surface are crucial ingredients
of its energy functional Helfrich (1973); Seifert (1997). The description of
membranes that exhibit a polar (or nematic) order requires not only the
knowledge of the principal curvatures, but that of the full local curvature
tensor, as couplings between the polar order parameter and the curvature
tensor are allowed by symmetry Helfrich and Prost (1988); MacKintosh and
Lubensky (1991). Moreover, active membranes such as the cell cortex or cell
monolayers can exert active tensions and active moments which generically
couple to the local curvature tensor Salbreux and Jülicher (2017). Such
couplings can for instance be a consequence of a preferential alignment of
filaments of the cytoskeleton along the axis of largest curvature. At the
macroscopic level, such systems can be described using tools from
(nonequilibrium) statistical mechanics Salbreux and Jülicher (2017) and from
differential geometry Deserno (2004), and the existence of a smooth manifold
and a well-defined curvature tensor is assumed. It is therefore crucial to
understand how the macroscopic continuum description in terms of local metric
and local curvature tensors can be obtained from the coarse-graining of
discrete objects such as molecules, colloids or cells.

Surfaces represented by discrete triangles or polygons are ubiquitous as they
arise from the discretization of a smooth surface for numerical purposes, or
from experimental data on real geometries. Defining curvature measures on such
polygonal or triangulated surfaces is a long-standing problem which has been
addressed extensively in mathematics Hug (2000); Schneider and Weil (2008);
Jensen and Kiderlen (2017). Despite these important contributions, that we
discuss more below, a robust and practical definition of a local curvature
tensor that converges to the continuum curvature tensor in the limit of
infinitesimal triangulation remains however elusive for such surfaces. Several
methods have been proposed to compute curvature tensors from triangulated
surfaces. Some of them rely on a local fit of the smooth surface to the
triangulated manifold Qingde Li and Griffiths (2004); Beale  _et al._ (2016).
The resulting curvature tensor depends on the method used and such approaches
are unsatisfactory in cases where the surface is intrinsically made of
triangles or polygons, as cellular tissues or foams, for instance. Other
methods are directly built on the triangulated surface, but their definitions
rely on arbitrary choices, leading to a proliferation of algorithms Taubin
(1995); Hameiri and Shimshoni (2003); Ramakrishnan  _et al._ (2010); Crane
(2020).

The development of integral geometry Hug (2000); Schneider and Weil (2008);
Jensen and Kiderlen (2017) has however provided physicists with new tools to
quantify shapes and curvature of surfaces, for instance by using Minkowski
functionals Mecke (2000); Schröder-Turk  _et al._ (2013). Minkowski
functionals are defined locally on open or closed surfaces and can take scalar
or tensorial values Hug (2000); Schneider and Weil (2008); Jensen and Kiderlen
(2017). They possess several properties that make them especially appealing:
they are additive, continuous, motion covariant (invariant for the scalars)
and span the space of scalar and tensor-valued valuations on convex shapes
Hadwiger (1951); Alesker (1999a, b); Schneider (2013); Hug and Schneider
(2014). The last property implies in particular that any extensive and motion
invariant scalar functional can be written as a linear combination of the
Minkowski scalars Mecke (2000). For a closed body embedded in a three-
dimensional space, there are only $d+1=4$ Minkowski scalars: the volume of
this body, its surface area, its integrated mean curvature and its integrated
Gaussian curvature (or Euler characteristic). Minkowski tensors are a
generalization of the Minkowski scalars to tensorial quantities Schröder-Turk
_et al._ (2011, 2013) and share their appealing properties. For the reasons
stated above, the Minkowski scalars have proven extremely robust to evaluate
the local curvature of surfaces Bian  _et al._ (2020); Jülicher (1996), while
the rank-2 Minkowski tensors have been used to analyze anisotropy in a wide
range of phenomena ranging from the shape of neuronal cells in the brain
Beisbart  _et al._ (2006) to the shape of galaxies in the universe Beisbart
_et al._ (2002).  Their robustness also permits to use them on cellular and
discretized structures, and a coarse-grained evaluation of Minkowski tensors
can even be obtained for pixelated images and polygonal surfaces defined on
grids Schröder-Turk  _et al._ (2013). Higher-rank Minkowski tensors have also
been used recently for shape reconstruction Kousholt and Kiderlen (2016);
Kousholt (2017), relying on the fact that a polytope in dimension $d$ with $m$
facets is uniquely determined by its surface tensors up to rank $m-d+2$
Kousholt (2017). Finally, we note that curvature-weighted tensorial measures
have also been used with success in density functional theories for fluids
made of arbitrarily-shaped convex hard particles Hansen-Goos and Mecke (2009,
2010).

In this paper, inspired by these concepts from integral geometry, we provide a
definition for a local integrated curvature tensor. This integrated curvature
tensor is defined as the sum of two local Minkowski tensors. It is well-
defined on polygonal surfaces and it is additive in the sense of Minkowski
functionals. Because of this property, we can use this tensor to define a
coarse-grained curvature tensor. We show that this coarse-grained curvature
tensor on arbitrary triangulated surfaces displays robust convergence to the
curvature tensor of the corresponding smooth surface.

The paper is organized as follows. We first define an integrated curvature
tensor for a smooth surface. Using a parallel body construction, we then
translate this definition to a triangulated surface, and discuss the
properties of this triangle-based curvature tensor. In Sec. II, we discuss the
convergence properties of this definition on a set of simple smooth surfaces.
We then present an algorithm to compute this curvature tensor on any
triangulated surface by extending the algorithm of Ref. Schröder-Turk  _et
al._ (2013). We illustrate this procedure on a few surfaces and compute
principal curvatures and directions of principal curvature of the coarse-
grained curvature tensor for these surfaces.

##  I Integrated curvature tensor

###  I.1 Differential geometry of curved surfaces

Figure 1: Notation used for a smooth surface (a) and a triangulated surface
(b). For the triangulated surface, the normal vector $\mathbf{n}$ is defined
everywhere on the surface (including on the sharp edges between triangles)
using a parallel body construction (see text).

We first consider the case of a smooth surface that can be described using
differential geometry. Let $\mathcal{S}$ be a two-dimensional surface
parameterized by two coordinates and let
$\mathbf{x}(s^{1},s^{2})\in\mathcal{S}$ be a point on this surface. At each
point on the surface we associate a local basis composed of two tangent
vectors $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$, and a normal vector
$\mathbf{n}$, defined as:

$\displaystyle\mathbf{e}_{1}=\frac{\partial\mathbf{x}}{\partial
s^{1}}\,,\quad\mathbf{e}_{2}=\frac{\partial\mathbf{x}}{\partial
s^{2}}\,,\quad\mathbf{n}=\frac{\mathbf{e}_{1}\times\mathbf{e}_{2}}{|\mathbf{e}_{1}\times\mathbf{e}_{2}|}\,,$
(1)

see Fig. 1 for illustration. The vectors $\mathbf{e}_{i}$ form the covariant
base of the tangent space. We furthermore define the contravariant base
$\mathbf{e}^{j}$ as $\mathbf{e}_{i}\cdot\mathbf{e}^{j}=\delta_{i}^{j}$, where
$\cdot$ denotes the scalar product, indices run from 1 to 2, and
$\delta_{i}^{j}$ is the Kronecker symbol. We also define the local metric
tensor with components $g_{ij}=\mathbf{e}_{i}\cdot\mathbf{e}_{j}$ in the local
basis $\mathbf{e}_{i}$. Indices can be raised or lowered by contraction with
the metric tensor: $a^{i}=g^{ij}a_{j}$ where a summation over repeated indices
is implied here and in the following. We denote by $\mathrm{d}\ell$ a line
element on the surface with
$\mathrm{d}\ell^{2}=g_{ij}\mathrm{d}s^{i}\mathrm{d}s^{j}$, and
$\mathrm{d}A=\sqrt{g}\mathrm{d}s^{1}\mathrm{d}s^{2}$ the components of an area
element, where $g={\rm det}\,g_{ij}$ is the determinant of the metric tensor.
We finally define the local curvature tensor with components $c_{ij}$ given
by:

$\displaystyle c_{ij}=-(\partial_{i}\partial_{j}\mathbf{x})\cdot\mathbf{n}\,,$
(2)

where $\partial_{i}=\partial/\partial s^{i}$. The coordinate-independent form
$\bm{c}$ of the curvature tensor is written as:

$\displaystyle\bm{c}=c_{ij}\mathbf{e}^{i}\odot\mathbf{e}^{j}\,,$ (3)

where $\bm{a}\odot\bm{b}\equiv(\bm{a}\otimes\bm{b}+\bm{b}\otimes\bm{a})/2$ is
the symmetric tensor product. The curvature tensor describes changes in the
local basis when moving on the surface as described by the Gauss–Weingarten
equations Deserno (2004):

$\displaystyle\nabla_{j}\mathbf{n}$ $\displaystyle=c_{ij}\mathbf{e}^{i}\,,$
(4a) $\displaystyle\nabla_{j}\mathbf{e}_{i}$
$\displaystyle=-c_{ij}\mathbf{n}\,,$ (4b)

where $\nabla_{j}$ denotes the covariant derivative Deserno (2004); Salbreux
and Jülicher (2017).

###  I.2 The integrated curvature tensor

Equation (3) is well-defined for smooth surfaces but difficult to define for
surfaces built from discrete components such as particles, molecules, polygons
or triangles. We therefore define a coarse-grained curvature measure defined
by integrating the curvature tensor over a finite surface patch $\mathcal{P}$.
This leads to the definition of the integrated curvature tensor $\bm{M}$:

$\displaystyle\bm{M}=\int_{\mathcal{P}}\mathrm{d}A\,\bm{c}=\int_{\mathcal{P}}\mathrm{d}A\,c_{ij}\mathbf{e}^{i}\odot\mathbf{e}^{j}\,.$
(5)

Using the Gauss–Weingarten equation (4a), the integrand of Eq. (5) can be
rewritten as

$\displaystyle
c_{ij}\mathbf{e}^{i}\odot\mathbf{e}^{j}=\nabla_{j}\mathbf{n}\odot\mathbf{e}^{j}$
$\displaystyle=\nabla_{j}(\mathbf{n}\odot\mathbf{e}^{j})-\mathbf{n}\odot\nabla_{j}\mathbf{e}^{j}$
$\displaystyle=\nabla_{j}(\mathbf{n}\odot\mathbf{e}^{j})+c_{i}^{i}\mathbf{n}\odot\mathbf{n}\,,$

where the second line has been obtained using Eq. (4b). Then, using the
divergence theorem on a curved surface Capovilla and Guven (2002); Salbreux
and Jülicher (2017), the integrated curvature tensor is rewritten as a sum of
a boundary contribution and a surface term:

$\displaystyle\bm{M}=\int_{\mathcal{P}}\mathrm{d}A\,c_{i}^{i}\mathbf{n}\odot\mathbf{n}+\oint_{\mathcal{C}}\mathrm{d}\ell\,\nu_{i}\mathbf{e}^{i}\odot\mathbf{n}\,.$
(6)

Here, $\mathcal{C}$ denotes a contour enclosing the surface patch
$\mathcal{P}$, and $\bm{\nu}=\nu_{i}\mathbf{e}^{i}$ is a unit vector, tangent
to $\mathcal{P}$, outward-pointing and normal to the contour $\mathcal{C}$
(see Fig. 1).

The definition (6) of the integrated curvature tensor has the form of a sum of
two Minkowski tensors. The surface term (first term in Eq. (6)) is
proportional to the local three-dimensional Minkowski tensor $W_{2}^{0,2}$
Schröder-Turk  _et al._ (2013). The boundary term (second term in Eq. (6)), is
formally a non-Euclidean two-dimensional Minkowski tensor embedded in three-
dimensional space. The representation of the integrated curvature tensor given
in Eq. (6) reduces for a closed surface $\mathcal{S}$ to a single rank-2
Minkowski tensor
$W_{2}^{0,2}=(1/6)\oint_{\mathcal{S}}\mathrm{d}A\,c_{i}^{i}\mathbf{n}\odot\mathbf{n}$
Schröder-Turk  _et al._ (2013), since the contour term vanishes in this case.
Importantly, both terms in Eq. (6) are additive and well-defined for a
triangulated surface, as we see in the next section.

Furthermore, the integrated mean curvature $M$ over the patch $\mathcal{P}$
can be directly obtained by taking the trace of the integrated curvature
tensor defined by Eq. (6):

$\displaystyle
M=\frac{1}{2}\int_{\mathcal{P}}\mathrm{d}A\,c_{i}^{i}=\frac{1}{2}{\rm
Tr}\,\bm{M}\,,$ (7)

since ${\rm Tr}\left(\mathbf{e}^{i}\odot\mathbf{n}\right)=0$ and ${\rm
Tr}\left(\mathbf{n}\odot\mathbf{n}\right)=1$.

Finally, the integrated curvature tensor can be normalized by the area
$A=\int_{\mathcal{P}}\mathrm{d}A$ of patch $\mathcal{P}$ to obtain the coarse-
grained curvature tensor $\bm{C}$:

$\displaystyle\bm{C}=\bm{M}/A\,,$ (8)

which will be crucial to define a local curvature tensor on polygonal
surfaces.

###  I.3 Integrated curvature tensor on triangulated surfaces

Figure 2: Definition and notation used for a triangulated surface and its
parallel body, following Ref. Schröder-Turk  _et al._ (2013). (a) Definition
of the edges using a doubly connected edge list data structure. Each edge
$\mathbf{e}$ is oriented and uniquely assigned to a triangle $T$, and has a
counter-oriented edge $\mathbf{e}^{\prime}$ assigned to the adjacent triangle
$T^{\prime}$. An oriented edge is also unambiguously assigned to the previous
edge $\mathbf{e}_{\rm previous}$ and to the next edge $\mathbf{e}_{\rm next}$.
(b) Cross-sectional view along the oriented edge $\mathbf{e}$. The parallel
surface of width $\varepsilon$ regularizes the triangulated surface and allows
the computation of integrated surface quantities. See main text for the
definitions of the quantities.

In the following, we consider triangulated surfaces for simplicity. Note that
our formalism can be applied to any polygonal surface. The parallel body
construction of a triangulated surface, that we describe below, will allow us
to define the integrated curvature tensor of surface patches represented by
planar triangles and even for a single triangle by using Eq. (6). For this
purpose, we first define some notation following Ref. Schröder-Turk  _et al._
(2013). On the triangulated surface, each edge $\mathbf{e}$ is oriented and
uniquely assigned to a triangle $T$, and has a counter-oriented edge
$\mathbf{e}^{\prime}$ assigned to the adjacent triangle $T^{\prime}$ (see Fig.
2). An oriented edge is also unambiguously assigned to the previous edge
$\mathbf{e}_{\rm previous}$ and the next edge $\mathbf{e}_{\rm next}$. This
description corresponds to a doubly connected edge list data structure. The
normal vector on a triangle $T$ is defined as $\mathbf{n}_{T}=(\mathbf{e}_{\rm
previous}\times\mathbf{e})/|\mathbf{e}_{\rm previous}\times\mathbf{e}|$.

The normal vectors $\mathbf{n}_{T}$ and $\mathbf{n}_{T^{\prime}}$ of triangles
$T$ and $T^{\prime}$ span the angle $\alpha_{\mathbf{e}}\in(-\pi,\pi]$ at
their common edge $\mathbf{e}$. Along each edge $\mathbf{e}$, we define the
normalized edge vector $\hat{\mathbf{e}}$, the mean normal vector
$\ddot{\mathbf{n}}_{\mathbf{e}}$ and its normal
$\dot{\mathbf{n}}_{\mathbf{e}}$, such that

$\displaystyle\hat{\mathbf{e}}=\frac{\mathbf{e}}{|\mathbf{e}|}\,,\quad\ddot{\mathbf{n}}_{\mathbf{e}}=\frac{\mathbf{n}_{T}+\mathbf{n}_{T^{\prime}}}{|\mathbf{n}_{T}+\mathbf{n}_{T^{\prime}}|}\,,\quad\dot{\mathbf{n}}_{\mathbf{e}}=\hat{\mathbf{e}}\times\ddot{\mathbf{n}}_{\mathbf{e}}\,$
(9)

form an orthonormal local basis.

By inflating the triangulated surface by an infinitesimal width $\varepsilon$,
we construct a parallel smooth surface that regularizes the sharp edges
between triangles by portions of cylinders and portions of spheres which allow
the computation of integrated surface quantities (see Fig. 2). We first focus
on the computation of the integrated curvature tensor of a single triangle, as
the integrated curvature of a patch of triangles is obtained by summing the
contribution of individual triangles (see below). The integrated curvature
tensor of a single triangle is defined by computing Eq. (6) over the triangle
and over the portions of the infinitesimal cylinders that connect it to its
neighbors and then by taking the limit $\varepsilon\to 0$. The area
contribution appearing in Eq. (6) for a single triangle reads:

$\displaystyle\int_{\mathcal{P}}\mathrm{d}A\,\mathbf{n}\odot\mathbf{n}c_{i}^{i}$
$\displaystyle=\lim_{\varepsilon\to 0}\sum_{\mathbf{e}\in
T}\int_{0}^{|\mathbf{e}|}\mathrm{d}\ell\int_{-\alpha_{\mathbf{e}}/2}^{\alpha_{\mathbf{e}}(a_{\mathbf{e}}-1/2)}\varepsilon\mathrm{d}\theta\mathbf{n}\odot\mathbf{n}\frac{1}{\varepsilon}\,,$
(10)

where the sum runs over the edges of triangle $T$ and where $\mathbf{n}$ is
the normal vector along the portion of cylinder and is parameterized by the
angle $\theta$ such that
$\mathbf{n}=\cos\theta\,\ddot{\mathbf{n}}_{\mathbf{e}}+\sin\theta\,\dot{\mathbf{n}}_{\mathbf{e}}$
(see App. A and Fig. 8 for details). Note that we integrate over a weighted
portion of the cylinder and have therefore introduced the weight fraction
$a_{\mathbf{e}}$, which we discuss in the following. Similarly, the boundary
contribution is given by:

$\displaystyle\oint_{\mathcal{C}}\mathrm{d}\ell\,\nu_{i}\mathbf{e}^{i}\odot\mathbf{n}=\lim_{\varepsilon\to
0}\sum_{\mathbf{e}\in
T}\int_{0}^{|\mathbf{e}|}\mathrm{d}\ell\,\bm{\nu}(\theta_{\rm
f})\odot\mathbf{n}(\theta_{\rm f})\,,$ (11)

where $\theta_{\rm f}=\alpha_{\mathbf{e}}(a_{\mathbf{e}}-1/2)$ is the final
angle of the portion of cylinder in the local basis defined by Eq. (9), and
$\bm{\nu}(\theta_{\rm f})=\mathbf{n}(\theta_{\rm f}+\pi/2)$
$=\hat{\mathbf{e}}\times\mathbf{n}(\theta_{\rm f})$ (see Fig. 8).

Evaluating and summing the contributions from Eqs. (10) and (11), we obtain
the integrated curvature tensor $\bm{M}_{T}$ of a single triangle as:

$\displaystyle\begin{split}\bm{M}_{T}&=\frac{1}{4}\sum_{\mathbf{e}\in
T}|\mathbf{e}|\bigl{[}(2a_{\mathbf{e}}\alpha_{\mathbf{e}}+\sin\alpha_{\mathbf{e}}+\sin(\alpha_{\mathbf{e}}-2a_{\mathbf{e}}\alpha_{\mathbf{e}}))\ddot{\mathbf{n}}_{\mathbf{e}}\odot\ddot{\mathbf{n}}_{\mathbf{e}}\\\
&+(2a_{\mathbf{e}}\alpha_{\mathbf{e}}-\sin\alpha_{\mathbf{e}}-\sin(\alpha_{\mathbf{e}}-2a_{\mathbf{e}}\alpha_{\mathbf{e}}))\dot{\mathbf{n}}_{\mathbf{e}}\odot\dot{\mathbf{n}}_{\mathbf{e}}\\\
&+4\cos(a_{\mathbf{e}}\alpha_{\mathbf{e}})\cos(\alpha_{\mathbf{e}}-a_{\mathbf{e}}\alpha_{\mathbf{e}})\ddot{\mathbf{n}}_{\mathbf{e}}\odot\dot{\mathbf{n}}_{\mathbf{e}}\bigr{]}\,.\end{split}$
(12)

Note that for concave triangles patches, the sign of $\alpha_{\mathbf{e}}$ in
Eq. (12) has to be changed (see App. A.3). The weight fraction
$a_{\mathbf{e}}$ determines the fraction of the integrated curvature
associated with bond $\mathbf{e}$ that is assigned to triangle $T$. Additivity
of the integrated curvature tensor imposes
$a_{\mathbf{e}}+a_{\mathbf{e}}^{\prime}=1$, where $a_{\mathbf{e}}^{\prime}$ is
the weight fraction associated to triangle $T^{\prime}$.

The integrated curvature tensor has the following properties. First, the
integrated curvature tensor $\bm{M}$ of a patch $P$ of triangles is directly
obtained by summing their individual contributions:

$\displaystyle\bm{M}=\sum_{T\in P}\bm{M}_{T}\,,$ (13)

where the sum runs over the triangles $T$ forming patch $P$. Furthermore,
similarly to Eq. (7) in the smooth case, the trace of the triangle-based
integrated curvature tensor reads:

$\displaystyle{\rm Tr}\,\bm{M}_{T}=\sum_{\mathbf{e}\in
T}|\mathbf{e}|a_{\mathbf{e}}\alpha_{\mathbf{e}}\,,$ (14)

and corresponds to the weighted sum of the integrated mean curvatures of each
edge. Note finally that $\bm{M}_{T}=0$ if the triangle $T$ is surrounded by
co-planar triangles (see App. A.2 for details).

###  I.4 Coarse-grained curvature tensor on triangulated surfaces

The integrated curvature tensor is a cumulative measure of curvature over a
surface patch. It is often convenient and useful to define coarse-grained
curvature measures, for example to define the curvature on a surface that is
formed by discrete particles. Starting from individual triangles, we define
the coarse-grained curvature tensor $\bm{C}_{T}$ of triangle $T$ as:

$\displaystyle\bm{C}_{T}=\bm{M}_{T}/A_{T}\,,$ (15)

where $A_{T}$ is the area of triangle $T$. Similarly, the coarse-grained
curvature tensor for a patch $P$ is given by Eq. (8) with $\bm{M}$ defined in
Eq. (13).

If a smooth surface is represented by a triangulation, we demand that the
coarse-grained curvature tensor of a triangle converges to the curvature of
the smooth surface in the limit where triangles become infinitesimal. This
requirement fixes the value of the weight $a_{\mathbf{e}}$ as follows. The
coarse-grained mean curvature $H_{T}=\sum_{\mathbf{e}\in T}H_{\mathbf{e}}$ of
triangle $T$ combines contributions
$H_{\mathbf{e}}=|\mathbf{e}|\alpha_{\mathbf{e}}a_{\mathbf{e}}/(2A_{T})$ from
all edges. Convergence to the curvature tensor of a smooth surface requires
that the mean curvature associated with a bond is the same for both adjacent
triangles $T$ and $T^{\prime}$. Thus,
$H_{\mathbf{e}}=H_{\mathbf{e}}^{\prime}$, which together with
$a_{\mathbf{e}}+a_{\mathbf{e}}^{\prime}=1$ uniquely determines

$\displaystyle a_{\mathbf{e}}=\frac{A_{T}}{A_{T}+A_{T^{\prime}}}\,.$ (16)

Equation (15), together with Eqs. (12) and (16), define the coarse-grained
curvature tensor, which is the main result of the paper. This coarse-grained
curvature tensor provides a robust definition of curvature on triangulated
surfaces, defined for a single triangle and its direct neighbors. The
definition (12) of the triangle-based integrated curvature tensor has several
appealing properties: (i) owing to the parallel body construction, it is well-
defined for triangulated surfaces, (ii) it is additive, and can therefore be
used for coarse-graining using Eq. (13), (iii) it displays robust convergence
when a triangulation approximates a smooth surface, as we show in the
following section.

##  II Approximation of smooth surfaces by triangulations

###  II.1 Convergence of the triangle-based integrated curvature tensor to a
continuum limit

(a) (b) (c)

Figure 3: Patch of triangles on a cylinder (a), on a sphere (b), and on an
ellipsoid (c). The triangle-based curvature tensor $\bm{C}_{T}$ is computed
for the central blue triangle $T$, using curvature information stemming form
its direct neighbors.

We have shown that the coarse-grained curvature tensor can be computed for a
triangulated surface. We now show that it converges to the curvature tensor of
a smooth surface in the limit where the triangle size becomes small. To this
end, we consider the curvature tensor of a smooth cylinder, sphere and
ellipsoid. For a smooth surface, the curvature tensor (3) is a $3\times 3$
matrix of rank 2, and has three eigenvalues. We denote by $\lambda_{1,2}$ the
two principal curvatures, and by $\lambda_{3}=0$ the last eigenvalue, whose
eigenvector is parallel to the surface normal. We compare the curvature tensor
of smooth surfaces with the coarse-grained curvature tensor $\bm{C}_{T}$ of a
single triangle $T$ surrounded by its three direct neighbors with triangle
vertices located on the surfaces (see Fig. 3). We determine the coarse-grained
curvature tensor of the central triangle using Eqs. (15) and (12) with an
area-weight (16). We denote by $\Lambda_{1,2}$ the principal curvatures of the
discrete curvature tensor that corresponds to $\lambda_{1,2}$, and by
$\Lambda_{3}$ the eigenvalue that converges to 0 and for which the associated
eigenvector is parallel to the surface normal as the areas of the triangles
become small.

Cylinder. We first consider an equilateral triangle of edge length $\ell$
surrounded by three equilateral triangles. All triangle vertices are located
on the surface of a cylinder of radius $R$. For the simple case where three of
the triangles are coplanars (see Fig. 3(a)), the eigenvalues read

$\displaystyle\begin{split}\Lambda_{1}&=\frac{\cos^{-1}\left(1-\frac{3\varepsilon^{2}}{8}\right)}{\sqrt{3}R\varepsilon}+\frac{\sqrt{2+3\varepsilon^{2}}}{2\sqrt{2}R}\,,\\\
\Lambda_{2}&=\frac{\cos^{-1}\left(1-\frac{3\varepsilon^{2}}{8}\right)}{\sqrt{3}R\varepsilon}-\frac{\sqrt{2+3\varepsilon^{2}}}{2\sqrt{2}R}\,,\\\
\Lambda_{3}&=0\,,\end{split}$ (17)

where $\varepsilon=\ell/R$. To lowest order in $\varepsilon$ we have:

$\displaystyle\Lambda_{1}=\frac{1}{R}+\frac{25\varepsilon^{2}}{64R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,\quad\Lambda_{2}=-\frac{23\varepsilon^{2}}{64R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,$
(18)

showing a convergence as $\varepsilon^{2}$ to the curvature tensor of a smooth
cylinder with $\lambda_{1}=1/R$ and $\lambda_{2}=0$. An alternative tiling,
where none of the triangles are co-planar, is provided in App. A.3. In this
case, we observe a slower convergence, linear in $\varepsilon$, to the same
eigenvalues.

(a) (b)

Figure 4:  Convergence of the triangle-based curvature tensor on an ellipsoid
with major axes $a=1,b=7,c=3/2$. (a) Polar angle dependence of the principal
curvatures of the smooth curvature tensor (solid lines) and of the discrete
one (black symbols). The triangle-based curvature tensor has been computed for
$\varphi=\pi/3$ and $\ell/a=10^{-3}$. (b) Convergence of the triangle-based
principal curvatures $\Lambda_{1,2}$ to the smooth principal curvatures
$\lambda_{1,2}$ for a patch of triangles located at $\theta=0$,
$\varphi=\pi/3$ (top) and $\theta=\pi/3$, $\varphi=\pi/3$ (bottom).

Sphere. Considering four equilateral triangles on a sphere of radius $R$, we
obtain the eigenvalues of the discrete curvature tensor $\bm{C}_{T}$, see Fig.
3(b). At lowest order in $\varepsilon=\ell/R$, we obtain:

$\displaystyle\Lambda_{1}=\frac{1}{R}+\frac{5\varepsilon^{2}}{18R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,\quad\Lambda_{2}=\frac{1}{R}+\frac{5\varepsilon^{2}}{18R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,\quad\Lambda_{3}=-\frac{5\varepsilon^{2}}{18R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,$
(19)

showing a convergence as $\varepsilon^{2}$ to the curvature tensor of a smooth
sphere with $\lambda_{1,2}=1/R$.

Ellipsoid. We now consider an ellipsoid defined by the equation

$\displaystyle x=a\sin\theta\cos\varphi\,,\quad
y=b\sin\theta\sin\varphi\,,\quad z=c\cos\theta\,.$ (20)

For simplicity, we generate a triangular patch on the ellipsoid by using four
equilateral triangles on the unit sphere and performing an affine
transformation of the sphere by rescaling the $x,y,z$ coordinates by the
factors $a,b,c$, respectively (see Fig. 3(c)). The curvature tensor on the
ellipsoid is a function of the polar and azimuthal angles $\theta$ and
$\varphi$. In Fig. 4(a), we show the comparison between the principal
curvatures of the discrete and of the smooth curvature tensor on an ellipsoid
with major axes $a=1,b=7,c=3/2$. We consider a patch of triangles with a
vertex of the central triangle located at $\theta=0,\varphi=\pi/3$. In this
case, the principal curvatures converge as $\varepsilon^{2}$ for small
$\varepsilon=\ell/a$ to the principal curvatures on the ellipsoid (see top
panel of Fig. 4(b)). The convergence is linear in $\varepsilon$ at
$\theta=\pi/3,\varphi=\pi/3$ (bottom panel of Fig. 4(b)). This difference in
convergence can be explained by the fact that for $\theta=0,\varphi=\pi/3$ all
triangles are equal in shape, while at $\theta=\pi/3,\varphi=\pi/3$ the
symmetry between triangles is broken. Note that the third eigenvalue
$\Lambda_{3}$ converges to 0 as $\varepsilon^{2}$ in both cases.

###  II.2 Directions of principal curvature on triangulated surfaces

Figure 5: Principal curvatures and principal directions of curvature of the
triangle-based curvature tensor for a Clifford torus (triangulated by 2880
triangles).

The definition (15) can be used to obtain the local principal directions of
the coarse-grained curvature tensor on triangulated surfaces. For this
purpose, we have extended the karambola package designed in Ref. Schröder-Turk
_et al._ (2013) to compute the triangle-based curvature tensor on arbitrary
surfaces. We provide the corresponding C++ code in App. B. As an illustration,
we have determined the principal curvatures and principal directions of
curvature of the triangle-based curvature tensor for three examples of closed
surfaces with different topology, see Figs. 5, 6 and 7.

Figure 5 displays the principal curvatures and principal directions of
curvature for a triangulated Clifford torus. The Clifford torus is an
axisymmetric torus with circular cross-section and a ratio $\sqrt{2}$ of the
radii of its two generating circles Willmore (1982). It has the important
property to minimize $G=\int_{\mathcal{S}}\mathrm{d}A\,(c_{i}^{i})^{2}$ for
toroidal topology. We indicate on each triangle the principal directions of
the triangle-based curvature tensors as bars with a color which indicates the
magnitude of the associated principal curvature. Note that one curvature is
constant and positive, while the other is positive (red) on the outer part of
the torus and negative (blue) on the inner part of the torus.

From this Clifford torus, we perform a special conformal transformation, which
consists of an inversion $\mathbf{R}^{\prime}=\mathbf{R}/R^{2}$, followed by a
translation by a vector $\bf t$ and a second inversion, such that a point
$\mathbf{R}$ is mapped to Jülicher (1996):

$\displaystyle\mathbf{R}^{\prime}=\frac{\mathbf{R}/R^{2}-\mathbf{t}}{(\mathbf{R}/R^{2}-\mathbf{t})^{2}}\,.$
(21)

For $\mathbf{t}=1.25r\,\mathbf{e}_{x}$, with $r$ the larger radius of the
Clifford torus and $\mathbf{e}_{x}$ a unit vector orthogonal to the axis of
symmetry of the torus, the transformation (21) applied to the Clifford torus
yields a non-axisymmetric torus with the same minimal value of $G$ (see Fig.
6). Indicated are again the principal curvatures and principal directions of
curvature of the triangle-based curvature tensor for each triangle.

The third example, displayed in Fig. 7, is the triangulated Lawson surface of
genus 2 Jülicher (1996); Lawson (1970). The Lawson surface minimizes $G$ for
genus-2 surfaces. The principal curvatures and principal directions of
curvature of the triangle-based curvature tensor are plotted on each triangle.
These examples show that the coarse-grained curvature tensor, defined on a
triangulation, provides a good approximation for the curvature tensor of the
underlying smooth surfaces.

Figure 6: Principal curvatures and principal directions of curvature of the
triangle-based curvature tensor for a Clifford torus after a special conformal
transformation defined in Eq. (21) (triangulated by 2880 triangles).  Figure
7: Principal curvatures and principal directions of curvature of the triangle-
based curvature tensor for the Lawson surface with genus 2 (triangulated by
24576 triangles).

## Conclusion

In this paper, starting from the definition of the integrated curvature tensor
for a smooth surface, we have used a parallel body construction to define an
integrated curvature tensor for polygonal and in particular for triangulated
surfaces. Parallel body construction is often used in the context of integral
geometry, and is at the heart of the definition and evaluation of Minkowski
functionals on polygonal manifolds Schröder-Turk  _et al._ (2013); Hug (2000);
Schneider and Weil (2008); Jensen and Kiderlen (2017). Our definition of the
integrated curvature tensor involves a sum of two Minkowski tensors (see Eq.
(6)). The first tensor is the surface contribution that is weighted by the
mean curvature of the surface patch on which it is evaluated and is commonly
used as a local curvature measure on its own. The second tensor is a boundary
contribution defined on the contour of the patch, and thus corresponds to a
two-dimensional Minkowski tensor embedded in three dimensional space. Such an
embedding of a Minkowski tensor to a higher dimensional space appears
naturally in our formalism and provides an extension of Minkowski tensors to
non-Euclidean geometries. Adding the surface and the boundary contributions
together provides the information to construct a local curvature tensor. Our
approach illustrates the importance of the boundary contribution to obtain a
robust convergence of the coarse-grained curvature tensor to that of smooth
surfaces. We note that to define the integrated curvature tensor, it is not
necessary to specify the value of the weight factor $a_{\mathbf{e}}$. However,
to define a coarse-grained curvature tensor, the weight factor takes a unique
value given by Eq. (16).

Importantly, our approach provides a definition of curvature on non-smooth
surfaces. This is relevant to build the continuum limit of the statistical
physics of curved manifolds starting from structures made of particles or
molecules. For example, in a fluid membrane molecular positions define
triangular meshes from which the coarse-grained curvature tensor introduced
here can be computed. Therefore, the differential geometry of a smooth surface
can emerge from a structure composed of discrete objects in a continuum limit.
The existence of such an explicit coarse-graining procedure is essential to
make sense and to provide explicit computational tools for theories where
couplings between the curvature tensor and an additional order parameter are
allowed by symmetry. This is the case for instance for equilibrium fluid
membranes that possess a nematic order parameter that can couple to the
surface geometry MacKintosh and Lubensky (1991); Lubensky and MacKintosh
(1993); Frank and Kardar (2008); Uchida (2002). It is also true for
nonequilibrium systems evolving on a curved manifold, where active couplings
between a polar or nematic order parameter and the curvature tensor are
generically allowed Salbreux and Jülicher (2017); Morris and Rao (2019);
Matoz-Fernandez  _et al._ (2020). Importantly, our coarse-grained curvature
tensor can also be used to compute the local curvature tensor of surfaces for
which keeping track of the discreteness of its constituents is relevant. It is
the case for instance of cell tissues. For such systems, triangle-based
methods have been developed to decompose the deformations of flat tissues into
cell contributions Merkel  _et al._ (2017) and have contributed to understand
the role of cellular processes in tissue patterning Etournay  _et al._ (2015).
For curved tissues, such method is still lacking and would require to also
quantify changes of the curvature tensor as the tissue deforms. Our coarse-
grained curvature tensor now provides the tool to perform such quantification.

Finally, a fast convergence of the coarse-grained curvature tensor to the
smooth curvature tensor has been demonstrated for various surfaces. Using the
karambola package designed in Ref. Schröder-Turk  _et al._ (2013) and the code
described in App. B, the coarse-grained curvature tensor of a triangulated
surface can be computed efficiently. Importantly, the principal directions of
curvature are computed with accuracy, as we have illustrated on Figs. 5, 6 and
7. An efficient tool is thus available to quantify curvature on triangulated
surfaces, that could be used in the context of numerical simulation or
experimental data processing, and also for computer vision and image
processing.

### Acknowledgments

We thank K. Ishihara for bringing Ref. Schröder-Turk  _et al._ (2013) to our
attention and for interesting discussions.

##  Appendix A Integrated curvature tensor

###  A.1 Computation of the integrated curvature tensor for a single triangle
and its neighbors

We recall here the definition of the integrated curvature tensor $\bm{M}$ for
a smooth surface as defined in Eq. (6):

$\displaystyle\bm{M}=\int_{\mathcal{P}}\mathrm{d}A\,c_{i}^{i}\mathbf{n}\odot\mathbf{n}+\oint_{\mathcal{C}}\mathrm{d}\ell\,\nu_{i}\mathbf{e}^{i}\odot\mathbf{n}\,.$
(22)

For a triangulated surface, one can compute the integrated curvature tensor
using a parallel body construction as discussed in the main text. One can then
first compute the surface integral term in Eq. (22). Since this term is
weighted by the mean curvature, the flat part of the triangle with zero
curvature does not contribute. The portions of sphere that regularize the
corners of the triangles yield a contribution of order $\varepsilon$ since
their area is of order $\varepsilon^{2}$ while their mean curvature goes as
$1/\varepsilon$. Therefore, contributions from triangle corners vanish, and
only the infinitesimal portions of cylinder regularizing the edges have non-
vanishing contributions in the limit $\varepsilon\to 0$, which read:

$\displaystyle\int_{\mathcal{P}}\mathrm{d}A\,\mathbf{n}\odot\mathbf{n}c_{i}^{i}$
$\displaystyle=\lim_{\varepsilon\to 0}\sum_{\mathbf{e}\in
T}\int_{0}^{|\mathbf{e}|}\mathrm{d}\ell\int_{-\alpha_{\mathbf{e}}/2}^{\alpha_{\mathbf{e}}(a_{\mathbf{e}}-1/2)}\varepsilon\mathrm{d}\theta\mathbf{n}\odot\mathbf{n}\frac{1}{\varepsilon}\,,$
(23)

where the sum runs over the edges of triangle $T$. Note that we integrate over
a weighted portion of the cylinder and have therefore introduced the weighting
factor $a_{\mathbf{e}}$.

The vector $\mathbf{n}$ is the normal vector along the portion of cylinder and
is then parameterized by the angle $\theta$ such that
$\mathbf{n}=\cos\theta\,\ddot{\mathbf{n}}_{\mathbf{e}}+\sin\theta\,\dot{\mathbf{n}}_{\mathbf{e}}$
(see Fig. 8). The surface integral thus reads:

$\displaystyle\begin{split}\int_{\mathcal{P}}\mathrm{d}A\,c_{i}^{i}\mathbf{n}\odot\mathbf{n}&=\frac{1}{4}\sum_{\mathbf{e}\in
T}|\mathbf{e}|\bigl{[}(2a_{\mathbf{e}}\alpha_{\mathbf{e}}+\sin\alpha_{\mathbf{e}}-\sin(\alpha_{\mathbf{e}}-2a_{\mathbf{e}}\alpha_{\mathbf{e}}))\ddot{\mathbf{n}}_{\mathbf{e}}\odot\ddot{\mathbf{n}}_{\mathbf{e}}\\\
&+(2a_{\mathbf{e}}\alpha_{\mathbf{e}}-\sin\alpha_{\mathbf{e}}+\sin(\alpha_{\mathbf{e}}-2a_{\mathbf{e}}\alpha_{\mathbf{e}}))\dot{\mathbf{n}}_{\mathbf{e}}\odot\dot{\mathbf{n}}_{\mathbf{e}}\\\
&-4\sin(a_{\mathbf{e}}\alpha_{\mathbf{e}})\sin(\alpha_{\mathbf{e}}-a_{\mathbf{e}}\alpha_{\mathbf{e}})\ddot{\mathbf{n}}_{\mathbf{e}}\odot\dot{\mathbf{n}}_{\mathbf{e}}\bigr{]}\,,\end{split}$
(24)

Second, the contour integral in Eq. (6) is evaluated by noting that the unit
vector $\bm{\nu}$, tangent to the regularized triangle and normal to its
contour lies at the end of the portions of cylinder, such that:

$\displaystyle\oint_{\mathcal{C}}\mathrm{d}\ell\,\nu_{i}\mathbf{e}^{i}\odot\mathbf{n}=\lim_{\varepsilon\to
0}\sum_{\mathbf{e}\in
T}\int_{0}^{|\mathbf{e}|}\mathrm{d}\ell\,\bm{\nu}(\theta_{\rm
f})\odot\mathbf{n}(\theta_{\rm f})\,,$ (25)

where $\theta_{\rm f}=\alpha_{\mathbf{e}}(a_{\mathbf{e}}-1/2)$ is the final
angle of the portion of cylinder in the local basis defined by Eq. (9), and
$\bm{\nu}(\theta_{\rm f})=\mathbf{n}(\theta_{\rm f}+\pi/2)$
$=\hat{\mathbf{e}}\times\mathbf{n}(\theta_{\rm f})$ (see Fig. 8). Note that
the portions of sphere that regularize the corners of the triangle do not
contribute to the previous equation as the contour length on each portion of
sphere is of order $\varepsilon$ and vanishes in the limit $\varepsilon\to 0$.
We thus obtain for the contour term:

$\displaystyle\oint_{\mathcal{C}}\mathrm{d}\ell\,\nu_{i}\mathbf{e}^{i}\odot\mathbf{n}$
$\displaystyle=\frac{1}{2}\sum_{\mathbf{e}\in
T}|\mathbf{e}|\bigl{[}\sin(\alpha_{\mathbf{e}}-2a_{\mathbf{e}}\alpha_{\mathbf{e}})\left(\ddot{\mathbf{n}}_{\mathbf{e}}\odot\ddot{\mathbf{n}}_{\mathbf{e}}-\dot{\mathbf{n}}_{\mathbf{e}}\odot\dot{\mathbf{n}}_{\mathbf{e}}\right)+2\cos(\alpha_{\mathbf{e}}-2a_{\mathbf{e}}\alpha_{\mathbf{e}})\ddot{\mathbf{n}}_{\mathbf{e}}\odot\dot{\mathbf{n}}_{\mathbf{e}}\bigr{]}\,.$
(26)

Summing Eqs. (24) and (26) yields Eq. (12) in the main text. For a patch $P$
of triangles, it is clear from Eqs. (23) and (25) that the integrated
curvature tensor $\bm{M}$ of the patch is given by adding individual triangles
contribution, and yields Eq. (13) in the main text.

Figure 8: Cross-sectional view of neighboring triangles in a convex (a), and
concave (b) configurations. The angles $\theta_{i}$ and $\theta_{f}$ denote
the lower and upper limits of integration, respectively.

###  A.2 Integrated curvature tensor for co-planar triangles

As expected, the integrated curvature tensor given by Eq. (12) vanishes for a
triangle $T$ surrounded by co-planar triangles. This fact is however not
obvious from Eq. (12) and we give details here. The first two terms in the sum
appearing in Eq. (12) vanish since $\alpha_{\mathbf{e}}=0$ for all edges for
co-planar triangles. In addition,
$\ddot{\mathbf{n}}_{\mathbf{e}}=\mathbf{n}_{T}$ is the triangle normal and is
the same for all edges, and
$\dot{\mathbf{n}}_{\mathbf{e}}=\bm{\nu}_{\mathbf{e}}$ is the unit vector
normal to each edge in the plane of the triangle. The last term in Eq. (12)
thus reduces to $\mathbf{n}_{T}\odot\sum_{\mathbf{e}\in
T}|\mathbf{e}|\bm{\nu}_{\mathbf{e}}$. One can show, for instance by taking the
scalar product with $\bm{\nu}_{\mathbf{e}^{\prime}}$ and using the law of
cosines, that $\sum_{\mathbf{e}\in T}|\mathbf{e}|\bm{\nu}_{\mathbf{e}}$
vanishes for any triangle, This shows that the integrated curvature tensor
vanishes for a triangle surrounded by co-planar triangles.

###  A.3 Concave decomposition

The derivation of Eq. (12) assumed a local convexity of the triangulated
surface. As illustrated in Fig. 8, the generalization to concave edges is
however straightforward as it simply requires to change the sign of
$\alpha_{\mathbf{e}}$ for concave edges.

A simple illustration of this procedure occurs when considering a cylinder
with triangles disposed in the direction orthogonal to the cylinder axis, see
Fig. 9. In this case, the central triangle and the purple triangle form a
concave edge for which the procedure highlighted above must be applied.

Similarly to the parallel tiling of the cylinder presented in the main text,
we can obtain the eigenvalues of the coarse-grained curvature tensor and
expand them in series of $\varepsilon=\ell/R$ to obtain:

$\displaystyle\Lambda_{1}=\frac{1}{R}+\frac{53\varepsilon^{2}}{288R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,\quad\Lambda_{2}=\frac{\varepsilon}{8\sqrt{3}R}-\frac{7\varepsilon^{2}}{96R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,,\quad\Lambda_{3}=-\frac{\varepsilon}{8\sqrt{3}R}-\frac{7\varepsilon^{2}}{96R}+\mathcal{O}\\!\left(\varepsilon^{3}\right)\,.$
(27)

Note that the convergence of the second and third eigenvalues is only linear
in $\varepsilon$, while it is quadratic in the case of a parallel tiling (see
Eq.(17)).

Figure 9: Patch of triangles on a cylinder. Tiling is orthogonal to the axis
of the cylinder.

##  Appendix B C++ code for computing the triangle-based integrated curvature
tensor

We have extended the karambola package
(http://www.theorie1.physik.fau.de/karambola) to include the computation of
the coarse-grained curvature tensor, its eigenvalues and eigenvectors. We
provide below the C++ function that computes the integrated curvature tensor
on a triangulated surface using Eq. (12). Similarly to the computation of
other quantities that can be determined using the karambola package, the
integrated curvature tensor can also be computed for labeled patches on the
surface. If all triangles are labeled uniquely, the integrated curvature
tensor is computed for each triangle. Note that in order to obtain Figs. 5, 6
and 7, we have computed the integrated curvature tensor of each triangle and
normalized it by the triangle area to obtain the coarse-grained curvature
tensor of each triangle.

⬇

CompWiseMatrixMinkValResultType calculate_curvature(const Triangulation&
surface){

//calculate c_int (integrated curvature)

CompWiseMatrixMinkValResultType c_int;

for (unsigned int l = 0; l<surface.n_triangles();l++){

double A_t = surface.area_of_triangle(l);

for (unsigned int k = 0; k<3;k++){

if(surface.ith_neighbour_of_triangle(l,k) != NEIGHBOUR_UNASSIGNED){

double A_tp= surface.area_of_triangle(surface.ith_neighbour_of_triangle(l,k));

double ar= A_t/(A_t+A_tp);

Vector
n2=surface.normal_vector_of_triangle(surface.ith_neighbour_of_triangle(l,k));

Vector n1= surface.normal_vector_of_triangle(l);

double c= dot(n1,n2)/(norm(n1)*norm(n2));

assert (c <= 1.000001);

if (c>= 1) c=1.;

double alpha= acos(c);

\parVector com_l = surface.com_of_triangle(l); //center of mass of triangle

Vector com_k =
surface.com_of_triangle(surface.ith_neighbour_of_triangle(l,k));

Vector convex=com_l+n1-(com_k+n2);

Vector concave=com_l-n1-(com_k-n2);

if(norm(convex) < norm(concave)) alpha=-alpha;

\pardouble e_norm= surface.get_edge_length(l,k);

Vector e_c1= surface.get_pos_of_vertex(surface.ith_vertex_of_triangle(l,k));

Vector e_c2=
surface.get_pos_of_vertex(surface.ith_vertex_of_triangle(l,(k+1)\parVector e=
(e_c2 - e_c1)/e_norm;

Vector n_a= (n1 + n2)/(norm(n1+n2));

Vector n_i= cross_product(e,n_a);

\pardouble local_value= 0;

for (unsigned int i = 0; i<3;i++){

for (unsigned int j = 0; j<=i;j++){

double term1=(2*ar*alpha+sin(alpha)+sin(alpha-2*ar*alpha)) * n_a[i]*n_a[j];

double term2=(2*ar*alpha-sin(alpha)-sin(alpha-2*ar*alpha)) * n_i[i]*n_i[j];

double term3=2*cos(ar*alpha)*cos(alpha-
ar*alpha)*(n_a[i]*n_i[j]+n_a[j]*n_i[i]);

local_value= e_norm*( term1 + term2 + term3 );

c_int[surface.label_of_triangle(l)].result_(i,j) += 1/4.*local_value;}

}

}

}

}

return c_int;

}

## References

  * Seifert (1997) U. Seifert, “Configurations of fluid membranes and vesicles,” Adv. Phys. 46, 13–137 (1997).
  * Bowick and Giomi (2009) M. J. Bowick and L. Giomi, “Two-dimensional matter: Order, curvature and defects,” Adv. Phys. 58, 449–563 (2009).
  * Shankar  _et al._ (2020) S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, “Topological active matter,” arXiv:2010.00364 (2020), arXiv:2010.00364 .
  * Paulose  _et al._ (2012) J. Paulose, G. A. Vliegenthart, G. Gompper, and D. R. Nelson, “Fluctuating shells under pressure,” Proc. Natl. Acad. Sci. U.S.A. 109, 19551–19556 (2012).
  * Mietke  _et al._ (2019) A. Mietke, F. Jülicher, and I. F. Sbalzarini, “Self-organized shape dynamics of active surfaces,” Proc. Natl. Acad. Sci. U.S.A. 116, 29–34 (2019).
  * Messal  _et al._ (2019) H. A. Messal, S. Alt, R. M. M. Ferreira, C. Gribben, V. M.-Y. Wang, C. G. Cotoi, G. Salbreux, and A. Behrens, “Tissue curvature and apicobasal mechanical tension imbalance instruct cancer morphogenesis,” Nature 566, 126–130 (2019).
  * Helfrich (1973) W. Helfrich, “Elastic Properties of Lipid Bilayers: Theory and Possible Experiments,” Z. Für Naturforschung C 28, 693–703 (1973).
  * Helfrich and Prost (1988) W. Helfrich and J. Prost, “Intrinsic bending force in anisotropic membranes made of chiral molecules,” Phys. Rev. A 38, 3065–3068 (1988).
  * MacKintosh and Lubensky (1991) F. C. MacKintosh and T. C. Lubensky, “Orientational order, topology, and vesicle shapes,” Phys. Rev. Lett. 67, 1169–1172 (1991).
  * Salbreux and Jülicher (2017) G. Salbreux and F. Jülicher, “Mechanics of active surfaces,” Phys. Rev. E 96, 032404 (2017).
  * Deserno (2004) M. Deserno,  _Notes on Differential Geometry_ (2004).
  * Hug (2000) D. Hug,  _Measures, Curvatures and Currents in Convex Geometry_ (Albert-Ludwigs-Universität Freiburg, Germany, 2000).
  * Schneider and Weil (2008) R. Schneider and W. Weil,  _Stochastic and Integral Geometry_ (Springer Science & Business Media, Berlin, Germany, 2008).
  * Jensen and Kiderlen (2017) E. B. V. Jensen and M. Kiderlen, eds.,  _Tensor Valuations and Their Applications in Stochastic Geometry and Imaging_, Lecture Notes in Mathematics (Springer International Publishing, Cham, 2017).
  * Qingde Li and Griffiths (2004) Qingde Li and J. Griffiths, “Least squares ellipsoid specific fitting,” in  _Geometric Modeling and Processing, 2004. Proceedings_ (IEEE, Beijing, China, 2004).
  * Beale  _et al._ (2016) D. Beale, Y.-L. Yang, N. Campbell, D. Cosker, and P. Hall, “Fitting quadrics with a Bayesian prior,” Comp. Visual Media 2, 107–117 (2016).
  * Taubin (1995) G. Taubin, “Estimating the tensor of curvature of a surface from a polyhedral approximation,” Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit. , 902–907 (1995).
  * Hameiri and Shimshoni (2003) E. Hameiri and I. Shimshoni, “Estimating the principal curvatures and the darboux frame from real 3-D range data,” IEEE Trans. Syst., Man, Cybern. B 33, 626–637 (2003).
  * Ramakrishnan  _et al._ (2010) N. Ramakrishnan, P. B. Sunil Kumar, and J. H. Ipsen, “Monte Carlo simulations of fluid vesicles with in-plane orientational ordering,” Phys. Rev. E 81, 041922 (2010).
  * Crane (2020) K. Crane,  _Discrete Differential Geometry: An Applied Introduction_ (2020).
  * Mecke (2000) K. R. Mecke, “Additivity, Convexity, and Beyond: Applications of Minkowski Functionals in Statistical Physics,” in  _Statistical Physics and Spatial Statistics_ (Springer, Berlin, Germany, 2000).
  * Schröder-Turk  _et al._ (2013) G. E. Schröder-Turk, W. Mickel, S. C. Kapfer, F. M. Schaller, B. Breidenbach, D. Hug, and K. Mecke, “Minkowski tensors of anisotropic spatial structure,” New J. Phys. 15, 083028 (2013).
  * Hadwiger (1951) H. Hadwiger, “Vorlesungen über Integralgeometrie,” Abhdl. Math. Sem. Hamburg 17, 69 (1951).
  * Alesker (1999a) S. Alesker, “Description of Continuous Isometry Covariant Valuations on Convex Sets,” Geom. Dedicata 74, 241–248 (1999a).
  * Alesker (1999b) S. Alesker, “Continuous Rotation Invariant Valuations on Convex Sets,” Annals of Mathematics 149, 977 (1999b).
  * Schneider (2013) R. Schneider, “Local tensor valuations on convex polytopes,” Monatsh. Math. 171, 459–479 (2013).
  * Hug and Schneider (2014) D. Hug and R. Schneider, “Local Tensor Valuations,” Geom. Funct. Anal. 24, 1516–1564 (2014).
  * Schröder-Turk  _et al._ (2011) G. E. Schröder-Turk, W. Mickel, S. C. Kapfer, M. A. Klatt, F. M. Schaller, M. J. F. Hoffmann, N. Kleppmann, P. Armstrong, A. Inayat, D. Hug, M. Reichelsdorfer, W. Peukert, W. Schwieger, and K. Mecke, “Minkowski Tensor Shape Analysis of Cellular, Granular and Porous Structures,” Adv. Mater. 23, 2535–2553 (2011).
  * Bian  _et al._ (2020) X. Bian, S. Litvinov, and P. Koumoutsakos, “Bending models of lipid bilayer membranes: Spontaneous curvature and area-difference elasticity,” Comput. Methods Appl. Mech. Eng. 359, 112758 (2020).
  * Jülicher (1996) F. Jülicher, “The Morphology of Vesicles of Higher Topological Genus: Conformal Degeneracy and Conformal Modes,” J. Phys. II France 6, 1797–1824 (1996).
  * Beisbart  _et al._ (2006) C. Beisbart, M. S. Barbosa, H. Wagner, and L. d. F. Costa, “Extended morphometric analysis of neuronal cells with Minkowski valuations,” Eur. Phys. J. B 52, 531–546 (2006).
  * Beisbart  _et al._ (2002) C. Beisbart, R. Dahlke, K. Mecke, and H. Wagner, “Vector- and Tensor-Valued Descriptors for Spatial Patterns,” in  _Morphology of Condensed Matter_ (Springer, Berlin, Germany, 2002).
  * Kousholt and Kiderlen (2016) A. Kousholt and M. Kiderlen, “Reconstruction of convex bodies from surface tensors,” Adv. Appl. Math. 76, 1–33 (2016).
  * Kousholt (2017) A. Kousholt, “Reconstruction of n-dimensional convex bodies from surface tensors,” Adv. Appl. Math. 83, 115–144 (2017).
  * Hansen-Goos and Mecke (2009) H. Hansen-Goos and K. Mecke, “Fundamental Measure Theory for Inhomogeneous Fluids of Nonspherical Hard Particles,” Phys. Rev. Lett. 102, 018302 (2009).
  * Hansen-Goos and Mecke (2010) H. Hansen-Goos and K. Mecke, “Tensorial density functional theory for non-spherical hard-body fluids,” J. Phys.: Condens. Matter 22, 364107 (2010).
  * Capovilla and Guven (2002) R. Capovilla and J. Guven, “Stresses in lipid membranes,” J. Phys. A 35, 6233–6247 (2002).
  * Willmore (1982) T. Willmore,  _Total Curvature in Riemannian Geometry_ (Ellis Horwood, Chicester, England, 1982).
  * Lawson (1970) H. B. Lawson, “Complete Minimal Surfaces in S 3,” Ann. Math. 92, 335 (1970).
  * Lubensky and MacKintosh (1993) T. C. Lubensky and F. C. MacKintosh, “Theory of “Ripple” Phases of Lipid Bilayers,” Phys. Rev. Lett. 71, 1565–1568 (1993).
  * Frank and Kardar (2008) J. R. Frank and M. Kardar, “Defects in nematic membranes can buckle into pseudospheres,” Phys. Rev. E 77, 041705 (2008).
  * Uchida (2002) N. Uchida, “Dynamics of orientational ordering in fluid membranes,” Phys. Rev. E 66, 040902 (2002).
  * Morris and Rao (2019) R. G. Morris and M. Rao, “Active morphogenesis of epithelial monolayers,” Phys. Rev. E 100, 022413 (2019).
  * Matoz-Fernandez  _et al._ (2020) D. A. Matoz-Fernandez, F. A. Davidson, N. R. Stanley-Wall, and R. Sknepnek, “Wrinkle patterns in active viscoelastic thin sheets,” Phys. Rev. Research 2, 013165 (2020).
  * Merkel  _et al._ (2017) M. Merkel, R. Etournay, M. Popović, G. Salbreux, S. Eaton, and F. Jülicher, “Triangles bridge the scales: Quantifying cellular contributions to tissue deformation,” Phys. Rev. E 95, 032401 (2017).
  * Etournay  _et al._ (2015) R. Etournay, M. Popović, M. Merkel, A. Nandi, C. Blasse, B. Aigouy, H. Brandl, G. Myers, G. Salbreux, F. Jülicher, and S. Eaton, “Interplay of cell dynamics and epithelial tension during morphogenesis of the Drosophila pupal wing,” eLife 4, e07090 (2015).


long burn in time where agents must spend large amounts of time initially
learning the results of their actions, the ability to transfer knowledge
learned to alternative environments, the right schedule for exploring versus
exploiting a learned policy, and the ability to learn from small data sets
efficiently. These areas are being actively researched today, and will go a
long way towards making reinforcement learning more robust for industry
applications.

###  7.4 Game Theory and Multiple Agents

Agents are almost never operating in an environment by themselves. Agents are
simultaneously competing and collaborating with other agents and entities to
meet a variety of goals. The study of acting intelligently in a shared
environment with and against other agents that have shared and conflicting
goals is called game theory. Companies, organizations, customers, and many
others can all be thought of as agents operating in a shared environment with
shared and conflicting goals. In this section, we aim to briefly outline the
main aspects game theory, and discuss how to frame problems in a way that
allows us to algorithmically answer how to act intelligently in this kind of
environment. For more extensive discussions of the topic, we point the reader
to the following resources [8] [62] [67].

####  7.4.1 Games

Game theory represents the natural evolution and extension of decision theory.
While decision theory seeks to understand optimal decision making in an
environment where only the decision makers decisions must be evaluated, game
theory seeks to understand the optimal decision making strategy when there are
multiple decision makers with shared or conflicting goals. This is of prime
importance to the development of algorithms that can prescribe intelligent
actions because decisions are rarely made in and isolated environment with one
decision maker. Games are framed in an analogous way to the original decision
making problem framing. The main components of games in game theory are,

  * •

Multiple Agents: A game has multiple agents, each seeking to take actions that
lead to states with maximal utility.

  * •

Actions: Each agent has it’s own, unique set of actions that can be taken
given the state of the agent and the environment. These actions are also
typically restricted given a well defined set of rules that are associated
with the overall game. These actions may also be informed through
communication with other agents.

  * •

States: States in game theory typically have a shared and a non-shared
component depending on the kind of game. The shared component is typically a
set of number representing aspects of the environment fully or partially
visible to all agents. The non-shared component of a state may be the
particular context of an individual agent within an environment described by a
state. Depending on the visibility of an agent, the context of the other
agents may be known, and factored into the decision making process. Other
times, this is not known which adds to the complexity of reasoning about
decision making.

  * •

Rewards: Rewards in game theory are similar to those of decision theory. They
still represent scale based utility functions. However, these functions can
become much more complicated when the utility of states for a specific agent
factor in not only the utility for that agent, but also for the set of other
agents interacting with the environment. Sometimes utility values of a
specific agent’s state are directly related to the utility values of one or
more opposing or collaborating agents.

Game theory also makes various assumptions about the games and agents acting
in a game to simplify the mathematical reasoning and modeling. Most of these
assumptions follow logically from the definitions of rational behavior of goal
oriented agents in general. The main set of assumptions are,

  * •

Assumes games have utility scales associated with their outcomes.

  * •

Assumes a given choice of strategies by opponents, each player chooses the
strategy to optimize expected utility.

  * •

Any two games with the same game trees are considered the same games.

Given this general description of the problem, game theory concerns itself
mainly with two challenges. The first is the challenge of designing optimal
agents. Here the goal is to understand how to develop intelligent agents that
can learn strategies to maximize their expected utility throughout a game. The
second challenge is called mechanism design. Here the goal is to develop games
such that, given the game and the rules of the environment, the collective
good of all agent’s will be maximized by each agent maximizing their own
utility. Given our interest in understanding how to prescribe optional actions
in an environment, we focus on the first challenge.

Given the goal of acting optimally in games, game theory further breaks down
games into a wide variety of types. Given this breakdown, the goal is to
determine ways to construct an optimal policy for acting. Some of the ways to
categorize games include,

  * •

Number of Agents: Games with variable numbers of agents are typically broken
down into two player games, and n-player games. Most of the classical analysis
of games has been performed under the view of two player games. While
simplifying analysis and providing insights into the basic challenges that
arise in decision making due to the introduction of a single new agent, two
player game analysis is limited. N-player games while general are much more
challenging to analyze than two player games because of the exponentially more
complex dynamics at play.

  * •

Visibility: Games in which the environment is fully visible are referred to as
’perfect information’ games. Games in which the environment is only partially
observable are referred to as ’imperfect information’ games.

  * •

Chance: Games which have stochastic elements are referred to as games with
chance.

  * •

Competitive/Collaborative: Games in which agents are all competing with each
other and have conflicting goals are called ’strictly competitive’ games, and
games that allow for shared goals are called ’non-strictly competitive’ games.

Other more specific categorizations exist such as static, dynamic, zero-sum,
all exist. For a more comprehensive taxonomy of methods, we refer the reader
to [25]. A full treatment of all aspects of game theory are beyond the scope
of this paper. However, we will mention a few central ideas from game theory
that apply to acting intelligently in complex, multi-agent environments.

The first central idea is that finding a policy to act intelligently can still
be viewed as a form of search problem. Given the space of all possible
policies, the problem is to efficiently search for the policy that will
maximize expected utility, conditioned on the predicted behavior of all other
agents and the dynamics of the environment. Therefore, some of the methods for
solving certain games use optimization and modified search algorithms. One key
algorithm to solve simple two person games is the ’minimax’ search algorithm.
This search algorithm suggests that the action an agent should take in the
current state should be the one that minimizes the maximum expected utility of
the opponents optimal action from the next state. We defer the entire
implementation here [58], and simply specify that this is an example of a
search algorithm that is used to determine an agent’s optimal policy.

The second major idea is that of equilibrium. An equilibrium is a condition
under which each agent is assigned a policy such that no agent can benefit by
switching its policy. Equilibrium are effectively local optimum where the
dynamics of the environment due to the agents is fixed because no change in
agent policies can cause an increase in utility. John Nash proved that every
game has at least one equilibrium, which is typically called the Nash
equilibrium. Equilibrium are important because most game theorists agree that
a solution must necessarily be a Nash equilibrium. A closely related idea is
Pareto Optimality. An outcome to a game is considered Pareto Optimal if there
is no other outcome that all players would prefer. The idea of optimality and
equilibrium are often combined when attempting to evaluate proposed solutions
to specific kinds of games.

Another major idea is the idea of n-player games and the evolving complexity
of decision making with the addition of more agents. The representation and
modeling of the dynamics of complex games (particularly N-player games) is
extremely challenging, and partially related to the study of Complex Systems.
Complex systems are systems which are composed of fully or partially
independent phenomena that are all collectively interacting with each other.
Complex systems include finite automota, biological systems, multi-agent
games, and many other kinds of environments [10]. The field draws upon the
mathematics from a wide variety of other fields that seek to model highly non-
linear systems composed of multiple smaller components. We mention field of
complex systems as a promising area for solving problems within game theory,
and for understanding how to act intelligently in realistically complex
environments.

Finally, we mention the growing research field of multi-agent reinforcement
learning (MARL). MARL is still in its early research stages, but algorithms
developed have already been applied to games and game theoretic environments.
While lots of research still needs to be done to make MARL agents ready to be
applied to industry problems, they offer a promising advancement in the
understanding and modeling of intelligently acting in realistic environments.

There are many other challenges and ideas in game theory that we haven’t
discussed here. The field of game theory is of crucial importance to the goal
of mathematically and algorithmically determining how to act intelligently,
and presents a burgeoning set of methods that data scientists can apply to
industry problems to provide value.

####  7.4.2 Challenges

Game theory and the study of acting intelligently in complex environments with
multiple agents is still a highly active area of research today. While still
in the early stages of exploration, the field even today offers an interesting
solution to the value proposition of how to act intelligently in real world
environments. Most of the game theory’s advancements come from its application
to the study of economics and economics based games. Major challenges still
exist as barriers to the widespread application of game theory methods to real
world problems. Many of these challenges can be categorized as either
complexity based challenges, or modeling based challenges. Computational
complexity for finding optimal policies arise due to the exponential or even
factorial sizes of game trees, unexpected emergent behavior arising in the
game as a whole, constantly changing and non-stationary dynamics, and all of
the same difficulties that arose in standard Reinforcement Learning are
multiplied when those same algorithms are applied to multi-agent environments.
Challenges to formalization arise because games in game theory typically need
a fixed set of rules, a fixed set of actions, and a complex multi-variate
utility function to perform learning. Real world problems rarely have a fixed
set of rules and actions or an easy way to specify the utility of the world.
However, even with these major challenges, the value proposition of
determining how to act in mathematically optimal ways means this area of study
is likely to have ever increasing value for data scientists in industry
settings.

###  7.5 Example Applications

####  7.5.1 Online and Marketing Applications

Reinforcement learning has been applied to many areas of online interactions
with customers including search, recommendation, and advertising paper. The
paper by Zhao, Xia, Tang, and Yin [73] surveys the increasing role of
reinforcement learning for the aforementioned use cases. The authors discuss
how reinforcement learning is used for determining the content to return (this
is the action) given a search query, the ranking of query results, web page
layout optimization, and session based search optimization where one query is
temporally and contextually related to a previous query. Much of the work
proposes small extensions to reinforcement learning algorithms for the
specific application. All are formulated as markov decision processes, where
the reinforcement learning problem is slightly modified either in optimization
metric, or optimization algorithm that takes into account some of the unique
factors of the problem. Reinforcement learning is used in the recommendation
case for optimizing recommendations that maximize long term user engagement,
show diversity but also remain within a fixed ’theme’, and incorporate
temporal context. Online advertising challenges lay along a similar vein as
the search and content recommendation challenges. Recommending advertisements
that are grouped within a theme and most likely to create a purchase is one of
the most challenging and most valuable business propositions for the online
environment.

In this 2015 paper [68], Theocharous, Thomas, and Ghavamzadeh outline an
algorithm developed to recommend ads that maximize ’life-time value’, or the
total number of times a user will click over multiple visits. Their paper
outlines a method for computing good policies in a scalable way, and a set of
off-policy techniques to check the validity of a policy before it is deployed
to a real environment. They implement Q-learning by estimating state action
values using a random forest regression model, weighting values using
importance sampling, and set Q values using a bounding mechanism. They show
that reinforcement learning can be an effective method for recommending ads
that optimize long term metrics such as life-time value. These are just a few
of the many ways in which reinforcement learning is being applied to online
and marketing business problems.

####  7.5.2 Customer Experience and Engagement

Reinforcement learning is also being increasingly applied to improve
engagements with customers for reasons other than marketing or ads. Cisco uses
reinforcement learning algorithms in the optimization of its digital
engagements with customers. Adobe’s marketing methods can also transfer to the
case of engagements with customers for the purpose of improving a customer’s
experience. This paper [61] describes how they use online click stream data
along with reinforcement learning to generate a value function, which can be
used as surrogate for measuring customer satisfaction. They argue that the
proxy based approach to measuring customer satisfaction using reinforcement
learning allows companies to do three things that can’t be done with standard
surveys. The first is to generate ratings by the use of nearly all customers,
whereas surveys are typically information gathered from only a few customer.
The second is that surveys are not tied to any kind of actual customer
behavior. If a customer interacts with lots of online content or uses their
device a lot more, this specific behavior cannot be tied back to their
responses on a survey. The third is that reinforcement learning can show the
value of particular sequences of click actions, while surveys naturally roll
up all previous customer sentiment into one point. The interesting value
proposition from this paper is that reinforcement learning is used as a tool
for understanding user behavior and customer experiences, going beyond just
the use of reinforcement learning as an automated method for taking actions.

Another paper [14] discusses the case of customer experience management from
the use case of digital programmable networks. They discuss how toady’s
adaptive programmable networks operate through the use of traditional rule-
based decision making, and how this may lead to less than optimal customer
experiences on the network. The authors do a good job of positioning both the
technical and the business use case for developing algorithms that can
automatically prescribe actions and make decisions that will increase metrics
such as ’Quality-of-Service’ and ’Quality-of-Experience’. The authors
introduce the idea of using an iterative learning algorithm like reinforcement
learning to attempt to solve the high level challenges with automatically
making technical decisions that are likely to improve customers experiences
when interacting with a network.

A third paper [24] describes how Microsoft used contextual bandits in an
attempt to implement a customer support chat bot for goal directed dialogue.
They utilized reinforcement learning to handle the diversity of customer
intents that exist, the variety of Microsoft products, and the combination of
both customer and business value metrics such as solving a customer problem
and increasing content click rate. Their work borrows from other applications
of reinforcement learning in the wild, and brings up some important points
about reinforcement learning systems that need to learn to interact in an
environment where no simulation exists initially. Some of their guiding
insights include starting with simplified action spaces, starting with
imitation learning, intelligent exploration, support for environment changes,
and a focus on reward shaping to name a few.

The customer experience and engagement case uncovers a major challenge of
applying reinforcement learning. It is challenging to discover the best
algorithm, parameters, and hyper-parameters to use for a specific engagement
case. It is common practice to use simulation environments to tune and
customize the parameters for nearly all reinforcement learning methods.
However a simulation environment typically doesn’t exist out of the box. A way
around this may be to attempt to create a simulated environment, but that can
be challenging because we often already assume the model of the environment is
too challenging to model, which is why we used the model free methods of RL in
the first place. In application, modern reinforcement learning creates a
chicken and egg problem. We need to pick an algorithm, parameters, and hyper-
parameters to apply reinforcement learning to an environment. However, we must
test and check our choices to determine the right combination for a particular
environment which leads to excessive testing. This can actually preclude
certain methods like deep reinforcement learning when a simulated environment
doesn’t exist to modify learning algorithms.

##  8 Case Study

In this section, we outline a common problem encountered by businesses, and
how to apply our methodology to maximize the value returned from this data
science project. Assume that a given company provides a service to customers
for which they must pay a monthly subscription. A well known question in this
domain is whether or not a particular customer or group of customers will
renew their subscription or not. This is typically referred to as predicting
churn, or churn risk prediction. It is closely related to sales and demand
forecasting mentioned in the prediction and pattern mining section. Many data
science projects focus simply on answering this question. We show below that
value can be gained beyond simply constructing a model that predicts the churn
risk of a customer by framing the problem in terms of the big three questions.

###  8.1 Question 1: What is happening or will happen?

As we’ve already stated, the main question that is typically posed to a team
of data scientists is ’Can we accurately predict which customers will renew
and which ones won’t?’ While this is a primary question asked by the business,
there are many other questions that fall into the area of prediction and
pattern mining including,

  1. 1.

How much revenue can we expect from renewals? What does the distribution look
like?

  2. 2.

What’s the upper/lower bound on the expected revenue predicted by the models?

  3. 3.

What are the similar attributes among customers likely to churn versus not
churn?

  4. 4.

What are the descriptive statistics for customers likely to churn vs not churn
collectively, in each label grouping, and in each unsupervised grouping?

Each of the above questions can be answered systematically by framing them as
problems either in prediction or pattern mining, and by using the wide variety
of mathematical methods found in the referenced materials from that section.
These are the questions and methods data scientists are most familiar, and
will most commonly be answered for a business.

###  8.2 Question 2: Why is this happening or going to happen?

Given this first question, the immediate next question is why. Why are
customers likely or not likely to churn? For each question that we can build a
model for, we can also perform a causal analysis. Thus, we can already
potentially double the value that a data science project returns by simply
adding on a causal analysis to each predictive model built. It’s important to
bring up again that this question is so important that most data scientists
are either answering it incorrectly, or are misrepresenting the information
from statistical associations.

Specifically, when a data scientist is asked the question of why a customer is
likely to churn, they almost exclusively turn to feature importance and local
models such as LIME [55] and others. These methods for describing the reason
for a prediction are almost always incorrect because there is a disconnect
between what the business stakeholder is asking for and what the data
scientist is providing because of two different interpretations of the term
’why’. Technically, one can argue that feature importance measures what
features are important to ’why’ a model makes a prediction, and this would be
correct. However, a business stakeholder is nearly always asking the question
of ’what is causing the metric itself’ and not ’what is causing the metric
prediction’. The business stakeholder wants to know the causal mechanisms for
why a metric is a particular number. This is something that feature importance
absolutely does not answer. The stakeholder wants to use the understanding of
the causal mechanisms to take an action to change the prediction to be more in
their favor. This requires a causal analysis. However, most data scientists
simply take the features with highest measured importance and present them to
the stakeholder as though they are the answer to their causal question. This
is objectively wrong, yet is time and again presented to stakeholders by
seasoned statisticians and data scientists.

The issue is compounded by the further confusion added by discussions around
’interpretable models’ and by the descriptions of feature importance analysis.
LIME goes so far as to say that is ’explaining what machine learning
classifiers (or models) are doing’. While still a technically correct
statement, these methods are being used to incorrectly answer causal
questions, leading stakeholders to take actions that can often have the
opposite effect of what they intended.

While we’ve outlined the main causal question, there are a number of questions
that can also be asked, and corresponding analysis that can be performed
including,

  1. 1.

How are variables correlated with each other and the churn label? (A non-
causal question)

  2. 2.

What are the important features for prediction in a model in general? (A non-
causal question)

  3. 3.

What are the most important features for prediction for an individual? Do
groupings of customers with locally similar relationships exist? (A non-causal
question)

  4. 4.

What are the possible confounding variables? (A causal question)

  5. 5.

After controlling for confounding variables, how do the predictions change? (A
non-causal question benefiting from causal methods)

  6. 6.

What does the causal bayes net structure look like? What are all of the
reasonable structures? (A causal question)

  7. 7.

What are the causal effect estimates between variables? What about between
variables and the class label? (A causal question).

Many of these questions can be answered in whole or in part by a thorough
causal analysis using the methods we outlined in the corresponding causal
inference section, and further multiply the value returned by a particular
data science project.

###  8.3 Question 3: How can we take action to make the metrics improve?

The third question to answer is ’what actions can a stakeholder take to
prevent churn?’ This is ultimately the most valueable of the three questions.
The first two question set the context for who to focus on and where to focus
efforts. Answering this question provides stakeholders with a directed and
statistically valid means to improve the metrics they care about given complex
environments. While still challenging given the methods available today (and
presented in the section on intelligent action), it provides one of the
greatest value opportunities. Some other questions that can be answered
related to intelligent action that stakeholders may be interested in include,

  1. 1.

What variables are likely to reduce churn risk if our actions could influence
them?

  2. 2.

What actions have the strongest impact on the variables that are likely to
influence churn risk, or to reduce churn risk directly?

  3. 3.

What are the important pieces of contextual information relevant for taking an
action?

  4. 4.

What are the new actions that should be developed and tested in an attempt to
influence churn risk?

  5. 5.

What actions are counter-productive or negatively impact churn risk?

  6. 6.

What does the diminishing marginal utility of an action look like? At what
point should an action no longer be taken?

The right method to use for prescribing intelligent action depends largely on
the problem and the environment. If the environment is complex, the risks
high, and there is not much chance for an automated system to be implemented,
then something like decision theory, influence diagrams, and game theory based
analysis are good candidates. However, if a problem and stakeholder are open
to the use of an automated agent to learn and prescribe intelligent actions,
then reinforcement learning may be a good choice. While possibly the most
valuable of the big three questions to answer, it also exists as one of the
most challenging. There still many open research questions related to
answering this question, but the value proposition means that it’s likely an
area that will see increased industry investment in the coming years.

##  9 Conclusions

The value of answering causal and intelligent action questions are real. They
are so important that current predictive and statistical methodologies are
miss-applied in an attempt to answer them. They are the natural next step
after any supervised or unsupervised prediction or pattern mining task. The
work that has been done in the fields of causal inference, decision theory,
and reinforcement learning will begin surface as a part of every project in
business applications. It’s likely that in the future, every data science
modeling or analysis project will have 3 distinct phases, the first being to
build a classification or regression model, the second being to perform causal
inference, and the third being to build decision networks, reinforcement
learning algorithms, or some other intelligent action system to improve
desired metrics.

The barriers to achieving this inevitable result include a lack of knowledge
and teaching about the methodologies, and possibly limitations with the
technologies and methods themselves. Causal inference requires an uncommon
statistical understanding, decision theory and influence diagrams require
formalizing a problem’s structure, and reinforcement learning is currently
sensitive to the problem formulation and can be data intensive. However,
regardless of the specific methods that are available today, the need to
understand and answer ’The Big Three’ questions in mathematically optimal ways
will remain. The value that Data Science, Machine Learning, and AI can bring
in the coming years will be in developing methods to provide answers to these
value oriented questions, as well as in the automation of the processes to
produce these answers.

## References

  * [1] Causal ml: A python package for uplift modeling and causal inference with ml.  https://github.com/uber/causalml, 2019.
  * [2] Dowhy: A python package for causal inference.  https://github.com/microsoft/dowhy, 2019.
  * [3] Charu C. Aggarwal.  Data Mining: The Textbook.  Springer Publishing Company, Incorporated, 2015.
  * [4] Charu C. Aggarwal.  Recommender Systems: The Textbook.  Springer Publishing Company, Incorporated, 1st edition, 2016.
  * [5] Andrew G. Barto and Thomas G. Dietterich.  Reinforcement Learning and Its Relationship to Supervised Learning, chapter 2, pages 45–63.  John Wiley & Sons, Ltd, 2012.
  * [6] Marianne Bertrand, Esther Duflo, and Sendhil Mullainathan.  How Much Should We Trust Differences-In-Differences Estimates?*.  The Quarterly Journal of Economics, 119(1):249–275, 02 2004.
  * [7] Christopher M. Bishop.  Pattern Recognition and Machine Learning (Information Science and Statistics).  Springer-Verlag, Berlin, Heidelberg, 2006.
  * [8] Giacomo Bonanno.  Game Theory.  North Charleston, SC, USA: CreateSpace Independent Publishing Platform, 2018.
  * [9] Ruichu Cai, Mei Liu, Yong Hu, Brittany L. Melton, Michael E. Matheny, Hua Xu, Lian Duan, and Lemuel R. Waitman.  Identification of adverse drug-drug interactions through causal association rule discovery from spontaneous adverse event reports.  Artificial Intelligence in Medicine, 76:7 – 15, 2017.
  * [10] Brian Castellani.  Map of the complexity sciences.  Art & Science Factory, 2018.
  * [11] David Chickering.  Optimal structure identification with greedy search.  Journal of Machine Learning Research, 3:507–554, 01 2002.
  * [12] David Maxwell Chickering.  Learning equivalence classes of bayesian-network structures.  Journal of machine learning research, 2(Feb):445–498, 2002.
  * [13] Gabriel Dulac-Arnold, Daniel J. Mankowitz, and Todd Hester.  Challenges of real-world reinforcement learning.  CoRR, abs/1904.12901, 2019.
  * [14] Haris Gacanin and Mark Wagner.  Artificial intelligence paradigm for customer experience management in next-generation networks: Challenges and perspectives.  CoRR, abs/1805.06254, 2018.
  * [15] Clark Glymour, Kun Zhang, and Peter Spirtes.  Review of causal discovery methods based on graphical models.  Frontiers in Genetics, 10:524, 2019.
  * [16] Carlos A. Gomez-Uribe and Neil Hunt.  The netflix recommender system: Algorithms, business value, and innovation.  ACM Trans. Manage. Inf. Syst., 6(4):13:1–13:19, December 2015.
  * [17] Sander Greenland, Michael P. Fay, Erica H. Brittain, Joanna H. Shih, Dean A. Follmann, Erin E. Gabriel, and James M. Robins.  On causal inferences for personalized medicine: How hidden causal assumptions led to erroneous causal claims about the d-value.  The American Statistician, 0(0):1–13, 2019.
  * [18] Justin Grimmer.  We are all social scientists now: How big data, machine learning, and causal inference work together.  PS: Political Science &; Politics, 48(1):80–83, 2015.
  * [19] Totte Harinen and Bonnie Li.  Using causal inference to improve the uber user experience, Jun 2019.
  * [20] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The Elements of Statistical Learning.  Springer Series in Statistics. Springer New York Inc., New York, NY, USA, 2001.
  * [21] Paul W. Holland.  Statistics and causal inference.  Journal of the American Statistical Association, 81(396):945–960, 1986.
  * [22] J. E. Matheson Howard, R. A.  Influence diagrams.  Readings on the Principles and Applications of Decision Analysis, 2(1):719–762, 1984.
  * [23] Guido W. Imbens and Donald B. Rubin.  Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction.  Cambridge University Press, 2015.
  * [24] Nikos Karampatziakis, Sebastian Kochman, Jade Huang, Paul Mineiro, Kathy Osborne, and Weizhu Chen.  Lessons from real-world reinforcement learning in a customer support bot.  CoRR, abs/1905.02219, 2019.
  * [25] Anthony Kelly.  Decision making using game theory: An introduction for managers.  Decision Making Using Game Theory: An Introduction for Managers, 01 2003.
  * [26] Samantha Kleinberg and George Hripcsak.  A review of causal inference for biomedical informatics.  Journal of Biomedical Informatics, 44(6):1102 – 1112, 2011.
  * [27] Michael Lechner.  The Estimation of Causal Effects by Difference-in-Difference Methods.  University of St. Gallen Department of Economics working paper series 2010 2010-28, Department of Economics, University of St. Gallen, September 2010\.
  * [28] Bonnie Li and Totte Harinen.  Mediation modeling at uber: Understanding why product changes work (and don’t work), Oct 2019.
  * [29] Yelin Li, Junjie Wu, and Hui Bu.  When quantitative trading meets machine learning: A pilot survey.  pages 1–6, 06 2016.
  * [30] Victor Lo.  The true lift model - a novel data mining approach to response modeling in database marketing.  SIGKDD Explorations, 4:78–86, 01 2002.
  * [31] Daniel Malinsky and David Danks.  Causal discovery algorithms: A practical guide.  Philosophy Compass, 13(1):e12470, 2018.  e12470 10.1111/phc3.12470.
  * [32] Colin McFarland, Michael Pow, and Julia Glick.  Quasi experimentation at netflix, May 2019.
  * [33] Thomas M. Mitchell.  Machine Learning.  McGraw-Hill, Inc., New York, NY, USA, 1 edition, 1997.
  * [34] Jovana Mitrovic, Dino Sejdinovic, and Yee Whye Teh.  Causal inference via kernel deviance measures.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 6986–6994. Curran Associates, Inc., 2018.
  * [35] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A. Rusu, Joel Veness, Marc G. Bellemare, Alex Graves, Martin Riedmiller, Andreas K. Fidjeland, Georg Ostrovski, Stig Petersen, Charles Beattie, Amir Sadik, Ioannis Antonoglou, Helen King, Dharshan Kumaran, Daan Wierstra, Shane Legg, and Demis Hassabis.  Human-level control through deep reinforcement learning.  Nature, 518(7540):529–533, February 2015.
  * [36] Ricardo Pio Monti, Kun Zhang, and Aapo Hyvärinen.  In UAI, 2019.
  * [37] Stephen L. Morgan and Christopher Winship.  Counterfactuals and Causal Inference: Methods and Principles for Social Research.  Analytical Methods for Social Research. Cambridge University Press, 2 edition, 2014.
  * [38] Kevin P. Murphy.  Machine Learning: A Probabilistic Perspective.  The MIT Press, 2012.
  * [39] Jorge Nocedal and Stephen J. Wright.  Numerical optimization.  Springer, 2006.
  * [40] Elizabeth L. Ogburn, Oleg Sofrygin, Ivan Diaz, and Mark J. van der Laan.  Causal inference for social network data, 2017.
  * [41] Nuno Oliveira, Paulo Cortez, and Nelson Areal.  The impact of microblogging data for stock market prediction: Using twitter to predict returns, volatility, trading volume and survey sentiment indices.  Expert Systems with Applications, 73, 12 2016.
  * [42] Barbara Osimani and Fiorenzo Mignini.  Causal assessment of pharmaceutical treatments: Why standards of evidence should not be the same for benefits and harms?  Drug Safety, 38, 12 2014.
  * [43] J. Pearl.  Bayesian Networks: A Model of Self-activated Memory for Evidential Reasoning.  Report. UCLA, Computer Science Department, 1985.
  * [44] J. Pearl, M. Glymour, and N.P. Jewell.  Causal Inference in Statistics: A Primer.  Wiley, 2016.
  * [45] Judea Pearl.  Causality: Models, Reasoning and Inference.  Cambridge University Press, New York, NY, USA, 2nd edition, 2009.
  * [46] Judea Pearl.  The foundations of causal inference.  Sociological Methodology, 40(1):75–149, 2010.
  * [47] Judea Pearl and Dana Mackenzie.  The Book of Why: The New Science of Cause and Effect.  Basic Books, Inc., New York, NY, USA, 1st edition, 2018.
  * [48] Erik Brynjolfsson Ryan Calo Oren Etzioni Greg Hager Julia Hirschberg Shivaram Kalyanakrishnan Ece Kamar Sarit Kraus Kevin Leyton-Brown David Parkes William Press AnnaLee Saxenian Julie Shah Milind Tambe Astro Teller Peter Stone, Rodney Brooks.  Artificial intelligence and life in 2030.  2016\.
  * [49] Jonas Peters, Dominik Janzing, and Bernhard Schlkopf.  Elements of Causal Inference: Foundations and Learning Algorithms.  The MIT Press, 2017.
  * [50] Bilal Qureshi, Faisal Kamiran, Asim Karim, Salvatore Ruggieri, and Dino Pedreschi.  Causal inference for social discrimination reasoning.  Journal of Intelligent Information Systems, 11 2019.
  * [51] Nicholas J. Radcliffe.  Using control groups to target on predicted lift: Building and assessing uplift model.  2007\.
  * [52] Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romero, and Clark Glymour.  A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images.  International journal of data science and analytics, 3(2):121–129, 2017.
  * [53] Microsoft Research.  Econml: A python package for ml-based heterogeneous treatment effects estimation.  https://github.com/microsoft/EconML, 2019.  Version 0.x.
  * [54] Michael D. Resnik.  Choices: An Introduction to Decision Theory.  University of Minnesota Press, ned - new edition edition, 1987.
  * [55] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin.  "why should I trust you?": Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016, pages 1135–1144, 2016.
  * [56] Paul R. Rosenbaum and Donald B. Rubin.  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70(1):41–55, 1983.
  * [57] Kenneth J. Rothman and Sander Greenland.  Causation and causal inference in epidemiology.  American Journal of Public Health, 95(S1):S144–S150, 2005.  PMID: 16030331.
  * [58] Stuart Russell and Peter Norvig.  Artificial Intelligence: A Modern Approach.  Prentice Hall Press, Upper Saddle River, NJ, USA, 3rd edition, 2009.
  * [59] Danni Ruthvan.  An illustration of the components of the Difference in Differences method for measuring causal impact. Self Published, Jun 2014.
  * [60] Dev Shah, Haruna Isah, and Farhana Zulkernine.  Stock market analysis: A review and taxonomy of prediction techniques.  International Journal of Financial Studies, 7(2):26, May 2019.
  * [61] Atanu R. Sinha, Deepali Jain, Nikhil Sheoran, Sopan Khosla, and Reshmi Sasidharan.  Surveys without questions: A reinforcement learning approach.  In AAAI, 2019.
  * [62] W. Spaniel.  Game Theory 101: the Complete Textbook.  CreateSpace Independent Publishing Platform, 2011.
  * [63] Peter Spirtes and Clark Glymour.  An algorithm for fast recovery of sparse causal graphs.  Social Science Computer Review - SOC SCI COMPUT REV, 9:62–72, 04 1991.
  * [64] Peter Spirtes, Clark Glymour, Richard Scheines, Stuart Kauffman, Valerio Aimale, and Frank Wimberly.  Constructing bayesian network models of gene expression networks from microarray data.  Proc. of the Atlantic Symposium on Computational Biology, Genome Information Systems & Technology, 02 2002.
  * [65] Eric Strobl, Kun Zhang, and Shyam Visweswaran.  Approximate kernel-based conditional independence tests for fast non-parametric causal discovery.  Journal of Causal Inference, 7, 02 2017.
  * [66] Richard S. Sutton and Andrew G. Barto.  Reinforcement Learning: An Introduction.  The MIT Press, second edition, 2018.
  * [67] S. Tadelis.  Game Theory: An Introduction.  Princeton University Press, 2013.
  * [68] Georgios Theocharous, Philip S. Thomas, and Mohammad Ghavamzadeh.  Personalized ad recommendation systems for life-time value optimization with guarantees.  In Proceedings of the 24th International Conference on Artificial Intelligence, IJCAI’15, pages 1806–1812. AAAI Press, 2015.
  * [69] Donald L. Thistlethwaite and D T Campbell.  Regression-discontinuity analysis: An alternative to the ex post facto experiment.  1960\.
  * [70] Sebastian Thrun and Anton Schwartz.  Issues in using function approximation for reinforcement learning.  1993\.
  * [71] Juan Pablo Usuga Cadavid, Samir Lamouri, and Bernard Grabot.  Trends in machine learning applied to demand & sales forecasting: A review.  In International Conference on Information Systems, Logistics and Supply Chain, Lyon, France, July 2018.
  * [72] Larry Wasserman.  All of Statistics: A Concise Course in Statistical Inference.  Springer Publishing Company, Incorporated, 2010.
  * [73] Xiangyu Zhao, Long Xia, Jiliang Tang, and Dawei Yin.  "deep reinforcement learning for search, recommendation, and online advertising: A survey" by xiangyu zhao, long xia, jiliang tang, and dawei yin with martin vesely as coordinator.  SIGWEB Newsl., (Spring):4:1–4:15, July 2019.