# Role of spatial patterns in fracture of disordered multiphase materials

Rajat Pratap Singh Parihar<EMAIL_ADDRESS>Department of Applied
Mechanics, Indian Institute of Technology Madras, Chennai-600036, India
Dhiwakar V. Mani  Department of Applied Mechanics, Indian Institute of
Technology Madras, Chennai-600036, India  Anuradha Banerjee<EMAIL_ADDRESS>Department of Applied Mechanics, Indian Institute of Technology Madras,
Chennai-600036, India  R.Rajesh<EMAIL_ADDRESS>The Institute of
Mathematical Sciences, C.I.T. Campus, Taramani, Chennai-600113, India  Homi
Bhabha National Institute, Training School Complex, Anushakti Nagar,
Mumbai-400094, India

###### Abstract

Multi-phase materials, such as composite materials, exhibit multiple competing
failure mechanisms during the growth of a macroscopic defect. For the
simulation of the overall fracture process in such materials, we develop a
two-phase spring network model that accounts for the architecture between the
different components as well as the respective disorders in their failure
characteristics. In the specific case of a plain weave architecture, we show
that any offset between the layers reduces the delocalization of the stresses
at the crack tip and thereby substantially lowers the strength and fracture
toughness of the overall laminate. The avalanche statistics of the broken
springs do not show a distinguishable dependence on the offsets between
layers. The power law exponents are found to be much smaller than that of
disordered spring network models in the absence of a crack. A discussion is
developed on the possibility of the avalanche statistics being those near
breakdown.

##  I Introduction

In biological materials, nature’s design exploits the beneficial aspects of
combining different materials and the spatial patterns between them to develop
lightweight material systems with significantly improved resistance to
fracture [1, 2, 3]. Synthetic or engineered composite materials, similarly,
have aimed to expand the strength to weight as well as toughness to weight
performances of available materials. Different spatial patterns, such as short
stiff fiber reinforcements randomly arranged in a compliant matrix, long
fibers unidirectionally arranged, or in a two- or three-dimensional woven
architecture, etc., have been developed successfully [4, 5, 6, 7, 8, 9]. The
more intricate the architecture, the higher the associated complexity of the
deformation and in particular, the fracture process, due to several competing
failure mechanisms, making predictive modeling significantly more challenging.
For better design and life assessment of heterogeneous composite materials, in
the present study, we develop a statistical approach to examine some of the
factors that influence the defect tolerance of spatially patterned
heterogeneous composite materials.

The effect of heterogeneity on the macroscopic mechanical performance of
composite materials, measured in terms of stiffness, strength, and toughness,
has been examined by several researchers in both statistical physics as well
as engineering. Depending on the relevant length scales of a given system,
many different theoretical frameworks have been adopted for analysis, ranging
from atomistic simulations for design of materials to continuum models for
design of components and structures [10, 11, 12, 13, 14, 15]. Continuum
description of composite material typically involves micro-mechanical modeling
of a representative volume to arrive at the effective anisotropic elastic and
plastic behavior. The homogenized anisotropic constitutive relations are then
combined with a damage initiation and evolution criterion for simulation of
progressive damage [16, 17, 18]. While these models have been shown to
accurately describe uniaxial responses in particular directions, the material
behavior under multiaxial loading is, however, not well reproduced. [16, 17,
18].

Brittle heterogeneous materials, when subject to mechanical loads, are prone
to micro-cracking at multiple sites. On the application of higher stresses,
these micro-cracks tend to interact and grow, resulting in multiple events
prior to final failure, as evident in the acoustic emission activity [19, 20].
Such complex fracture processes, which typically occur over large process
zones, are stochastic in nature, and have been investigated comprehensively
using lattice-based models for single-phase systems with disorder. In these
lattice-based models, the continuum is discretized into a collection of
interconnected particles, and a simple inter particle interaction is defined
to account for the mechanical and fracture properties of the material [21, 22,
23]. The scaling properties observed in acoustic emission experimental data
during fracture of heterogeneous materials have been of significant interest
in the field of statistical physics. Attempts have been made to interpret the
scaling laws of energy release and fluctuations in temporal statistics via
concepts of criticality and phase transitions [24, 22, 19, 25]. Qualitative
features of the role of disorder in crack paths, and macroscopic response have
also been investigated in detail because arbitrarily oriented, multiple crack
paths, typically observed in materials like concrete, particulate composites,
ceramics, cortical bones, etc., are easy to simulate using lattice-based
models [26, 27, 28, 29, 30, 31].

Brittle multiphase heterogeneous systems are comparatively less well studied.
Investigating the role of elastic heterogeneity in a two-phase system, studied
using a network of hard and soft springs, Urabe et al. [32] showed the
possibility of designing material systems tougher and stronger than individual
components. The two-phase network was also utilized in the translaminar
fracture simulation, under mode I and mixed mode, of woven composites [33]. In
these studies, while the relative proportion of hard and soft bonds and their
orientations was maintained to ensure equivalence with the macroscopic elastic
modulus of the composite material, the details of the meso-structural patterns
of the woven reinforcements were ignored, as the heterogeneity was randomly
distributed in the domain.

The role of spatial patterns in fracture of heterogeneous systems was analyzed
in a series of investigations by Buehler and co-workers [34, 35, 36]. Spatial
patterning in a two-dimensional heterogeneous spring bead network was
represented by introducing a correlation length in the spatial distribution of
the Young’s moduli modeled as a Gaussian process, such that any significant
variation in the modulus occurs only over the correlation length. Toughening
in the heterogeneous solid with disordered elastic fields was attributed to
the “distribution-of-weakness” mechanism, resulting from the observed crack
arrest and stress delocalization [34]. Further, some of these features could
be reproduced in solvable models of rectangular blocks, rods, and plates [35,
36].

How does the fracture process differ between laminates that are elastically
equivalent but differ in the spatial patterns between their hard and soft
phases? Under what relative fracture properties of the phases does the elastic
heterogeneity result in enhancement of toughness? Do avalanche statistics
depend on spatial patterning? In this paper, we explore the role of meso-
structural patterns by considering a representative system of woven fiber-
reinforced laminate. To incorporate the composite’s mesoscopic spatial
pattern, we develop a multiphase two-dimensional computational model using
spring networks, the elastic properties of each phase being estimated using
standard homogenization technique and the rule of mixture. Within the model,
we study in detail the dependence of the elastic stress distribution and
fracture properties such as strength, toughness, crack paths, and avalanche
statistics on the spatial patterning of the laminate. In particular, we show
that toughness is enhanced with increased elastic heterogeneity.

##  II Modeling and Simulation

In this section, we develop a two-dimensional model for the analysis of plain
weave laminates. Four steps are involved. First, the three-dimensional plain
weave laminate is mapped onto a two-dimensional geometric pattern based on the
local volume fraction of the yarns (Sec. II.1). Second, from the two-
dimensional geometric pattern, a discrete element model in the form of a two-
dimensional spring network model made up of bonds of multiple types is
constructed (Sec. II.2). Third, a methodology of assigning local elastic
parameters to springs is developed so that the experimentally observed
macroscopic elastic behavior is reproduced (Sec. II.3). Fourth, the algorithm
for simulating the resultant spring network model using molecular dynamics is
discussed (Sec. II.4).

###  II.1  Geometric Modeling

The representative material system being modeled is a soft matrix embedded
with a harder reinforcement with a plain weave architecture. We focus on the
most common reinforcement patterns in which yarns of fibers are interlaced in
a textile-like geometric pattern. A typical pattern of a single layer or
lamina is shown in Fig. 1a. The pattern consists of two types of fibers: those
oriented in the x-direction denoted as $x$-fibers and those oriented in the
orthogonal $y$-direction denoted as $y$-fibers. When sectioned in the
thickness direction, a well defined geometric pattern emerges where, as shown
in Fig. 1b, the cross-section consists of $y$-fibers appearing as isolated
elliptic regions while the $x$-fibers are continuous and have an oscillatory
pattern.

Figure 1: Schematic diagrams of (a) a single lamina or layer with fibers woven
in perpendicular directions, (b) cross-sectional view of the lamina, at the
sectioning plane (shown in red [translucent cutting plane] in (a)), (c) a
laminate in which the constituent laminae are stacked randomly, and (d) cross-
sectional view of the laminate at the sectioning plane.

Typical composite laminates consist of 6 to 8 such layers in the thickness
direction ($z$-direction), with each layer having possible offsets in the $xy$
plane (Fig. 1c). Cross-sections of these laminates have periodically patterned
rows that are phase shifted depending on the relative positioning of the
layers with respect to the sectioning plane, as shown in Fig. 1d. Note that
depending on the relative positioning of any lamina with respect to the
sectioning plane, the cross-sectional thickness of the $x$-fiber may vary from
its maximum thickness (first layer in Fig. 1d) to zero thickness (last layer
in Fig. 1d, where the grey tube is missing).

For a typical single lamina, as well as the laminate, as shown in Fig. 1,
based on the pattern observed in the given cross-section, the volume fractions
of $x$\- and $y$-fibers along the section’s length are shown in Fig. 2a-2b. In
the cross-section of a single lamina, the volume fraction of the $y$-fibers
varies between its maximum at the center of the elliptic cross-section and
zero in the region between the ellipses as seen in Fig. 2a. Correspondingly,
the $x$-fibers have a well-defined pattern as well. In a laminate, as shown in
Fig. 2b, the volume fraction of both $y$-fibers as well as $x$-fibers when
averaged over the layers has comparatively less variation between the maximum,
$\phi_{\rm max}$, and minimum, $\phi_{\rm min}$.

(a)  (b)

Figure 2: Variation of the volume fractions of x-fibers and y-fibers at a
cross-section of (a) a single lamina and (b) a representative laminate of six
layers. The horizontal dotted line $\phi_{c}$ is a threshold as defined in
Table 1.

A representative area element of a single lamina is shown in Fig. 3a. It is
composed of 4 interconnected regions: regions which have no fibers (referred
as R1 and represented by the color green [grey] in Fig. 3b), regions which
have fibers only in the $x$-direction (R2, red [dark grey]), regions which
have fibers only in the $y$-direction (R3, yellow [light grey]) and regions
which have fibers in both directions (R4, blue [black]). The distinction
between the fiber-rich regions and the matrix-rich (no fibers) regions becomes
less apparent when several layers are stacked together in a laminate. For a
simplified 2-dimensional representation of a laminate, the spatial
distribution of volume fraction of fibers can be discretized by introducing a
threshold cut-off $\phi_{c}$, such that it is approximated to be composed of
fiber-rich and matrix-rich regions. For the present study, the criteria for
allocating different regions based on the threshold and the volume fractions
of $x$-fibers and $y$-fibers, was taken to be as per the rule in Table 1.

Figure 3: (a) A representative area element of a single lamina. (b) The corresponding allocation of regions based on the configuration of $x$-fibers, $y$-fibers and matrix.  Table 1: Criteria for allocating different regions based on the threshold, $\phi_{c}$ and the volume fractions of $x$-fibers and $y$-fibers  Regions |  $\phi$ (x-fiber) |  $\phi$ (y-fiber)
---|---|---
Region rich in $y$-fiber and $x$-fiber (R4) |  $\geq\phi_{c}$ |  $\geq\phi_{c}$
$y$-fiber rich region (R3) |  $<\phi_{c}$ |  $\geq\phi_{c}$
$x$-fiber rich region (R2) |  $\geq\phi_{c}$ |  $<\phi_{c}$
Matrix rich region (R1) |  $<\phi_{c}$ |  $<\phi_{c}$

Thus, the choice of the threshold volume fraction would decide whether the
laminate behavior is dominated by the matrix-rich (softer) regions or the
fiber-rich (harder) regions. If $\phi_{c}$ is chosen to be close to $\phi_{\rm
min}$, then most of the laminate would be of type R4 (fiber-rich hard
material), while as we increase $\phi_{c}$, the percentage of softer region
would increase. Finally, when $\phi_{c}$ is close to $\phi_{\rm max}$, then
most of the material would be matrix-rich. Examples of the decomposition for
different choices of $\phi_{c}$ are shown in Fig. 4.

(a)  (b)

(c)

Figure 4: The different regions (R1-R4) when the ratio of the threshold
$\phi_{c}$ to the average $\phi$ (y-fiber) is (a) 0.70, (b) 1.0 and (c) 1.36.
(The color mapping is as in Fig. 3)

We use a threshold $\phi_{c}$, primarily to simplify the heterogeneity in
elastic as well as fracture properties within the composite laminate, such
that role of the hard and soft phases in the growth of damage can be
differentiated easily. The simplification to a two-phase system also brings
out the comparative differences between the behavior of spatially patterned
two phase networks and the randomly distributed two phase networks reported in
the existing literature [32, 33]. In simulation of actual experimental data,
however, a more continuous variation in stiffness and fracture properties
would probably be desirable. Predictions of such a model would have a stronger
dependence on the lattice size.

###  II.2 Random spring network model (RSNM)

The geometric pattern consisting of four regions, obtained from the laminate
using the procedure in Sec. II.1 is now further discretized in space to
convert it into a discrete element model. We discretize a rectangular domain
of the laminate using a square lattice of size $105\times 119$ with lattice
spacing $a=0.5\rm mm$, as shown in Fig. 5a. Each lattice point, representative
of an elemental area, is connected to its nearest and next-nearest neighbors
with linear elastic extensional springs, as shown in Fig. 5b. The resistance
to shear deformation of the spring network is provided by torsional springs
that resist changes to the angle, $\theta_{ijk}$, subtended between the
lattice point and its two adjacent neighbors, as shown in Fig. 5c. For
simulating fracture, a macroscopic crack is modeled by deleting the spring
connections between the lattice points at the top and bottom surface of the
crack which is located at mid-height whose length is approximately $0.28$
times the width of the system.

Figure 5: (a) RSNM with a notch or crack. (b) Magnified view of spring network
showing the spring connections. (c) 2-dimensional square lattice-springs
connection with neighboring and next neighboring lattice points.

When the spring network is deformed from its initial configuration, the total
potential energy stored, $\Phi$, has both extensional component, $\Phi_{\rm
elongation}$ and rotational component, $\Phi_{\rm rotation}$:

$\Phi=\Phi_{\rm elongation}+\Phi_{\rm rotation},$ (1)

For a spring network with $N$ lattice points, the extensional energy is

$\Phi_{\rm
elongation}=\sum_{i=1}^{N}\sum_{j=1}^{4}\frac{1}{2}k_{ij}\\{|\overrightarrow{r_{j}}-\overrightarrow{r_{i}}|-a_{ij}\\}^{2},$
(2)

where $\overrightarrow{r_{i}}$, $\overrightarrow{r_{j}}$ are the position
vectors of lattice sites $i$ and $j$, $a_{ij}$ is the initial equilibrium
distance between them, and $k_{ij}$ is the elastic stiffness of spring joining
$i$ and $j$. The rotational component is given by

$\Phi_{\rm rotation}=\sum_{i=1}^{N}\sum_{\langle
ijk\rangle}\frac{1}{2}c_{ijk}\left(\theta_{ijk}-\frac{\pi}{4}\right)^{2},$ (3)

where the second sum is over the $8$ angles $\theta_{ijk}$ as shown in Fig. 5c
and $c_{ijk}$ is the torsional stiffness of the spring resisting the change of
the angle $\theta_{ijk}$ from $\pi/4$. Using the equivalence of strain energy
density of the spring network and that of a linear elastic isotropic continuum
[37], the elastic constants of the continuum, Young’s modulus $E$ and
Poisson’s ratio $\nu$, are related to the local spring constants as [37]:

$\displaystyle E$ $\displaystyle=$
$\displaystyle\dfrac{8k(k+c/a^{2}_{0})}{3k+c/a^{2}_{0}},$ (4)
$\displaystyle\nu$ $\displaystyle=$
$\displaystyle\dfrac{k-c/a^{2}_{0}}{3k+c/a^{2}_{0}}.$ (5)

In the above RSNM, Poisson’s ratio that can be simulated lie in the range
$(-1,1/3)$.

In the above discussion, the relation between spring constants and the elastic
constants are for a homogeneous material. To extend the RSNM to multi-phase
heterogeneous material, in particular a plain weave laminate, springs
connected to neighbors 1,2,3 and 4 of a lattice point have been assigned as
hard or soft springs based on the region of that location (Fig. 6). The
regions which are rich in $y$-fibers (R3) have hard vertical and hard left
adjacent diagonal springs, while other two springs are soft, and all others
are soft springs. Similarly, regions which rich in $x$-fibers (R2) have hard
horizontal and left adjacent diagonal springs while other springs are soft,
matrix rich regions have all soft springs (R1) and regions rich in $x$\- and
$y$-fibers have all hard springs (R4) respectively.

Figure 6: Different types of springs based on the different regions of the
laminate.

###  II.3 Elastic Modeling

To estimate the elastic behavior of hard and soft springs, we proceed as
follows. The elastic behavior of hard and soft springs depends on the matrix,
the reinforcements, on the choice of the critical threshold $\phi_{c}$, and
the spatial distribution of $\phi$ (x-fiber) and $\phi$ (y-fiber). The choice
of $\phi_{c}$, other than being constrained to be between $\phi_{\rm max}$ and
$\phi_{\rm min}$ of y-fiber, is arbitrary. However, the macroscopic elastic
behavior of the composite cannot depend on the choice of $\phi_{c}$. For this
purpose, elastic moduli of hard regions ($E_{h}$), represented by hard
springs, and soft regions ($E_{s}$), represented by soft springs, of the
laminate at different $\phi_{c}$ are calculated using the rule of mixtures. A
given choice of $\phi_{c}$ decides the fraction $f$ of the laminate which will
be designated as hard. The fraction $f$ of the laminate for a chosen
$\phi_{c}$ is evaluated as

$f=\dfrac{\sum\Delta w_{i}^{h}}{\sum\Delta w_{i}^{h}+\sum\Delta w_{i}^{s}},$
(6)

where $\Delta w_{i}^{h}$ and $\Delta w_{i}^{s}$ denote the hard and soft
phases respectively as defined in Fig. 7. To ensure effective elastic behavior
of the laminate to be independent of the choice of $\phi_{c}$, the elastic
modulus of the laminate, $E_{\rm laminate}$, is then given by the rule of
mixtures to be

$E_{\rm laminate}=E_{h}f+E_{s}(1-f),$ (7)

We assume that a quadratic dependence well captures the dependence of $E_{h}$
and $E_{s}$ on the fraction $f$:

$\displaystyle E_{h}(f)$ $\displaystyle=$ $\displaystyle
a_{0}+a_{1}f+a_{2}f^{2},$ (8) $\displaystyle E_{s}(f)$ $\displaystyle=$
$\displaystyle b_{0}+b_{1}f+b_{2}f^{2}.$ (9)

Figure 7: Schematic diagram showing the discretization of the laminate into
hard (shaded region) and soft phases based on the choice of threshold,
$\phi_{c}$ for $\phi_{\rm min}<\phi_{c}<\phi_{\rm max}$.

The six constants $a_{i}$, $b_{i}$ can be determined as follows. Eq. (7) now
expresses $E_{\rm laminate}$ as a polynomial of degree three in $f$. Since
$E_{\rm laminate}$ cannot depend on the choice of $\phi_{c}$ (and hence $f$),
each of the polynomial coefficients, other than the constant term, must be
zero. This gives $b_{0}=E_{\rm laminate}$, $a_{0}=b_{0}-b_{1}$,
$a_{1}=b_{1}-b_{2}$, $a_{2}=b_{2}$. Thus, only two constants need to be
determined. It is easy to check that the solution is consistent with the
constraint that when $f=0$, i.e. $\phi_{c}=\phi_{\rm max}$, then $E_{s}=E_{\rm
laminate}$. Likewise for the constraint when $f=1$, i.e. $\phi_{c}=\phi_{\rm
min}$, then $E_{h}=E_{\rm laminate}$.

Two more conditions obeyed by the constraint may be obtained by examining the
limits $f\to 0$, and $f\to 1$ more closely. As the volume fraction $\phi$ of
the y-fibers vary in space, the local elastic modulus of the laminate, using
the rule of mixtures, varies between its maximum and minimum values as

$\displaystyle E_{\rm max}$ $\displaystyle=$ $\displaystyle\phi_{\rm
max}E_{\rm fibers}+(1-\phi_{\rm max})E_{\rm matrix},$ (10) $\displaystyle
E_{\rm min}$ $\displaystyle=$ $\displaystyle\phi_{\rm min}E_{\rm
fibers}+(1-\phi_{\rm min})E_{\rm matrix}.$ (11)

When $\phi_{c}\to\phi_{\rm min}$, i.e. $f\to 1$, the material represented by
soft regions corresponds to the laminate part that has fibers at $\phi_{\rm
min}$. This implies that

$E_{s}(f=1)=E_{\rm min}.$ (12)

Applying the same argument to the limit $\phi_{c}\to\phi_{\rm max}$, i.e.
$f\to 0$, we obtain

$E_{h}(f=0)=E_{\rm max}.$ (13)

Eqs. (12) and (13) provide us with the two extra conditions that allow us to
determine the constants $a_{i}$, $b_{i}$. Thus, $E_{h}$ and $E_{s}$ are
determined in terms of $E_{\rm fibers}$, $E_{\rm matrix}$, and $E_{\rm
laminate}$.

We assume that both hard and soft regions have the same Poisson’s ratio $\nu$.
This, in turn, allows us to determine the values of the spring constants $k$
and $c$ for both hard and soft springs using Eqs. (4) and (5). The numerical
values for the elastic constants are taken to be the same as reported for a
typical plain weave laminate (for example, see Ref. [38]). These are tabulated
in Table 2.

Table 2:  Effective material properties for a typical plain weave laminate (taken from Ref. [38]). Properties |  Values
---|---
$E_{\rm matrix}$ (GPa) |  3.12
$E_{\rm fibers}$ (GPa) |  136.5
$E_{\rm laminate}$ (GPa) |  46.7
$\nu$ (laminate) |  0.0654

The solution from the above quadratic interpolation of $E_{h}$ and $E_{s}$ at
different $\phi_{c}$ for representative laminate is plotted in Fig. 8. That
this scheme works well and reproduces the elastic response correctly can be
seen from the macroscopic stress-strain curves shown later [see Fig. 17].

Figure 8: Elastic moduli of hard region, $E_{h}$, and soft region, $E_{s}$, as
a function of the threshold $\phi_{c}$ [see Eqs. (8) and (9)]. $E_{h}$ and
$E_{s}$ satisfy the constraint that $E_{\rm laminate}$ is constant for any
choice of $\phi_{c}$. Inset: The ratio $E_{h}/E_{s}$ has a maximum at
$\phi_{c}=\langle\phi\rangle$ (of y-fiber).  Figure 9: A schematic diagram of
the loading condition. The left bottom corner is pinned while the other
lattice points on the bottom corner can only displace in the horizontal
direction.

###  II.4 Simulation procedure

The lattice points have two translational degrees of freedom. To simulate mode
I fracture, lattice points at the top edge of the domain was given upward
displacement while those at the bottom edge was restricted to translate only
in the x-direction (see Fig. 9). In the simulation, the total macroscopic
strain of 0.1 was applied incrementally in 400 steps. In each step, the top
row of lattice particles were displaced as per the macroscopic strain. The
resulting deformation of the spring network was evaluated iteratively. In the
iterative procedure, the equation of motion

$\overrightarrow{a_{i}}=-\nabla_{\overrightarrow{r_{i}}}\Phi-\gamma\frac{d{\overrightarrow{r_{i}}}}{dt},$
(14)

where $\overrightarrow{a_{i}}$ is the acceleration of lattice site $i$,
$\gamma=0.25$, is the damping coefficient which drives the system to
equilibrium and the mass has been set to $1$, is integrated using the
velocity-verlet algorithm:

$\displaystyle\overrightarrow{r_{i}}(t+\Delta t)$ $\displaystyle=$
$\displaystyle\overrightarrow{r_{i}}(t)(2-\gamma\Delta
t)-\overrightarrow{r_{i}}(t-\Delta t)(1-\gamma\Delta t)$ (15)
$\displaystyle+\overrightarrow{a_{i}}^{\rm bond}(\Delta t)^{2},$

for each incremental time step, $\Delta t\approx 10^{-2}\gamma^{-1}$. We note
that strains are applied quasistatically. The dissipiation term in Eq. (14) is
to ensure relaxation to equilibrium. The time step in these numerical
integration is taken to be much smaller than $1/\gamma$ to ensure time step
independence.

Critical strain and stress approach are used for failure analysis, i.e.,
spring breaks if strain/stress in the respective spring becomes greater than
the critical strain/stress value. If a spring breaks, the system is again
brought into equilibrium until no further spring breakage occurs. To speed up
computation, we implement a parallelized version of the algorithm.

##  III Results

###  III.1 Elastic properties

We first study the effects of the choice of the cut-off threshold, $\phi_{c}$,
and heterogeneity arising from the meso-structure of the laminate on the
elastic stress distribution in the presence of a pre-existing crack. To study
the effects of the choice of $\phi_{c}$, we characterize the spatial stress
distribution in a representative laminate, in which the offsets in the $x$\-
and $y$-directions between layers are randomly chosen, for three choices of
$\phi_{c}$, namely $\phi_{c}/\langle\phi\rangle=0.70,1.0,1.36$. The contours
of the stress $\sigma_{yy}$ in the $y$-direction (when the crack is in the
$x$-direction) are shown in Fig. 10 for a fixed macroscopic strain. These
contours were obtained from the discrete displacement field using two-
dimensional linear interpolation based on the moving least square method, as
discussed in Ref. [39]. In this method, the value of displacement at an
arbitrary location, $P_{i}$ (shown in Fig. 11) is computed from the weighted
least square fit of the displacements of the lattice sites contained in a
circular region $R^{i}$ encircling $P_{i}$. The point $P_{i}$ is then moved
over the entire surface. In Fig. 10a, the stress distribution is shown for the
homogeneous case. This refers to the case when the Young’s modulus is equal to
the effective modulus of the laminate, and is constant throughout the domain
such that the system has discrete translational symmetry. The stress pattern
seen resembles the standard pattern that is observed in a homogeneous system
with a crack, for example, see Ref. [40]. When heterogeneity is introduced
[see Fig. 10b-10d], the stresses no longer have a smooth transition from
highly intensified stresses near the crack tip to far-field lower stresses.
Rather, the stress contours exhibit discontinuous behavior. Also, the highly
stressed zones close to the crack tip, seen in the homogeneous case, become
delocalized when heterogeneity is introduced. When the stress patterns are
compared with the hard and soft bond distributions, as shown in Fig. 4, it is
evident that the hard and soft regions are distinguishable as they develop
significantly different stresses.

(a)  (b)

(c)  (d)

Figure 10: The spatial distribution of the stress $\sigma_{yy}$ for a laminate
whose layers are offset from each other in the $xy$-plane by a random amount.
The macroscopic strain is the same in all the panels which correspond to
different cut-off thresholds $\phi_{c}$. (a) Homogeneous model. (b)
$\phi_{c}=0.7\langle\phi\rangle$. (c) $\phi_{c}=1.0\langle\phi\rangle$. (d)
$\phi_{c}=1.36\langle\phi\rangle$. The data are for a single realization.
Figure 11: A schematic diagram illustrating the moving least square method
approximation to determine displacement of $P_{i}$ based on the weighted least
square fit of the displacements of the lattice sites (shown by crosses) in a
circular region $R^{i}$ encircling $P_{i}$. Another example is shown for point
$P_{i+1}$.

The sharpest contrast between neighboring regions nearer to crack tip is
observed when $\phi_{c}$ is chosen to be $\langle\phi\rangle$. Increasing or
decreasing $\phi_{c}$ away from the mean reduces the contrast and the contours
are closer to the homogeneous continuous patterns. These differences are
quantified in terms of the stress normal to the plane of the initial crack
ahead of the crack tip, as seen in Fig. 12. Ahead of the crack tip, in the
homogeneous solution, the opening stress decays rapidly in a smooth manner. In
the heterogeneous combinations, the maximum deviations from the homogeneous
solution are seen when $\phi_{c}=\langle\phi\rangle$, which is a consequence
of $E_{h}/E_{s}$ being the highest for a given $E_{\rm eff}$. As the threshold
is increased, even though the hard phase develops significantly higher
stresses, there are fewer locations that belong to the hard phase and the
stresses in the softer phase are, thus, closer to the homogeneous solution.
The same effect is also seen for low threshold as very few locations now
belong to the softer phase. In the absence of any compelling reason to choose
otherwise, as well as to enhance the effects of heterogeneity, we will
consider $\phi_{c}=\langle\phi\rangle$ in the remainder of the paper.

Figure 12: The variation of the stress $\sigma_{yy}$ with distance ahead of
the crack tip ($r$) for different choices of the cut-off threshold. The data
are for a laminate whose offsets between layers are randomly chosen. The
results for the homogeneous model are shown for comparison.

For the same effective elastic modulus, the relative offset between the layers
can also have a strong effect on the spatial distribution of the volume
fractions of $x$\- and $y$-fibers as well as its maximum and minimum limits.
To establish the role of relative offset in comparison to the homogeneous
material, we consider two extreme configurations which we call as zero offset
and staggered offset configurations. In the zero offset configuration, the
fibers are perfectly aligned across the thickness, as shown in Fig. 13a, such
that, along the length of the laminate, the fiber volume fraction varies from
maxima to zero (similar to a lamina). In the staggered offset configuration,
each adjacent layer has an offset of exactly half the wavelength of the
undulation in both $x$\- and $y$-directions, as shown in Fig. 13b, resulting
in fiber volume fraction having double the frequency of the undulations.
Contours of the opening stress field $\sigma_{yy}$ are presented in Fig. 14
for different types of laminate configurations in comparison to that for a
homogeneous material. In the zero-offset configuration shown in Fig. 14b,
meso-structure patterns dominate the stress distribution, accompanied by
significant delocalization of the high stresses near the crack tip. However,
in the staggered configuration, in Fig. 14c, even though the stress
distribution has patterns of the heterogeneity, contour shapes are much closer
to the homogeneous model. Higher disorder in the elastic field, as seen in the
case of the zero-offset configuration contributes to its observed increase in
strength and toughness, as we will see later.

(a)  (b)

Figure 13: Cross-sectional view of the two extreme configurations of the
laminate. The panels correspond to meso-structures with (a) zero offset
between layers and (b) staggered offset between layers.

(a)  (b)

(c)

Figure 14: The spatial distribution of the stress $\sigma_{yy}$ for fixed cut-
off $\phi_{c}=\langle\phi\rangle$ and different laminate configurations. The
macroscopic strain is the same in all the panels which correspond to (a)
homogeneous model (b) zero offset meso-structure and(c) staggered offset meso-
structure. The data are for a single realization.

The distribution of $\sigma_{yy}$, for the different configurations
considered, is presented in Fig. 15. The single peak observed in the
homogeneous case corresponds to the lattice points in the region away from the
crack tip. On the introduction of heterogeneity, two peaks form as the strains
are compatible in the neighborhood but two different material behavior result
in differences in the $\sigma_{yy}$. In the zero offset configuration, the
ratio between the elastic moduli of the two phases is largest, and thus the
peaks are furthest away from each other, leading to higher stresses for a
larger number of lattice points, whereas, in staggered offset, peaks are much
closer. For random-offset configuration, we find that the peaks are in between
(not plotted for the sake of clarity).

Figure 15: Distribution of stress $\sigma_{yy}$ of lattice points for
homogeneous system and system when heterogeneity is included. All the data
correspond to the same macroscopic strain for single realization.

###  III.2 Fracture properties

We now study the effect of elastic heterogeneity, inherent in the laminates,
on the fracture behavior. In the existing studies of fracture behavior of
heterogeneous materials using network models, two main simplifications have
been adopted in relation to the failure criteria for inter-particle
interaction: a common threshold failure strain [34] or a common failure
threshold stress [32], as shown in Fig. 16. Since both elastic and failure
properties of the springs affect the macroscopic response significantly, in
the present work, we employ both the criteria to establish a comprehensive
understanding of the role of spatial patterns in heterogeneity on its fracture
process. In a common stress threshold criteria, we choose failure criteria for
the bonds to be such that for all springs, the stress threshold, $\sigma^{*}$,
is taken to be the same, as shown in Fig. 16a. This results in softer bonds
requiring the highest energy for failure while harder bonds require the least,
and the bonds of the effective homogeneous medium require energy in between.
For a common strain threshold criteria, as the failure strain threshold is
taken to be the same ($=\epsilon^{*}$) for all phases, softer bonds require,
thus, the least energy for failure compared to the harder bonds, as evident in
Fig. 16b. We do simulations for laminates with zero and staggered offsets, and
for comparison, we also simulate fracture response of an effective homogeneous
laminate.

(a)  (b)

Figure 16: The stress-strain diagram for springs depicting their failure
criteria based on (a) common stress threshold and (b) common strain threshold.

In addition to the elastic heterogeneity, we also account for the inherent
disorder in the fracture properties, arising from microscopic variations and
defects, by assigning the failure threshold strain for each spring
independently from a Gaussian distribution with mean at $\sigma^{*}/E$ or
$\epsilon^{*}$ as shown in Fig. 16 and a standard deviation taken to be within
the range of zero to 10% in all the simulations of fracture. To show that the
comparative responses of the three laminate meso-structures differ from each
other, depending on the type of failure criteria chosen, we take an initial
crack of length $a_{0}=0.28w$ and perform 25 realizations each for 3 different
meso-structures: effective homogeneous, zero offset and staggered offset. We
assume the disorder in the failure strain threshold to have a standard
deviation of $5\%$ of the mean. We choose the failure strain threshold for the
effective homogeneous laminate arbitrarily to be $\epsilon^{*}=0.025$. From
the responses based on a common failure stress threshold shown in Fig. 17a,
the staggered configuation has a response very similar to the homogeneous
case, while the configuration with zero offset exhibits significantly enhanced
strength. The area under the response curve, a measure of the material’s
toughness, is also highest for the zero offset configuration.

(a)  (b)

Figure 17: The force-displacement curve of laminates with different
configuration for a crack of length $a_{0}=0.28w$. The disorder in the failure
threshold corresponds to standard deviation $5\%$ of the mean. The panels
correspond to data for (a) common stress threshold and (b) common strain
threshold.

Even though the meso-structure with zero-offset exhibits toughening in both
the common stress and strain threshold criteria, there are differences in
dissipative mechanisms when considering the relative contribution of hard
bonds and soft bonds towards the total energy that is dissipated during crack
growth. The number of extensional springs failing during crack growth is
primarily driven by the energy costs involved. In the common stress threshold
criterion, since the soft bonds are tougher, fewer soft bonds break compared
to hard bonds. This is true for both zero-shift as well as staggered
configuration, as is evident from Fig. 18a and 18b. Similarly, the common
strain threshold criterion implies the harder bonds are tougher resulting in
fewer hard bonds breaking compared to soft bonds, as can be seen for both
meso-structures in Fig. 18c and 18d, respectively. However, for both the
criteria, the total number of failed springs are higher for the zero-shift
configuration than the effective homogeneous counterpart, as seen in Fig. 18a
and 18c, resulting in its enhanced toughness.

(a)  (b)

(c)  (d)

Figure 18: The cumulative number of broken or failed bonds with increasing
displacement. The panels (a) and (c) are for the zero-offset meso-structure
while the panels (b) and (d) are for the staggered offset meso-structure. The
panels (a) and (b) correspond to a common stress threshold failure criteria
while the panels (c) and (d) correspond to a common constant strain threshold
criteria.

The spatial location of the broken bonds is shown in Fig. 19 for
representative realizations of zero-offset and staggered offset meso-
structures for both common stress as well as strain threshold failure
criteria. The broken hard, soft, and unbroken bonds are shown in orange
[grey], black [black] and light green [light grey] colors respectively. The
crack path is most tortuous for zero-offset configuration with a common stress
failure threshold, as can be seen in Fig. 19a. This is because the crack
preferably avoids breaking soft bonds as they cost larger energy. It is also
observed that a few hard bonds that are not directly connected to the crack
path also fail. Thus, the damage is spread over a larger area, presumably
resulting in higher toughness. The crack path for zero-offset configuration
for a common strain failure criterion has a comparatively less tortuous path,
though as expected, there is a clear preference for breakage of softer bonds
(see Fig. 19c). Also, all the broken bonds are connected to the main crack. On
the other hand, the staggered configuration, for both choices of failure
criteria, exhibits a very localized path of the propagating crack along the
initial crack plane, as shown in Figs. 19b and 19d. This is because
$E_{h}/E_{s}$ is closer to unity and the differentiation between soft and hard
bonds is minimal.

(a)  (b)

(c)  (d)

Figure 19: The final geometry of the crack for different meso-structures and
failure criteria. The orange [grey], black [black] and light green [light
grey] colors represent broken hard, soft, and unbroken bonds respectively. The
panels (a) and (c) are for the zero-offset meso-structure while the panels (b)
and (d) are for the staggered offset meso-structure. The panels (a) and (b)
correspond to a common stress threshold failure criterion while the panels (c)
and (d) correspond to a common constant strain threshold criterion.

To gain an insight into the role of disorder in failure threshold on fracture
for different laminate configurations with their inherent elastic
heterogeneity, we simulate the fracture response for a range of disorder by
varying the standard deviation of the Gaussian distribution for failure
threshold from $0\%$ \- $10\%$. When the hard and soft phase have a common
stress threshold, the zero-offset configuration develops the highest strength
as well as toughness for the entire range of disorder considered, as seen in
Fig. 20a-20d. With slightest increase in disorder, strength of the elastically
homogeneous solid drops sharply and is nearly a constant thereafter. The
effect of increasing disorder has a marginal effect on the strength of
staggered configuration and, except for very low disorder, the failure
threshold the strength is very close to the effective elastically homogeneous
laminate.

(a)  (b)

(c)  (d)

Figure 20: Effect of disorder in failure criteria on the peak force and
fracture toughness of laminate. The different panels correspond to (a), (b):
common stress threshold criterion and (c), (d): common strain threshold
criterion. The data are for a crack of length $a_{0}=0.28w$, and have been
averaged over 25 realizations (50 for homogeneous).

We now study the effect of elastic heterogeneity and the failure threshold
criteria on the critical stress intensity factor or fracture toughness
$K_{Ic}$. As per the linear elastic fracture theory, the failure stress
$\sigma_{\rm max}$ depends on the crack length $a_{0}$ as

$\frac{1}{\sigma^{2}_{\rm max}}=\frac{\beta^{2}a_{0}}{K^{2}_{Ic}},$ (16)

where the geometric correction factor $\beta$ for a given height to width
ratio of the specimen geometry is known [41]. In presence of a fracture
process zone, it has been proposed that Eq. (16) is modified to

$\frac{1}{\sigma^{2}_{\rm max}}=\frac{\beta^{2}(a_{0}+\xi)}{K^{2}_{Ic}},$ (17)

where $\xi$ is another additional length length scale [42, 43, 44]. In
presence of strong disorder, it has been argued that the crack length $a_{0}$
in Eq. (17) is dependent also on the disorder [45, 46]. In our simulation, the
disorder is small, with damage limited to near crack, and we will therefore
compare our data with Eq. (17). To confirm the linear relation between
$1/\sigma^{2}_{\rm max}$ and $a_{0}$, as well as to measure $K_{Ic}$, we
simulate the fracture response at $5\%$ standard deviation in failure strength
for a range of initial crack lengths, $a_{0}$, as shown in Fig. 21. It is
evident that, in our simulations, $1/\sigma^{2}_{\rm max}$ is linearly
proportional to $a_{0}$ for different elastic heterogeneities as well as
different failure threshold criteria. Interestingly, for the constant stress
failure threshold [see Fig. 21a], we find that $K_{Ic}$ is larger for zero-
offset laminate than the homogeneous or staggered laminate, thus increasing
fracture toughness. For a common strain threshold failure criterion [see Fig.
21b], we find that on an average $K_{Ic}$ is independent of the elastic
heterogeneity. However, on a more local scale, there are deviations from
linearity for the zero-offset laminate, and $K_{Ic}$ depends on the initial
location of the tip of the crack front.

(a)  (b)

Figure 21: The variation of the failure stress $\sigma_{\rm max}$ with crack
length $a_{0}$. The results are for laminates with staggered offset and zero
offset, as well as the homogeneous model, with $5\%$ standard deviation in the
failure threshold. The panels correspond to failure thresholds based on (a)
common stress criterion and (b) common strain criterion. Each data point has
been averaged over 25 realizations (50 for homogeneous).

(a)  (b)

Figure 22: The variation of the avalanche size distribution, $P(s)$, with
avalanche size, $s$, for different laminate configurations, as well as for the
homogeneous model. The panels represent data for (a) a common stress failure
thresholds, and (b) a common strain failure thresholds. The data are for
$10\%$ standard deviation in failure threshold and averaged over 250
realizations.

Another commonly used probe to characterize the fracture process is statistics
of acoustic emissions during fracture, or equivalently avalanche statistics in
simulations. An avalanche is defined as the number of bonds that are broken in
one increment of applied strain. Let $P(s)$ denote the probability density
function for avalanches of size $s$. The distribution is usually insensitive
to the nature of the material undergoing fracture and follows a power law
distribution $P(s)\sim s^{-\alpha}$ for large $s$. The avalanche size
distribution, for $10\%$ standard deviation in failure strength, is shown in
Fig. 22 for the different laminate configurations for both choices of the
critical thresholds and compared to that for the homogeneous model. For
intermediate $s$, we see that the data for the different meso-structure follow
power law distributions for both common stress theshold (see Fig. 22a) as well
as common strain threshold (see Fig. 22b). In this range of s, there is no
distinguishable effect of meso-structure on the exponent $\alpha$. However,
the range of the power law is limited, and a more definite statement on the
universality requires simulation of much larger system sizes. We conclude that
exponent $\alpha$ is atmost weakly dependent on the meso-structure as well as
the choice of threshold. Note that the power law distribution is quite
different from the avalanche distributions obtained from RSNM simulations, in
the absence of a pre-existing notch or crack, which have exponents close to
2.5 [22]. However, close to breakdown, it is known from exact solution of the
fiber bundle model as well as simulations of the random fuse model that the
exponent for critical avalanches could be very different from the exponents of
avalanche distribution away from breakdown [25]. For instance, in the fiber
bundle model, the exponent changes from $5/2$ to $3/2$ near breakdown, while
in a network of electric fuses the exponent decreases from 3 to 2 near the
catastrophic failure. The corresponding exponent for critical avalanches in
RSNM, while not known, is expected to decrease from 2.5 [47, 25]. In our
simulations, where there is a pre-existing crack in the system as well as when
the standard deviation of the failure threshold distributions is very small,
the system can be thought of as one close to breakdown and the avalanche
distribution is presumably closer to near failure condition, and thus in
analogy with other models it is expected to have an exponent significantly
lower than 2.5 (between 1.1 and 1.3), as shown in Fig. 22.

We note that the avalanche distribution is not stationary, i.e., it is not
independent of the value of the strain. This is clear since there are no
avalanches for small strains. For the choice of disorder, the integrated
$P(s)$ is dominated by avalanche near breakdown.

##  IV Conclusion

Fracture processes in composite materials strongly depend not only on the
properties of the constituents but also equally on their architecture. In this
paper, we develop a methodology for constructing a discrete element model for
composites, particularly a plain weave laminate, which incorporates the
inherent spatial patterns of the meso-structure. Though the macroscopic
elastic behavior does not depend on the spatial patterning of its
constituents, the spatial distribution of the stresses for a composite is
distinctly different from that of an elastically equivalent homogeneous system
as the high stresses near the crack tip are delocalized significantly.

For studying fracture behavior, two failure threshold criteria were used: a
common stress threshold and a common strain threshold. Interestingly, the
enhancement in strength and fracture toughness is observed for all initial
crack lengths when the failure criteria have a common stress threshold. The
enhancement is most pronounced for the composite laminate with zero offset.
The dominant mechanism is observed to be the breakage of bonds that are
energetically less dissipative. Consequently, when energetically more
dissipative bonds are encountered in the path of a growing crack, it tends to
deflect from its plane, resulting in a more tortuous path and a larger number
of total bonds broken. In contrast, for a common strain threshold, while there
is no enhancement in toughness on an average, occasional enhancement in
toughness is observed depending on the initial position of the crack tip.

We also studied the effect of spatial patterning on the avalanche size
distribution. Within the simulation data, the avalanche distribution is
independent of the elastic heterogeneity as well as the failure criteria.
However, the value of the avalanche exponent is much smaller than the value
reported for RSNM. In the presence of a crack as well as very little disorder,
we argued that the avalanche resemble critical avalanche close to breakdown.
In such case, it is known for other models like fiber bundle models and random
fuse model that the critical avalanche exponent could be significantly lower
than off-critical exponents. Determining this exponent more accurately for
RSNM is a promising area for future study.

The simulations provide valuable insights into the role of elastic
heterogeneity in a two-phase network in the enhancement of toughness for the
two specific scenarios of comparative failure thresholds. In particular, a
common stress threshold criterion in which the harder phase is weaker and the
softer phase is tougher is a common comparative material failure response.
However, for the realistic approximation of actual experimental data, the
imposition of a constraint such as a common failure threshold would be
inapplicable. Equivalent statistical signatures of two-phase behavior, as
observed in the count of broken bonds, if found in experimental data such as
acoustic emission activity, can provide vital inputs in establishing a failure
criterion for the two phases.

## References

  * Gao  _et al._ [2003] H. Gao, B. Ji, I. L. Jäger, E. Arzt, and P. Fratzl, Proceedings of the National Academy of Sciences 100, 5597 (2003).
  * Gupta  _et al._ [2009] H. Gupta, S. Krauss, J. Seto, W. Wagermaier, M. Kerschnitzki, G. Benecke, P. Zaslansky, P. Boesecke, S. Funari, H. Kirchner, and P. Fratzl, Bone 44, S33 (2009).
  * Sen and Buehler [2011] D. Sen and M. J. Buehler, Scientific Reports 1, 35 (2011).
  * Piggott [1994] M. R. Piggott, Journal of Composite Materials 28, 588 (1994).
  * Puck and Schürmann [2002] A. Puck and H. Schürmann, Composites Science and Technology 62, 1633 (2002).
  * Mishnaevsky and Brøndsted [2009] L. Mishnaevsky and P. Brøndsted, Computational Materials Science 44, 1351 (2009).
  * Daniel  _et al._ [2009] I. M. Daniel, J.-J. Luo, P. M. Schubel, and B. T. Werner, Composites Science and Technology 69, 764 (2009).
  * Wicks  _et al._ [2014] S. S. Wicks, W. Wang, M. R. Williams, and B. L. Wardle, Composites Science and Technology 100, 128 (2014).
  * Das  _et al._ [2018] S. Das, K. Kandan, S. Kazemahvazi, H. Wadley, and V. Deshpande, International Journal of Solids and Structures 136-137, 137 (2018).
  * Tan  _et al._ [2000] P. Tan, L. Tong, G. Steven, and T. Ishikawa, Composites Part A: Applied Science and Manufacturing 31, 259 (2000).
  * Carol  _et al._ [2001a] I. Carol, E. Rizzi, and K. Willam, International Journal of Solids and Structures 38, 491 (2001a).
  * Carol  _et al._ [2001b] I. Carol, E. Rizzi, and K. Willam, International Journal of Solids and Structures 38, 519 (2001b).
  * Luccioni and Oller [2003] B. Luccioni and S. Oller, Computer Methods in Applied Mechanics and Engineering 192, 1119 (2003).
  * Ansar  _et al._ [2011] M. Ansar, W. Xinwei, and Z. Chouwei, Composite Structures 93, 1947 (2011).
  * Zhang  _et al._ [2012] T. Zhang, X. Li, S. Kadkhodaei, and H. Gao, Nano Letters 12, 4605 (2012).
  * Camanho  _et al._ [2013] P. Camanho, M. Bessa, G. Catalanotti, M. Vogler, and R. Rolfes, Mechanics of Materials 59, 36 (2013).
  * Vogler  _et al._ [2013] M. Vogler, R. Rolfes, and P. Camanho, Mechanics of Materials 59, 50 (2013).
  * Tan and Falzon [2016] W. Tan and B. G. Falzon, Composites Science and Technology 126, 60 (2016).
  * Rosti  _et al._ [2009] J. Rosti, X. Illa, J. Koivisto, and M. J. Alava, Journal of Physics D: Applied Physics 42, 214013 (2009).
  * Baró  _et al._ [2013] J. Baró, A. Corral, X. Illa, A. Planes, E. K. H. Salje, W. Schranz, D. E. Soto-Parra, and E. Vives, Phys. Rev. Lett. 110, 088702 (2013).
  * Curtin and Scher [1990] W. A. Curtin and H. Scher, Journal of Materials Research 5, 535 (1990).
  * Alava  _et al._ [2006] M. J. Alava, P. K. Nukala, and S. Zapperi, Advances in Physics 55, 349 (2006).
  * Pan  _et al._ [2018] Z. Pan, R. Ma, D. Wang, and A. Chen, Engineering Fracture Mechanics 190, 382 (2018).
  * Ray [2006] P. Ray, Computational Materials Science 37, 141 (2006).
  * Pradhan  _et al._ [2006] S. Pradhan, A. Hansen, and P. C. Hemmer, Phys. Rev. E 74, 016122 (2006).
  * Karihaloo  _et al._ [2003] B. Karihaloo, P. Shao, and Q. Xiao, Engineering Fracture Mechanics 70, 2385 (2003).
  * Berton and Bolander [2006] S. Berton and J. E. Bolander, Computer Methods in Applied Mechanics and Engineering 195, 7172 (2006).
  * Grassl and Bažant [2009] P. Grassl and Z. P. Bažant, Journal of Engineering Mechanics 135, 85 (2009).
  * Mayya  _et al._ [2016] A. Mayya, P. Praveen, A. Banerjee, and R. Rajesh, Journal of the Royal Society Interface 13, 20160809 (2016).
  * Mayya  _et al._ [2017] A. Mayya, A. Banerjee, and R. Rajesh, Phys. Rev. E 96, 053001 (2017).
  * Mayya  _et al._ [2018] A. Mayya, A. Banerjee, and R. Rajesh, Journal of the Mechanical Behavior of Biomedical Materials 83, 108 (2018).
  * Urabe and Takesue [2010] C. Urabe and S. Takesue, Phys. Rev. E 82, 016106 (2010).
  * Boyina  _et al._ [2015] D. Boyina, T. Kirubakaran, A. Banerjee, and R. Velmurugan, Mechanics of Materials 91, 64 (2015).
  * Dimas  _et al._ [2014] L. S. Dimas, T. Giesa, and M. J. Buehler, Journal of the Mechanics and Physics of Solids 63, 481 (2014).
  * Dimas  _et al._ [2015a] L. S. Dimas, D. Veneziano, T. Giesa, and M. J. Buehler, Journal of Applied Mechanics 82, 011003 (2015a).
  * Dimas  _et al._ [2015b] L. S. Dimas, D. Veneziano, T. Giesa, and M. J. Buehler, Journal of the Mechanics and Physics of Solids 84, 116 (2015b).
  * Monette and Anderson [1994] L. Monette and M. P. Anderson, Modelling and Simulation in Materials Science and Engineering 2, 53 (1994).
  * Boyina  _et al._ [2014] D. Boyina, A. Banerjee, and R. Velmurugan, Composites Part B: Engineering 60, 21 (2014).
  * Liu and Long [2015] W. Liu and R. Long, Journal of Applied Mechanics 83, 011006 (2015).
  * Chona  _et al._ [1982] R. Chona, G. R. Irwin, and A. Shukla, The Journal of Strain Analysis for Engineering Design 17, 79 (1982).
  * Hammond and Fawaz [2016] M. J. Hammond and S. A. Fawaz, Engineering Fracture Mechanics 153, 25 (2016).
  * Bažant [1984] Z. P. Bažant, Journal of Engineering Mechanics 110, 518 (1984).
  * Bažant [1996] Z. P. Bažant, Advanced Cement Based Materials 4, 128 (1996).
  * Bažant [2004] Z. P. Bažant, Proceedings of the National Academy of Sciences 101, 13400 (2004).
  * Alava  _et al._ [2008] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi, Phys. Rev. Lett. 100, 055502 (2008).
  * Papanikolaou  _et al._ [2019] S. Papanikolaou, P. Shanthraj, J. Thibault, C. Woodward, and F. Roters, Materials Theory 3, 5 (2019).
  * Pradhan  _et al._ [2005] S. Pradhan, A. Hansen, and P. C. Hemmer, Phys. Rev. Lett. 95, 125501 (2005).


# Better Approximate Inference for Partial Likelihood Models with a Latent
Structure

Amrith Setlur
School of Computer Science
Carnegie Mellon University
Pittsburgh, PA 15213
<EMAIL_ADDRESS>
&Barnabás Póczós
School of Computer Science
Carnegie Mellon University
Pittsburgh, PA 15213
<EMAIL_ADDRESS>

###### Abstract

Temporal Point Processes (TPP) with partial likelihoods involving a latent
structure often entail an intractable marginalization, thus making inference
hard. We propose a novel approach to Maximum Likelihood Estimation (MLE)
involving approximate inference over the latent variables by minimizing a
tight upper bound on the approximation gap. Given a discrete latent variable
$Z$, the proposed approximation reduces inference complexity from
$\mathcal{O}(|Z|^{c})$ to $\mathcal{O}(|Z|)$. We use convex conjugates to
determine this upper bound in a closed form and show that its addition to the
optimization objective results in improved results for models assuming
proportional hazards as in Survival Analysis.

##  1 Introduction

Temporal Point Processes (TPPs) provide a formal framework to model the
occurrences of discrete events in time (like failures or financial
transactions). Recent work ((Linderman and Adams, 2014) (Snoek et al., 2013))
on modelling TPPs with latent factors have showcased their ability to capture
correlations such as inhibitory relationships & a dichotomy of classes of
neurons in neural spike recordings. Although, there have been several advances
in non-parametric Bayesian inference ((Samo and Roberts, 2015)) most models
are parametric ((Cox, 1955)) where parameter estimation is done by maximizing
the likelihood of observed point values. Survival analysis is the problem of
estimating survival times for entities (like nodes in a machine) and it has
largely relied on TPPs to estimate survival times in the presence of censored
observations. Semi-parametric methods like the Cox Proportional Hazards (CPH)
(Cox, 1955) allow parametric estimation using a partial likelihood objective
without estimating the baseline hazard. Therefore, we propose an approximate
inference strategy for latent variable models with a partial likelihood
objective. We introduce an inference method for models where the normalization
factor includes interactions over log-linear factors. Such models are common
in TPPs assuming proportional hazards ((Rosen and Tanner, 1999)) or in latent
Conditional Random Fields (CRFs) where the normalization involves a sum over
finite potential functions ((Sutton et al., 2012)). (Rosen and Tanner, 1999)
introduce an inference strategy for CPH which is similar to our proposed
method, but they fail to identify cases where the approximation fails.
Although our inference strategy is applicable to the full likelihood in a TPP,
we focus on its impact in the case of partial likelihoods since the objective
there is closely related to the MLE objectives observed in latent CRFs, thus
making our work applicable to a broader class of problems.

Inspired by (Jebara and Choromanska, 2012) we introduce a distribution
agnostic closed form tight upper bound on likelihood estimations for TPPs
resembling (Diggle, 2005). The upper bound can be minimized via standard
gradient descent based iterative methods (Ruder, 2016). Finally, we prove a
tight upper bound on the Jensen inequality for strictly convex polynomial
functions on $\mathbb{R}_{++}$.

##  2 The Inference Problem

Given a compact set $S$ equipped with Borel $\sigma$-algebra $\mathbb{B}(S)$,
$X(t):\Omega\to\Gamma$ is a TPP if $X(t)$ is a measurable transform from the
probability space $(\Omega,\mathfrak{F},\mathbb{P})$ to the space of counting
measures $\Gamma$ on S. Given a series of events
$\\{(\delta_{i},x_{i},t_{i})\\}_{N}$: for a given entity $x_{i}$ the event of
interest occurred $(\delta_{i}=1)$ at time $t_{i}$ or the observation was
censored $(\delta_{i}=0)$. The risk set for $x_{i}$ is given by
$R(t_{i})=\\{x_{j}:t_{j}\geq t_{i}\\}$. Under a Poisson TPP with intensity
$\lambda(x)$, the likelihood of the event $\delta_{i}=1$ for $x_{i}$ at
$t_{i}$ given that the event hasn’t occurred till $t_{i}$ is given by
$\mathbb{P}(N(B_{a})=1)/\mathbb{P}(N(B_{b})=0)$.
$B_{a},B_{b}\in\mathbb{B},B_{a}=[t_{i},t_{i}+\delta),B_{b}=[0,t_{i})$ where
$N(B)$ follows a Poisson distribution with mean $\int_{B}\lambda(x)dx$.

We modify the formulation by adding latent variables $z$ (see figure 1) and
now the intensity function in the TPP is function of parameters $\beta$, input
$x_{i}$ and latent variable $z_{i}\sim p(.|x_{i},\theta)$. (Rosen and Tanner,
1999) and (Diggle, 2005) used partial likelihood models to efficiently compute
the MLE estimates for the parameters in an inhomogenous Poisson process.
Partial likelihood was first introduced by (Cox, 1955) with the aim of
identifying variables that impact survival analysis without worrying about the
baseline hazard. For the same reasons, we choose to maximize the partial
likelihood of an event $(\delta_{i},x_{i},t_{i})$ conditioned on the risk set
$R(t_{i})$. Thus the denominator in eq. 1 now involves a sum over a finite set
of factors from $R(t_{i})$ (closely resembling latent CRFs (Quattoni et al.,
2007) in likelihood estimation).

$\displaystyle\vbox{\hbox{ \leavevmode\hbox to105.29pt{\vbox
to71.51pt{\pgfpicture\makeatletter\hbox{\hskip
51.31305pt\lower-58.56853pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{ }
\pgfsys@beginscope\pgfsys@invoke{
}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{}\pgfsys@moveto{7.90042pt}{0.0pt}\pgfsys@curveto{7.90042pt}{4.36333pt}{4.36333pt}{7.90042pt}{0.0pt}{7.90042pt}\pgfsys@curveto{-4.36333pt}{7.90042pt}{-7.90042pt}{4.36333pt}{-7.90042pt}{0.0pt}\pgfsys@curveto{-7.90042pt}{-4.36333pt}{-4.36333pt}{-7.90042pt}{0.0pt}{-7.90042pt}\pgfsys@curveto{4.36333pt}{-7.90042pt}{7.90042pt}{-4.36333pt}{7.90042pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{-6.02776pt}{-1.74998pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\lambda(t)$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}{}{}{{{{}}}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{0.875,0.875,0.875}\pgfsys@color@gray@fill{0.875}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\definecolor{pgffillcolor}{rgb}{0.875,0.875,0.875}\pgfsys@color@gray@fill{0.875}\pgfsys@invoke{
}{}\pgfsys@moveto{50.91307pt}{0.0pt}\pgfsys@curveto{50.91307pt}{3.86601pt}{47.77911pt}{6.99997pt}{43.9131pt}{6.99997pt}\pgfsys@curveto{40.04709pt}{6.99997pt}{36.91313pt}{3.86601pt}{36.91313pt}{0.0pt}\pgfsys@curveto{36.91313pt}{-3.86601pt}{40.04709pt}{-6.99997pt}{43.9131pt}{-6.99997pt}\pgfsys@curveto{47.77911pt}{-6.99997pt}{50.91307pt}{-3.86601pt}{50.91307pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{43.9131pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{41.91277pt}{-1.50693pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$x$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}{}{}{{{{}}}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{}\pgfsys@moveto{6.99997pt}{-43.9131pt}\pgfsys@curveto{6.99997pt}{-40.04709pt}{3.86601pt}{-36.91313pt}{0.0pt}{-36.91313pt}\pgfsys@curveto{-3.86601pt}{-36.91313pt}{-6.99997pt}{-40.04709pt}{-6.99997pt}{-43.9131pt}\pgfsys@curveto{-6.99997pt}{-47.77911pt}{-3.86601pt}{-50.91307pt}{0.0pt}{-50.91307pt}\pgfsys@curveto{3.86601pt}{-50.91307pt}{6.99997pt}{-47.77911pt}{6.99997pt}{-43.9131pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{-43.9131pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{-1.78159pt}{-45.42003pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$z$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}{}{}{{{{}}}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{}\pgfsys@moveto{-36.91312pt}{0.0pt}\pgfsys@curveto{-36.91312pt}{3.86601pt}{-40.04707pt}{6.99997pt}{-43.91309pt}{6.99997pt}\pgfsys@curveto{-47.7791pt}{6.99997pt}{-50.91306pt}{3.86601pt}{-50.91306pt}{0.0pt}\pgfsys@curveto{-50.91306pt}{-3.86601pt}{-47.7791pt}{-6.99997pt}{-43.91309pt}{-6.99997pt}\pgfsys@curveto{-40.04707pt}{-6.99997pt}{-36.91312pt}{-3.86601pt}{-36.91312pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{-43.91309pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{-45.89276pt}{-1.74998pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\beta$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}{}{}{{{{}}}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,1,1}\pgfsys@color@gray@fill{1}\pgfsys@invoke{
}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}{}\pgfsys@moveto{-36.01266pt}{-43.91296pt}\pgfsys@curveto{-36.01266pt}{-40.04695pt}{-39.14662pt}{-36.913pt}{-43.01263pt}{-36.913pt}\pgfsys@curveto{-46.87865pt}{-36.913pt}{-50.0126pt}{-40.04695pt}{-50.0126pt}{-43.91296pt}\pgfsys@curveto{-50.0126pt}{-47.77898pt}{-46.87865pt}{-50.91293pt}{-43.01263pt}{-50.91293pt}\pgfsys@curveto{-39.14662pt}{-50.91293pt}{-36.01266pt}{-47.77898pt}{-36.01266pt}{-43.91296pt}\pgfsys@closepath\pgfsys@moveto{-43.01263pt}{-43.91296pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{-44.65567pt}{-46.34349pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\theta$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}}{{{{}}}}{{{{}}}}{{{{}}}}{{{{}}}}{{}}{{{{}}}}{{{{}}}}{{{{}}}}{{{{}}}}{{}}{{{{}}}}{{{{}}}}{{{{}}}}{{{{}}}}{}
{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{0.72563pt}{-21.50618pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{\minipage[t]{59.37328pt}\endminipage}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}}{{}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{}
{{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{37.17664pt}{-51.67322pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{\footnotesize{$N$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{ {}{}{}}{{{
{}{}{}}}}{{{{}{}{}}}}{{{ {}{}{}}}}{{{ {}{}{}}}}{ {}{}{}}{{{
{}{}{}}}}{{{{}{}{}}}}{{{ {}{}{}}}}{{{ {}{}{}}}}{}
{{}}{{}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{
{}{}} {{}{{}}}{{}{}}{}{{}{}}{}{}{}{}{}
{{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{{}{}{{}}}{{}{}{{}}}{}{}{}\pgfsys@moveto{49.57594pt}{12.5418pt}\pgfsys@lineto{-6.56389pt}{12.5418pt}\pgfsys@curveto{-8.77306pt}{12.5418pt}{-10.56389pt}{10.75096pt}{-10.56389pt}{8.5418pt}\pgfsys@lineto{-10.56389pt}{-54.16853pt}\pgfsys@curveto{-10.56389pt}{-56.3777pt}{-8.77306pt}{-58.16853pt}{-6.56389pt}{-58.16853pt}\pgfsys@lineto{49.57594pt}{-58.16853pt}\pgfsys@curveto{51.78511pt}{-58.16853pt}{53.57594pt}{-56.3777pt}{53.57594pt}{-54.16853pt}\pgfsys@lineto{53.57594pt}{8.5418pt}\pgfsys@curveto{53.57594pt}{10.75096pt}{51.78511pt}{12.5418pt}{49.57594pt}{12.5418pt}\pgfsys@closepath\pgfsys@moveto{-10.56389pt}{-58.16853pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7}{0.0}{0.0}{0.7}{6.76376pt}{-22.81337pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{\minipage[t]{42.12096pt}\endminipage}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} ;
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope {{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{{
{\pgfsys@beginscope\pgfsys@setdash{}{0.0pt}\pgfsys@miterjoin{{}{}} {} {{}{}}
\pgfsys@moveto{-6.09349pt}{2.7351pt}\pgfsys@lineto{0.34999pt}{0.0pt}\pgfsys@lineto{-6.09349pt}{-2.7351pt}\pgfsys@closepath\pgfsys@fillstroke\pgfsys@endscope}}
}{}{}{{}}\pgfsys@moveto{38.76527pt}{-5.14767pt}\pgfsys@lineto{6.11926pt}{-37.7937pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{-0.7071}{-0.7071}{0.7071}{-0.7071}{6.11926pt}{-37.7937pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} ; ; {{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{}{}{}{{}}\pgfsys@moveto{36.63304pt}{0.0pt}\pgfsys@lineto{9.55438pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{-1.0}{0.0}{0.0}{-1.0}{9.55438pt}{0.0pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} ; ; {{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{}{}{}{{}}\pgfsys@moveto{-36.63301pt}{0.0pt}\pgfsys@lineto{-9.55437pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-9.55437pt}{0.0pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} ; ; {{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{}{}{}{{}}\pgfsys@moveto{-35.73256pt}{-43.91283pt}\pgfsys@lineto{-8.6539pt}{-43.91286pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-8.6539pt}{-43.91286pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} ; ; {{}}{}{{}}
{{{{{}}{}{}{}{}{{}}}}}{}{{{{{}}{}{}{}{}{{}}}}}{{}}{}{}{}{}{}{}{}{{}}\pgfsys@moveto{0.0pt}{-36.63301pt}\pgfsys@lineto{0.0pt}{-9.55437pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.0}{1.0}{-1.0}{0.0}{0.0pt}{-9.55437pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} ; ; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\quad\quad\mathbb{P}(\delta_{i}|x_{i},R(t_{i}),\beta)$
$\displaystyle=\mathbb{E}_{z_{i}\sim
p(.|\theta,x_{i})}\left[\frac{\exp{\beta_{z_{i}}^{T}x_{i}}}{\sum\limits_{j\in
R(t_{i})}\exp{\beta_{z_{j}}^{T}x_{j}}}\right]$ (1)

##  3 Approximate Inference Solution

In this section we provide a computationally tractable approximation for the
maximum-likelihood estimation of the semi-parametric latent variable model
defined in section 2 and in section 4 we show the conditions under which the
approximation is tight. Assuming $z_{i}\perp\\!\\!\\!\perp z_{j,j\neq
i}\,|\,x_{i},x_{j}$ we can define positive random variables (R.V.)
$\alpha_{z_{i}}$ and $\eta_{z_{-i}}$ which are functions of $\beta,R(t_{i})$ &
$z$. This assumption implies that
$\alpha_{z_{i}}\perp\\!\\!\\!\perp\eta_{z_{-i}}\,|\,x_{i},x_{j}$. In the rest
of the paper (unless stated otherwise), the expectation $\mathbb{E}$ is over
the distribution $z_{i}\sim p(z_{i}|\theta,x_{i})$. Using this re-formulation
and the Taylor series expansion we can re-write eq. 1 as,

$\displaystyle\alpha_{z_{i}}=\exp{\beta_{z_{i}}^{T}x_{i}}\quad\eta_{z_{-i}}=\sum\limits_{j\in
R(t_{i})}\exp{\beta_{z_{j}}^{T}x_{j}}\quad
P(\delta_{i}|x_{i},R(t_{i}))=\mathbb{E}\left(\frac{\alpha_{z_{i}}+\eta_{z_{-i}}}{\alpha_{z_{i}}}\right)^{-1}$
(2)
$\displaystyle\mathbb{E}\left[\left(\frac{\alpha_{z_{i}}+\eta_{z_{-i}}}{\alpha_{z_{i}}}\right)^{-1}\right]=\mathbb{E}\left[\sum\limits_{p=1}^{\infty}(-1)^{p}\mathbb{E}(\eta_{z_{-i}}^{p})\mathbb{E}(\alpha_{z_{i}}^{-p})\right]\approx\left[\sum\limits_{p=1}^{K}(-1)^{p}\mathbb{E}(\eta_{z_{-i}}^{p})\mathbb{E}(\alpha_{z_{i}}^{-p})\right]$
(3)

######  Lemma 1.

If we assume $\alpha_{z_{i}}$, $\eta_{z_{-i}}$ to have moments of order
$\mathcal{H}_{\alpha_{z_{i}}},\mathcal{H}_{\eta_{z_{-i}}}$ respectively, then
their ratio distribution will have moments of order
$\frac{\mathcal{H}_{\alpha_{z_{i}}}\mathcal{H}_{\eta_{z_{-i}}}}{\mathcal{H}_{\alpha_{z_{i}}}+\mathcal{H}_{\eta_{z_{-i}}}}$.
(Cedilnik et al., 2006)

Using the Mellin Transform theory for ratio distributions of positive
independent random variables (R.V.); we have
$\mathbb{E}((X/Y)^{p})=\mathbb{E}(X^{p})\mathbb{E}(Y^{-p})$. Based on lemma 1
we limit the expansion in eq. 3 to a finite $K$. At this point, we are
computing expectations over convex functions $x^{p}$ and $x^{-p}$ with $p>0$
(defined on $\mathbb{R}_{++}$). Since for a convex $f$, $E(f(x))\geq f(E(x))$
(Jensen inequality) we can further approximate eq. 3 with eq. 4. Once again we
use the Taylor series approximation to finally arrive at a tractable maximum-
likelihood objective (eq. 5). For each data point the inference complexity
under the original objective is $\mathcal{O}(|z|^{|R(t_{i})|})$ whereas under
the proposed marginalization the complexity reduces to $\mathcal{O}(|z|)$.

$\displaystyle\mathbb{E}\left[\left(\frac{\alpha_{z_{i}}+\eta_{z_{-i}}}{\alpha_{z_{i}}}\right)^{-1}\right]\approx\sum\limits_{p=1}^{K}(\mathbb{E}(\eta_{z_{-i}}))^{p}(\mathbb{E}(\alpha_{z_{i}}))^{-p}\quad\mathbb{P}(\delta_{i}|x_{i},R(t_{i}),\beta)$
$\displaystyle=\left[\frac{\mathbb{E}(\exp{\beta_{z_{i}}^{T}x_{i}})}{\sum\limits_{j\in
R(t_{i})}\mathbb{E}(\exp{\beta_{z_{j}}^{T}x_{j}})}\right]$ (4)

$\displaystyle\mathbb{P}(\delta_{i}|x_{i},R(t_{i}),\beta)$
$\displaystyle=\left[\frac{\mathbb{E}_{z_{i}\sim
p(.|\theta,x_{i})}(\exp{\beta_{z_{i}}^{T}x_{i}})}{\sum\limits_{j\in
R(t_{i})}\mathbb{E}_{z_{j}\sim
p(.|\theta,x_{j})}(\exp{\beta_{z_{j}}^{T}x_{j}})}\right]$ (5)

The crux of the approximation lies in the Jensen inequality. Therefore, we
spend the following section on identifying a tight distribution independent
bound on the inequality gap. If this inequality gap is reduced then we know
that 5 is a good approximation for 1.

##  4 Bounding the Approximation [Analysis]

We identify the conditions under which the approximation 3 is feasible and
provide a closed form bound for it. In order to simplify the statements in the
rest of the paper, we introduce some notations and assumptions here. We assume
that the R.V. $z_{i}\sim p_{\theta}(.|x_{i})$ has mean $\mu(x_{i})$ and
$z_{i}-\mu(x_{i})$ is sub-Gaussian with parameter $\sigma(x_{i})$. This
assumption is fairly common in latent variable models where the true posterior
$p(z_{i}|x_{i})$ is approximated by a Gaussian distribution
$N(\mu(x_{i}),\sigma(x_{i}))$. For a continuous function
$\beta_{z_{i}}=\beta(z_{i})$, with $z_{i}$ lying in a closed, bounded set
(with probability $\delta$), one can bound the values attained by
$\exp{\beta_{z_{i}}^{T}x_{i}}$. Thus for a R.V. $z_{i}$ we obtain reasonable
probabilistic bounds on $\alpha_{z_{i}},\eta_{z_{-i}}$, formalized by the
following statement:
$\alpha_{z_{i}}\in[\mathcal{L}(\alpha_{z_{i}}),\mathcal{U}(\alpha_{z_{i}})],\eta_{z_{-i}}\in[\mathcal{L}(\eta_{z_{-i}}),\mathcal{U}(\eta_{z_{-i}})]$
with probability $\delta$ (defined in theorem 1)).

######  Theorem 1.

Although this is stated for $\alpha_{z_{i}}$, the statement for
$\eta_{z_{-i}}$ is similar.
With probability
$\delta=erf\left(\frac{1}{2}\frac{\mathcal{U}(\alpha_{z_{i}})-\exp{\beta_{z_{i}}^{T}\mu(x_{i})}}{\sqrt{2}\sigma(x_{i})(\exp{\beta_{z_{i}}^{T}\mu(x_{i})})\|\beta_{z_{i}}\|_{2}}\right)-erf\left(\frac{1}{2}\frac{-\mathcal{L}(\alpha_{z_{i}})-\exp{\beta_{z_{i}}^{T}\mu(x_{i})}}{\sqrt{2}\sigma(x_{i})(\exp{\beta_{z_{i}}^{T}\mu(x_{i})})\|\beta_{z_{i}}\|_{2}}\right)$,

$\displaystyle\mathbb{E}(\phi_{p}(\alpha_{z_{i}}))-\phi_{p}(E(\alpha_{z_{i}}))$
$\displaystyle\leq\kappa\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))+(1-\kappa)\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))-\phi_{p}(\kappa\mathcal{U}(\alpha_{z_{i}})+(1-\kappa)\mathcal{L}(\alpha_{z_{i}}))$
(6) $\displaystyle\kappa$
$\displaystyle=\nabla\phi*\left[\frac{\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))-\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))}{\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}})}\right]$
(7)

The conjugate $\phi_{p}^{*}(y)=\sup\limits_{x}y^{T}x-\phi_{p}(x)$,
$\nabla\phi_{p}^{*}(y)=\\{x:y^{T}x-\phi_{p}(x)=\phi_{p}^{*}(y)\\}$, which is a
singleton set for strictly convex functions.

##  5 Joint Objective Function

Since gradients for conjugate functions $\phi_{p}^{*}(x)$ are well defined in
our case (see Appendix A.2), we can show that the approximation in eq. 5 is
good when the joint objective is minimized (eq. 8). The joint objective
enforces the model to find optimal $(\theta^{*},\beta^{*})$ that maximizes the
likelihood in eq. 1, while ensuring proximity to the true objective. Eq. 8 can
also be viewed from the perspective of a regularized objective where the model
learns to enforce additional constraints on the variance of
$\alpha_{z_{i}},\eta_{z_{-i}}$ and thus ends up with distributions of $z_{i}$
with rapidly decaying Gaussian tails.

$\displaystyle(\beta^{*},\theta^{*})=\mathop{\mathrm{argmin}}_{\beta,\theta}\sum_{i=1}^{N}\left[\frac{\mathbb{E}_{z_{i}\sim
p(.|\theta,x_{i})}(\exp{\beta_{z_{i}}^{T}x_{i}})}{\sum\limits_{j\in
R(t_{i})}\mathbb{E}_{z_{j}\sim
p(.|\theta,x_{j})}(\exp{\beta_{z_{j}}^{T}x_{j}})}+\sum\limits_{j\in
R(t_{i})}L(x_{j},\kappa_{j},\theta,\beta)\right]$ (8)

Here $L(x_{j},\kappa_{j},\theta,\beta)$ is obtained by using theorem 1 which
bounds the Jensen’s inequality for each data point $(\delta_{i},x_{i},t_{i})$,
via a sum over gradients computed for the functions $\phi_{p}^{*}$.

##  6 Results

We analyze two types of results: (1) we evaluate our combined objective (eq.
8) on a proportional hazards (CPH) model and show an improvement in the
concordance-index (table 1), (2) we compare our proposed distribution agnostic
bound against a standard bound for the Jensen inequality (Dragomir, 1999).

Figure 2: Comparisons of the upper bounds on the Jensen inequality established
by our model and the baseline (Dragomir, 1999). 95% confidence intervals are
shown when values of $\mathcal{U}(\alpha_{z_{i}})\sim N(\mu_{u},\sigma^{2})$
and $\mathcal{L}(\alpha_{z_{i}})\sim N(\mu_{l},\sigma^{2})$.
Figure 3: Visualization of our distribution agnostic bound on Jensen’s
inequality for strictly convex functions.

###  6.1 Survival Analysis

Given a discrete $z$ we model the distribution $z_{i}\sim p(,|\theta,x_{i})$
to be a multinomial. The final layer of the input encoder network is a softmax
operation ensuring that the distribution over the latent space of $z$ is a
valid one. We compare our models: Latent Variable CPH via Hard/Soft Gating
(LV-CPH-HG/LV-CPH-SG) against three popular baselines: Cph (Cox, 1955), Rsf
(Ishwaran and Lu, 2007), DeepSurv (Katzman et al., 2016) on common datasets in
survival analysis: METABRIC (Yao, 2014), ROTTERDAM-GBSG (Schumacher et al.,
1994), SUPPORT (Knaus et al., 1995).

For the discrete case, it is easy to see that the regularizer
$L(x_{i},\kappa_{i},\theta,\beta)$ in eq. 8 is minimized when $z_{i}$ has low
variance (or entropy). We enforce a low entropy distribution by gating
(soft/hard) the predictions $p(.|\theta,x_{i})$ obtained from the softmax
layer. Since for the discrete case, low entropy ($\mathbb{H}_{z_{i}}$) on
$z_{i}$ $\Longleftrightarrow$
$|\mathcal{L}(\alpha_{z_{i}})-\mathcal{U}(\alpha_{z_{i}})|+|\mathcal{L}(\eta_{z_{-i}})-\mathcal{U}(\eta_{z_{-i}})|<C$
and $C\to 0$ as $\mathbb{H}_{z_{i}}\to 0$, one can instead minimize
$\mathbb{H}_{z_{i}}$ to effectively reduce the upper bound in theorem 1.
Therefore, we conclude that optimizing for the joint objective function in eq.
8 instead of the mere approximation in eq. 5 leads to an improved concordance-
index for CPH models.

###  6.2 Tightness of the proposed bound

Figure 2 compares the bound computed by (Dragomir, 1999) [baseline] (Appendix
A.1) against our tight bound (Appendix A.2), by sampling
$\mathcal{U}(\alpha_{z_{i}}),\mathcal{U}(\eta_{z_{-i}}),\mathcal{L}(\alpha_{z_{i}}),\mathcal{L}(\eta_{z_{-i}})$
from distinct fixed normal distributions. Our bound is much tighter for
smaller values of $p$, and it converges to the baseline’s value for large p.
Looking at figure 3 it is easy to verify that the bound we propose on Jensen’s
inequality is stronger than the baseline.

Model | Metabric | Rotterdam-Gbsg | Support
---|---|---|---
Cph | $0.6306\pm 0.004$ | $0.6578\pm 0.004$ | $0.5828\pm 0.002$
DeepSurv | $0.6434\pm 0.004$ | 0.$6684\pm 0.003$ | $0.6183\pm 0.002$
Rsf | $0.6243\pm 0.004$ | $0.6512\pm 0.003$ | $0.6130\pm 0.002$
LV-CPH-SG | $\mathbf{0.6585\pm 0.003}$ | $0.6752\pm 0.002$ | $\mathbf{0.6196\pm 0.001}$
LV-CPH-HG | $0.6349\pm 0.003$ | $\mathbf{0.6866\pm 0.002}$ | $0.5706\pm 0.001$
Table 1: Results of hard and soft linear gating networks and their comparison
with relevant baselines (95% bootstrap CI ). (Nagpal et al., 2019)

##  7 Discussion

We propose an approximation for the MLE objective in TPPs involving a partial
likelihood function with latent factors. We also show that the MLE
approximation can be bounded by minimizing a joint objective which includes an
upper bound on the approximation gap. We have shown this to be theoretically
and empirically better for the partial likelihood estimation (in survival
analysis). Future work on this inference method would be to further exploit
the tractable closed form approximation gap by directly optimizing for it with
iterative methods like ADAM. Yet another direction would be to extend this
work to latent variable models for semi-parametric models like Gaussian
processes for survival analysis (Fernández et al., 2016).

## References

  * Cedilnik et al. [2006] Anton Cedilnik, Katarina Kosmelj, and Andrej Blejec.  Ratio of two random variables: A note on the existence of its moments.  _Metodoloski zvezki_ , 3(1):1, 2006.
  * Cox [1955] David R Cox.  Some statistical methods connected with series of events.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 17(2):129–157, 1955.
  * Diggle [2005] Peter J Diggle.  A partial likelihood for spatio-temporal point processes.  _Biostats_ , 2005.
  * Dragomir [1999] Silvestru S Dragomir.  A converse result for jensen’s discrete inequality via grüss’ inequality and applications in information theory.  _An. Univ. Oradea Fasc. Mat_ , 1999.
  * Fernández et al. [2016] Tamara Fernández, Nicolás Rivera, and Yee Whye Teh.  Gaussian processes for survival analysis.  In _Advances in Neural Information Processing Systems_ , pages 5021–5029, 2016.
  * Ishwaran and Lu [2007] Hemant Ishwaran and Min Lu.  Random survival forests.  _Wiley StatsRef: Statistics Reference Online_ , pages 1–13, 2007\.
  * Jebara and Choromanska [2012] Tony Jebara and Anna Choromanska.  Majorization for crfs and latent likelihoods.  In _Advances in Neural Information Processing Systems_ , pages 557–565, 2012.
  * Katzman et al. [2016] Jared L Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger.  Deep survival: A deep cox proportional hazards network.  _stat_ , 1050:2, 2016.
  * Knaus et al. [1995] William A Knaus, Frank E Harrell, Joanne Lynn, Lee Goldman, Russell S Phillips, Alfred F Connors, Neal V Dawson, William J Fulkerson, Robert M Califf, Norman Desbiens, et al.  The support prognostic model: objective estimates of survival for seriously ill hospitalized adults.  _Annals of internal medicine_ , 122(3):191–203, 1995.
  * Linderman and Adams [2014] Scott Linderman and Ryan Adams.  Discovering latent network structure in point process data.  In _International Conference on Machine Learning_ , pages 1413–1421, 2014.
  * Nagpal et al. [2019] Chirag Nagpal, Rohan Sangave, Amit Chahar, Parth Shah, Artur Dubrawski, and Bhiksha Raj.  Nonlinear semi-parametric models for survival analysis.  _arXiv preprint arXiv:1905.05865_ , 2019.
  * Quattoni et al. [2007] Ariadna Quattoni, Sybor Wang, Louis-Philippe Morency, Michael Collins, and Trevor Darrell.  Hidden conditional random fields.  _IEEE Transactions on Pattern Analysis & Machine Intelligence_, pages 1848–1852, 2007.
  * Rosen and Tanner [1999] Ori Rosen and Martin Tanner.  Mixtures of proportional hazards regression models.  _Statistics in Medicine_ , 18(9):1119–1131, 1999\.
  * Ruder [2016] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  _arXiv preprint arXiv:1609.04747_ , 2016.
  * Samo and Roberts [2015] Yves-Laurent Kom Samo and Stephen Roberts.  Scalable nonparametric bayesian inference on point processes with gaussian processes.  In _International Conference on Machine Learning_ , pages 2227–2236, 2015.
  * Schumacher et al. [1994] M Schumacher, G Bastert, H Bojar, K Huebner, M Olschewski, W Sauerbrei, C Schmoor, C Beyerle, RL Neumann, and HF Rauschecker.  Randomized 2 x 2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients. german breast cancer study group.  _Journal of Clinical Oncology_ , 12(10):2086–2093, 1994.
  * Simic [2009] Slavko Simic.  On an upper bound for jensen’s inequality.  _Journal of Inequalities in Pure and Applied Mathematics_ , 2009.
  * Snoek et al. [2013] Jasper Snoek, Richard Zemel, and Ryan P Adams.  A determinantal point process latent variable model for inhibition in neural spiking data.  In _Advances in Neural Information Processing Systems_ , pages 1932–1940, 2013.
  * Sutton et al. [2012] Charles Sutton, Andrew McCallum, et al.  An introduction to conditional random fields.  _Foundations and Trends® in Machine Learning_ , 4(4):267–373, 2012.
  * Yao [2014] Cindy Q Yao.  _Methods for the Identification of Biomarkers in Prostate and Breast Cancer_.  PhD thesis, University of Toronto (Canada), 2014.

##  Appendix A Appendix

###  A.1 Bounding Jensen’s Inequality [Dragomir – Loose Bound]

[Simic, 2009] propose multiple distribution agnostic upper bounds for Jensen’s
inequality in the case of generic continuous convex functions defined on a
compact set $[\mathcal{L},\mathcal{U}]$. One of the popular bounds in this
regime is the Dragomir’s inequality proposed in [Dragomir, 1999]. Given
$z_{i}$ is bounded with probability $\delta$, section 4 bounds
$\alpha_{z_{i}},\eta_{z_{-i}}$ with
$(\alpha_{z_{i}}\in[\mathcal{L}(\alpha_{z_{i}}),\mathcal{U}(\alpha_{z_{i}})],\eta_{z_{-i}}\in[\mathcal{L}(\eta_{z_{-i}}),\mathcal{U}(\eta_{z_{-i}})])$.
By Dragomir’s [Dragomir, 1999] inequality for a convex function $f$,

$\displaystyle\mathbb{E}[f(\alpha_{z_{i}})]-f(\mathbb{E}(\alpha_{z_{i}}))\leq\frac{1}{4}(\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}}))(f^{\prime}(\mathcal{U}(\alpha_{z_{i}}))-f^{\prime}(\mathcal{L}(\alpha_{z_{i}})))$
(9)

For $f\in\\{f_{p}\\}_{p=2}^{\infty},f_{p}(x)=x^{p}$ and
$f\in\\{f_{p}\\}_{p=2}^{\infty},f_{p}(x)=x^{-p}$, the bounds are given by eq.
10 and eq. 11 respectively.

$\displaystyle\frac{p}{4}(\mathcal{U}(\eta_{z_{-i}})-\mathcal{L}(\eta_{z_{-i}}))(\mathcal{U}(\eta_{z_{-i}})^{p-1}-\mathcal{L}(\eta_{z_{-i}})^{p-1})$
(10)
$\displaystyle\frac{p}{4}(\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}}))(\mathcal{L}(\alpha_{z_{i}})^{-(p+1)}-\mathcal{U}(\alpha_{z_{i}})^{-(p+1)})$
(11)

This bound is easy to compute and can be visualized via the gap shown in
figure 3. This is also quite naive and generic since it only uses the first
order conditions for convex functions to arrive at an upper bound. In the
following section we provide a tighter upper bound under the stronger
assumptions of strict convexity.

###  A.2 Bounding Jensen’s Inequality [Convex Conjugate – Tight Bound]

This section provides the proof for theorem 1 in the main paper. We
investigate bounds under the special case of strictly convex functions
$\\{\phi_{p}\\}_{|p|=2}^{\infty}$ with
$\phi_{p}:\mathbb{R}_{++}\to\mathbb{R}_{++}$ and $\phi_{p}(x)=x^{p}$. We also
show visually (figure 3) that our bound is the tightest possible distribution
agnostic bound for the given set of functions. With
$\mathbb{G}(\alpha_{z_{i}},\phi_{p})$ as the bound of interest,

$\displaystyle\mathbb{G}(\alpha_{z_{i}},\phi_{p})\leq\max_{\kappa\in[0,1]}\kappa\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))+(1-\kappa)\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))-\phi_{p}(\kappa\mathcal{U}(\alpha_{z_{i}})+(1-\kappa)\mathcal{L}(\alpha_{z_{i}}))$
(12)

Figure 3 depicts geometrically the maximization problem in the RHS of eq. 12
which we solve via the convex conjugate of $\phi_{p}$. The optimization
problem in eq. 12 involves identifying
$\alpha_{z_{i}}\in[\mathcal{L}(\alpha_{z_{i}}),\mathcal{U}(\alpha_{z_{i}})]$,
such that the line in figure 3 denoted by
$y=\phi^{\prime}(\alpha_{z_{i}}^{*})\alpha_{z_{i}}+c$ where
$\alpha_{z_{i}}^{*}$ is farthest away from $\phi(\alpha_{z_{i}})$ at the
optimal point.

Given that $\phi_{p}(x)=(1/x^{p})$ is a closed proper strict convex function
on $[\mathcal{L}(\alpha_{z_{i}}),\mathcal{U}(\alpha_{z_{i}})]$, (the case when
we consider $\phi_{p}(x)=x^{p}$ on
$[\mathcal{L}(\eta_{z_{-i}}),\mathcal{U}(\eta_{z_{-i}})]$ is similar) we can
define the convex conjugate of $\phi_{p}(x)$. The form in eq. 13 is similar to
what we need (distance between $\phi_{p}(x)$ and a line with slope defined by
$\frac{\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))-\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))}{\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}})}$

$\displaystyle\phi_{p}^{*}(y)=\sup\limits_{x}yx-\phi_{p}(x)\quad\partial\phi_{p}^{*}(y)=\\{x:yx-\phi_{p}(x)=\phi_{p}^{*}(y)\\}$
(13)
$\displaystyle\partial\phi_{p}^{*}\left(\frac{\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))-\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))}{\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}})}\right)=\mathop{\mathrm{argmax}}_{w\in[\mathcal{L}(\alpha_{z_{i}}),\mathcal{U}(\alpha_{z_{i}})]}w.\frac{\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))-\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))}{\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}})}-\phi_{p}(w).$
(14)

Since $(1/x^{p})$ is strictly convex the subgradient of $\phi_{p}^{*}$ at $y$
is a singleton set and is exactly equal to the $\partial\phi_{p}^{*}(y)$,
giving us a unique maximizer $w^{*}$ in eq. 14.

$\displaystyle\partial\phi_{p}^{*}(y)=\left(\frac{-p}{y}\right)^{\frac{1}{p+1}}\quad
w^{*}=\left(p.\frac{\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}})}{\mathcal{L}(\alpha_{z_{i}})^{-p}-\mathcal{U}(\alpha_{z_{i}})^{-p}}\right)^{\frac{1}{p+1}}$
(15)

Eq. 15 tells us that $w^{*}$ (and thus the bound in Eq. 12) for a given sample
$(x_{i},t_{i},R(t_{i}))$ is a function of $\beta,x_{i}$ & $x_{k}\in R(t_{i})$.
Given $w*$ it is easy to compute $\mathbb{G}(\alpha_{z_{i}},\phi_{p})$ by,

$\displaystyle
w^{*}=\kappa^{*}\mathcal{U}(\alpha_{z_{i}})+(1-\kappa^{*})\mathcal{L}(\alpha_{z_{i}})$
(16)
$\displaystyle\kappa^{*}=\frac{w^{*}-\mathcal{L}(\alpha_{z_{i}})}{\mathcal{U}(\alpha_{z_{i}})-\mathcal{L}(\alpha_{z_{i}})}$
(17)
$\displaystyle\mathbb{G}(\alpha_{z_{i}},\phi_{p})=\kappa^{*}\phi_{p}(\mathcal{U}(\alpha_{z_{i}}))$
$\displaystyle+(1-\kappa^{*})\phi_{p}(\mathcal{L}(\alpha_{z_{i}}))-\phi_{p}(\kappa^{*}\mathcal{U}(\alpha_{z_{i}})+(1-\kappa^{*})\mathcal{L}(\alpha_{z_{i}}))$
(18)

In section 6 we empirically compared this bound with eqs. 10, 11 from appendix
A.1.