# Influence Estimation and Maximization via Neural Mean-Field Dynamics

Shushan He Department of Mathematics and Statistics, Georgia State University,
Atlanta, Georgia 30303, USA. (she4@student.gsu.edu). Hongyuan Zha School of
Data Science, Shenzhen Research Institute of Big Data, The Chinese University
of Hong Kong, Shenzhen, Guangdong, China, 518172. (zhahy@cuhk.edu.cn).
Xiaojing Ye Department of Mathematics and Statistics, Georgia State
University, Atlanta, Georgia 30303, USA. (xye@gsu.edu).

###### Abstract

We propose a novel learning framework using neural mean-field (NMF) dynamics
for inference and estimation problems on heterogeneous diffusion networks. Our
new framework leverages the Mori-Zwanzig formalism to obtain an exact
evolution equation of the individual node infection probabilities, which
renders a delay differential equation with memory integral approximated by
learnable time convolution operators. Directly using information diffusion
cascade data, our framework can _simultaneously_ learn the structure of the
diffusion network and the evolution of node infection probabilities.
Connections between parameter learning and optimal control are also
established, leading to a rigorous and implementable algorithm for training
NMF. Moreover, we show that the projected gradient descent method can be
employed to solve the challenging influence maximization problem, where the
gradient is computed extremely fast by integrating NMF forward in time just
once in each iteration. Extensive empirical studies show that our approach is
versatile and robust to variations of the underlying diffusion network models,
and significantly outperform existing approaches in accuracy and efficiency on
both synthetic and real-world data.

Keywords— Diffusion networks, influence estimation, Mori-Zwanzig formalism,
influence maximization

##  1 Introduction

Continuous-time information diffusion on heterogenous networks is a prevalent
phenomenon [2, 44, 48]. News spreading on social media [13, 15, 56], viral
marketing [26, 27, 58], computer malware propagation, and epidemics of
contagious diseases [1, 42, 48, 53] are all examples of diffusion on networks,
among many others. For instance, a piece of information (such as a tweet) can
be retweeted by users (nodes) with followee-follower relationships (edge) on
the Twitter network. We call a user _infected_ if she retweets, and her
followers see her retweet and can also become infected if they retweet in
turn, and so on. Such information diffusion mimics the epidemic spread where
an infectious virus can spread to individuals (human, animal, or plant) and
then to many others upon their close contact. The study of heterogeneous
diffusion networks only emerged in the past decade and is considered very
challenging, mainly because of the extremely large scale of modern networks,
the heterogeneous inter-dependencies between the nodes, and the randomness
exhibited in cascade data.

In the remainder of this section, we provide the mathematical formulations of
the inference, influence estimation, and influence maximization problems on an
arbitrary diffusion network. Throughout this paper, we use boldfaced lower
(upper) letter to denote vector (matrix) or vector-valued (matrix-valued)
function, and $(\cdot)_{k}$ (or $(\cdot)_{ij}$) for its $k$th component (or
$(i,j)$-th entry). All vectors are column vectors unless otherwise noted. We
follow the Matlab syntax and use $[\bm{x};\bm{y}]$ to denote the vector that
stacks $\bm{x}$ and $\bm{y}$ vertically, and $\bm{x}\cdot\bm{y}$ or
$\bm{x}^{\top}\bm{y}$ for the inner product. Time is denoted by $t$ in either
continuous ($t\in[0,T]$) or discrete case ($t=0,1,\dots,T$) for some time
horizon $T\in\mathbb{R}_{+}$ ($\mathbb{N}$ in discrete case). Derivative ′ is
with respect to $t$, and gradient $\nabla_{\bm{x}}$ is with respect to
$\bm{x}$. Probability is denoted by $\mathrm{Pr}(\cdot)$, and expectation with
respect to $X$ (or its distribution function $p_{X}$) is denoted by
$\mathbb{E}_{X}[\,\cdot\,]$. The $n$-vectors
$\bm{1}_{n},\bm{0}_{n}\in\mathbb{R}^{n}$ stand for the vectors of ones and
zeros respectively, and we often omit the subscript $n$ when their dimensions
are obvious from the context.

###  1.1 Diffusion network models

Consider a diffusion network model, which consists of a network (directed
graph) $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with node set
$\mathcal{V}=[n]\mathrel{\mathop{\ordinarycolon}}=\\{1,\dots,n\\}$ and edge
set $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$, and a _diffusion model_
that describes the distribution $p(t;\alpha_{ij})$ of the time $t$ that an
infected node $i$ takes to infect her healthy neighbor
$j\in\\{j^{\prime}\mathrel{\mathop{\ordinarycolon}}(i,j^{\prime})\in\mathcal{E}\\}$.
Here $\alpha_{ij}$ is the infection rate of $i$ on $j$ which _vary across
different edges_. That is, $t_{ij}$ is a random variable following the density
function $p(t;\alpha_{ij})$ for each $(i,j)\in\mathcal{E}$. We assume that the
infection is _progressive_ , i.e., a node will not be infected again nor
recover once infected, since generalization to the case with recovery is
straightforward. Then, given a source set $\mathcal{S}$ (a subset of
$\mathcal{V}$) of nodes that are infected at time $0$, they will infect their
healthy neighbors with random infection times described above; and the
infected neighbors will then infect their healthy neighbors, and so on. As
such, the infection initiated by $\mathcal{S}$ at time $0$ propagates to other
nodes of the network. We call one course of such propagation a _cascade_. For
simplicity, it is common to assume that the infection times across different
edges are independent, known as the _continuous-time independent cascade_
(CIC) model [21, 13, 19].

$t_{1}$$t_{2}$$t_{3}$$t_{4}$$t_{5}$$t_{6}$$t_{7}$$t_{8}$ Figure 1: Example of
a sample cascade on a diffusion network. The cascade was originated from the
source set $\mathcal{S}=\set{1}$ and gradually propagates to other nodes
through their directed edge connections. The time line below the network shows
the wall-clock time $t_{i}$ that each node $i$ was infected during the cascade
with $t_{1}=0$. The orange edges indicate whom each node got infection from,
and $t_{ij}\mathrel{\mathop{\ordinarycolon}}=t_{j}-t_{i}$ is the time that
node $i$ took to infect node $j$.

In Figure 1, we illustrate one such cascade originated from a singleton source
set $\mathcal{S}=\set{1}$, which spreads to other nodes during the
propagation. The orange edges indicate whom a node got infection from, for
example, node 4 succeeded in infecting node 6 before node 1 did. The time line
below the network indicates the wall-clock time $t_{i}$ of each node $i$ got
infected in this cascade. In particular, $t_{1}=0$. Moreover,
$t_{ij}\mathrel{\mathop{\ordinarycolon}}=t_{j}-t_{i}$ is the time node $i$
took to infect node $j$. Note that this is one sample cascade of
$\mathcal{S}=\set{1}$, and a different sample cascade of the same source
$\mathcal{S}$ may yield different infected nodes and infection times due to
the randomness of $t_{ij}$.

The standard diffusion model with exponential distribution $p(t;\alpha)=\alpha
e^{-\alpha t}$ is mostly widely used in the literature. That is, $t_{ij}\sim
p(t;\alpha_{ij})$ for each $(i,j)\in\mathcal{E}$. Note that the parameter
$\alpha_{ij}>0$ in the exponential distribution indicates the _strength_ of
impact node $i$ has on $j$—the expectation of $t_{ij}\sim p(t;\alpha_{i}j)$ is
$1/\alpha_{ij}$—and the larger $\alpha_{ij}$ is, the sooner node $j$ will be
infected by $i$ on expectation. We focus on the diffusion model with
exponential distribution in this work. Other distributions, such as Rayleigh
and general Weibull distributions, are also experimented in our empirical
studies in this work.

###  1.2 Cascade data

Observation data $\mathcal{D}$ of a diffusion network are often in the form of
sample cascades
$\mathcal{D}\mathrel{\mathop{\ordinarycolon}}=\\{\mathcal{C}_{k}=(\mathcal{S}_{k},\bm{\tau}_{k})\in\mathcal{V}\times\mathbb{R}_{+}^{n}\mathrel{\mathop{\ordinarycolon}}k\in[K]\\}$,
where the $k$th cascade $\mathcal{C}_{k}$ records its source set
$\mathcal{S}_{k}\subset\mathcal{V}$ and the time $(\bm{\tau}_{k})_{i}\geq 0$
which indicates when node $i$ was infected (if $i$ was not infected during
$\mathcal{C}_{k}$ then $(\bm{\tau}_{k})_{i}=\infty$). See Figure 1 for one of
such sample cascades, where we have
$\bm{\tau}=\\{t_{1},\dots,t_{8},\infty,\dots,\infty\\}$ if no other nodes were
infected in this cascade. Cascade data are collected from historical events
for training purposes.

###  1.3 Network inference and influence estimation

Suppose $\mathcal{G}=(\mathcal{V},\mathcal{E})$ is a diffusion network with
transmission matrix $\bm{A}$, where $(\bm{A})_{ji}=\alpha_{ij}$ is the
parameter of $p(t;\alpha_{ij})$ for edge $(i,j)$. Then the goal of _infection
probability estimation_ (_influence estimation_ , or _influence prediction_ ,
for short) is to compute

$\bm{x}(t;{\bm{\chi}}_{\mathcal{S}})=[x_{1}(t;{\bm{\chi}}_{\mathcal{S}}),\dots,x_{n}(t;{\bm{\chi}}_{\mathcal{S}})]^{\top}\in[0,1]^{n}$
(1)

for all time $t>0$ and any given source set $\mathcal{S}\subset\mathcal{V}$.
In (1), $x_{i}(t;{\bm{\chi}}_{\mathcal{S}})$ is the probability of node $i$
being infected at time $t$ given the source set $\mathcal{S}$, and
${\bm{\chi}}_{\mathcal{S}}\in\set{0,1}^{n}$ indicates the identities of
$\mathcal{S}$, i.e., $({\bm{\chi}}_{\mathcal{S}})_{i}=1$ if $i\in\mathcal{S}$
and $0$ otherwise. Note that we use ${\bm{\chi}}_{\mathcal{S}}$ and
$\mathcal{S}$ interchangeably hereafter. The probability
$\bm{x}(t;{\bm{\chi}}_{\mathcal{S}})$ can also be used to compute the
_influence_ function
$\sigma(t;{\bm{\chi}}_{\mathcal{S}})\mathrel{\mathop{\ordinarycolon}}=\bm{1}_{n}^{\top}\bm{x}(t;{\bm{\chi}}_{\mathcal{S}})$,
the expected number of infected nodes at time $t$. Our method can be readily
applied to influence functions defined with uneven weights (rather than 1’s on
all nodes) if the severity of infection varies at different nodes, but we
focus on the even weight case for the sake of conciseness.

Most influence estimation problems do not assume knowledge of $\bm{A}$.
Instead, only cascade data $\mathcal{D}$ are available. In this case, _network
inference_ is often needed. Network inference refers to the problem of
uncovering $\mathcal{E}$ and $\bm{A}$ from cascade data $\mathcal{D}$, and is
of independent research interests in the literature. Now influence estimation
can be tackled by a _two-stage_ approach, where a network inference is
performed first to learn the network structure $\mathcal{E}$ and the diffusion
model parameters $\bm{A}$, and then an influence estimation is used to compute
the influence of the source set $\mathcal{S}$. However, both influence
estimation and network inference problems are very challenging, and the
approximation errors and biases in the two stages will certainly accumulate.
Alternately, one can use a _one-stage_ approach to directly estimate
$\bm{x}(t;{\bm{\chi}}_{\mathcal{S}})$ of any $\mathcal{S}$ from the cascade
data $\mathcal{D}$, which is more versatile and less prone to diffusion model
misspecification. Our method is a such kind of one-stage method. Additionally,
it allows knowledge of $\mathcal{E}$ and/or $\bm{A}$, if available, to be
integrated for further performance improvement.

###  1.4 Influence maximization

Given a budget size $n_{0}\in\\{1,\dots,n-1\\}$, the goal of _influence
maximization_ is to find the source set $\mathcal{S}$ which generates the
maximal influence $\sigma(t;{\bm{\chi}}_{\mathcal{S}})$ at a prescribed time
$t$ among all source sets of size $n_{0}$. Influence maximization can be
formulated as follows:

$\max_{{\bm{\chi}}_{\mathcal{S}}}\
\sigma(t;{\bm{\chi}}_{\mathcal{S}}),\quad\mathrm{s.t.}\quad{\bm{\chi}}_{\mathcal{S}}\in\\{0,1\\}^{n},\quad\bm{1}_{n}^{\top}{\bm{\chi}}_{\mathcal{S}}=n_{0}.$
(2)

There are two main ingredients of an influence maximization method for solving
(2): an influence estimation subroutine that evaluates the influence
$\sigma(t;{\bm{\chi}}_{\mathcal{S}})$ for any given source set $\mathcal{S}$,
and an (approximate) combinatorial optimization solver to find the optimal set
$\mathcal{S}$ of (2) that repeatedly calls the subroutine. The combinatorial
optimization problem is NP-hard and is often approximately solved by greedy
algorithms with guaranteed sub-optimality when
$\sigma(t;{\bm{\chi}}_{\mathcal{S}})$ is submodular in
${\bm{\chi}}_{\mathcal{S}}$. In this work, we propose a variational method
based on the continuous relaxation of (2), and show that it can tackle this
challenging problem very efficiently using our solution framework.

###  1.5 Summary of contribution

In this paper, we develop a comprehensive framework, called neural mean-field
(NMF) dynamics, for simultaneous influence estimation and network inference
from cascade data on a diffusion network. We substantially extend our
preliminary work [25] which first proposed NMF for influence estimation and
network influence with a discrete-time setting. The novelty and contribution
of the present work in contrast to existing ones, including [25], are
summarized as follows:

  1. 1.

We extend the NMF dynamics developed in [25] to the continuous-time setting
which is more suitable for real-world applications of diffusion networks. We
show that the continuous-time NMF dynamics can be naturally incorporated into
the likelihood function of the corresponding point process, which in turn
plays the role of loss function, whereas [25] directly uses cross-entropy of
the observed discrete-time data as the loss function.

  2. 2.

We prove rigorously the connections between parameter learning in continuous-
time NMF and optimal control, where the NMF parameter serves as the time
invariant control. The derivations lead to a fast algorithm based on numerical
ordinary differential equation (ODE) solver that is easy to implement. Unlike
the standard deep residual network training used in [25], we prove rigorously
that the gradients in continuous-time NMF training can be efficiently computed
by solving the ODE defined by NMF forward in time and an augmented co-state
ODE backward in time.

  3. 3.

Based on our continuous-time NMF framework, we develop a highly efficient
algorithm for the very challenging influence maximization problem. In each
iteration, our algorithm only requires solving one augmented ODE based on NMF
dynamics forward in time and one quadratic program, both of which can be
computed very efficiently.

All the theoretical and algorithm developments mentioned above are supported
by extensive empirical studies in this work. The numerical results show that
our approach is robust to the variation of the unknown underlying diffusion
models, and it also significantly outperforms existing approaches on both
synthetic and real-world diffusion networks.

###  1.6 Paper outline

The remainder of this paper is organized as follows. In Section 2, we develop
the proposed neural mean-field dynamics for network inference and influence
estimation on diffusion networks, as well as an optimal control formulation
for parameter training. We show that our new solution framework leads to an
efficient influence maximization algorithm in Section 3. We demonstrate the
performance of the proposed method on influence estimation and maximization on
a variety of synthetic and real-world networks in Section 4. A comprehensive
review of related work in the literature is provided in Section 5. Section 6
concludes the paper.

##  2 Neural Mean-Field Dynamics

###  2.1 Modeling diffusion by stochastic jump processes

We begin with the jump process formulation of network diffusion. Given a
source set ${\bm{\chi}}_{\mathcal{S}}$, let
$X_{i}(t;{\bm{\chi}}_{\mathcal{S}})$ denote the infection status of the node
$i$ at time $t$. Namely, $X_{i}(t)=1$ if node $i$ is infected by time $t$, and
$0$ otherwise. Then $\\{X_{i}(t)\mathrel{\mathop{\ordinarycolon}}i\in[n]\\}$
are a set of $n$ coupled jump processes, such that
$X_{i}(t;{\bm{\chi}}_{\mathcal{S}})$ jumps from $0$ to $1$ when the node $i$
is infected at $t$. Let $\lambda_{i}^{*}(t)$ be the conditional intensity of
$X_{i}(t;{\bm{\chi}}_{\mathcal{S}})$ given history
$\mathcal{H}(t)=\\{X_{i}(s;{\bm{\chi}}_{\mathcal{S}})\mathrel{\mathop{\ordinarycolon}}s\leq
t,\,i\in[n]\\}$, i.e.,

$\lambda_{i}^{*}(t)\mathrel{\mathop{\ordinarycolon}}=\lim_{\tau\to
0^{+}}\frac{\mathbb{E}[X_{i}(t+\tau;{\bm{\chi}}_{\mathcal{S}})-X_{i}(t;{\bm{\chi}}_{\mathcal{S}})|\mathcal{H}(t)]}{\tau}.$
(3)

In influence estimation, our goal is to compute the probability
$\bm{x}(t;{\bm{\chi}}_{\mathcal{S}})=[x_{i}(t;{\bm{\chi}}_{\mathcal{S}})]$ in
(1), which is the expectation of $X_{i}(t;{\bm{\chi}}_{\mathcal{S}})$
conditioning on $\mathcal{H}(t)$:

$x_{i}(t;{\bm{\chi}}_{\mathcal{S}})=\mathbb{E}_{\mathcal{H}(t)}[X_{i}(t;{\bm{\chi}}_{\mathcal{S}})|\mathcal{H}(t)].$
(4)

To this end, we adopt the following notations (for notation simplicity we
temporarily drop ${\bm{\chi}}_{\mathcal{S}}$ in this subsection as the source
set $\mathcal{S}$ is arbitrary but fixed):

$x_{I}(t)=\mathbb{E}_{\mathcal{H}(t)}\mathinner{\bigl{[}\textstyle\prod\nolimits_{i\in
I}X_{i}(t;{\bm{\chi}}_{\mathcal{S}})\big{|}\mathcal{H}(t)\bigr{]}},\quad
y_{I}(t)=\textstyle\prod\nolimits_{i\in I}x_{i}(t),\quad
e_{I}(t)=x_{I}(t)-y_{I}(t)$ (5)

for all $I\subset[n]$ and $|I|\geq 2$. Then we can derive the evolution of
$\bm{z}\mathrel{\mathop{\ordinarycolon}}=[\bm{x};\bm{e}]$. Here
$\bm{x}(t)\in[0,1]^{n}$ is the _resolved_ variable whose value is of interests
and samples can be observed in cascade data $\mathcal{D}$, and
$\bm{e}(t)=[\cdots;e_{I}(t);\dots]\in\mathbb{R}^{2^{n}-n-1}$ is the
_unresolved_ variable. The evolution of $\bm{z}$ is given by the following
theorem, and the proof is provided in Section A.1.

######  Theorem 1.

The evolution of $\bm{z}(t)=[\bm{x}(t);\bm{e}(t)]$ follows the nonlinear
differential equation:

$\bm{z}^{\prime}=\bar{\bm{f}}(\bm{z}),\quad\mbox{where}\quad\bar{\bm{f}}(\bm{z})=\bar{\bm{f}}(\bm{x},\bm{e})=\begin{bmatrix}\bm{f}(\bm{x};\bm{A})-(\bm{A}\odot\bm{E})\bm{1};\
\cdots,f_{I}(\bm{x},\bm{e});\cdots\end{bmatrix},$ (6)

with initial value $\bm{z}_{0}=[{\bm{\chi}}_{\mathcal{S}};\bm{0}]$,
$\bm{E}=[e_{ij}]\in\mathbb{R}^{n\times n}$, and

$\displaystyle\bm{f}(\bm{x};\bm{A})$
$\displaystyle=\bm{A}\bm{x}-\operatorname{diag}(\bm{x})\bm{A}\bm{x},$ (7)
$\displaystyle f_{I}(\bm{x},\bm{e})$ $\displaystyle=\sum_{i\in I}\sum_{j\notin
I}\alpha_{ji}(y_{I}-y_{I\cup\\{j\\}}+e_{I}-e_{I\cup\\{j\\}})-\sum_{i\in
I}y_{I\setminus\\{i\\}}\sum_{j\neq i}\alpha_{ji}(x_{j}-y_{ij}-e_{ij}).$ (8)

We remark that the dimension of $\bm{z}$ is $2^{n}-1$ which grows
exponentially fast in $n$ and hence renders the computation infeasible in
practice. To overcome this issue, we employ the Mori-Zwanzig formalism [6] to
derive a reduced-order model of $\bm{x}$ that has dimensionality $n$ only, as
shown in the next subsection.

###  2.2 Mori-Zwanzig memory closure

We employ the Mori-Zwanzig (MZ) formalism[6] that allows to introduce a
generalized Langevin equation (GLE) of the $\bm{x}$ part of the dynamics. The
GLE of $\bm{x}$ is derived from the original equation (6) describing the
evolution of $\bm{z}=[\bm{x};\bm{e}]$, while maintaining the effect of the
unresolved part $\bm{e}$. This is particularly useful in our case, as we only
need $\bm{x}$ for infection probability and influence estimation.

Define the Liouville operator $\mathcal{L}$ such that
$\mathcal{L}[g](\bm{z})\mathrel{\mathop{\ordinarycolon}}=\bar{\bm{f}}(\bm{z})\cdot\nabla_{\bm{z}}g(\bm{z})$
for any real-valued function $g$ of $\bm{z}$. Let $e^{t\mathcal{L}}$ be the
Koopman operator associated with $\mathcal{L}$ such that
$e^{t\mathcal{L}}g(\bm{z}(0))=g(\bm{z}(s))$ where $\bm{z}(t)$ solves (6). Then
$\mathcal{L}$ is known to satisfy the semi-group property for all $g$, i.e.,
$e^{t\mathcal{L}}g(\bm{z})=g(e^{t\mathcal{L}}\bm{z})$. Now consider the
projection operator $\mathcal{P}$ as the truncation such that
$(\mathcal{P}g)(\bm{z})=(\mathcal{P}g)(\bm{x},\bm{e})=g(\bm{x},0)$ for any
$\bm{z}=(\bm{x},\bm{e})$, and its orthogonal complement as
$\mathcal{Q}=I-\mathcal{P}$ where $I$ is the identity operator. The following
theorem describes the _exact_ evolution of $\bm{x}(t)$, and the proof is given
in Section A.2.

######  Theorem 2.

The evolution of $\bm{x}$ specified in (6) can also be described by the
following GLE:

$\bm{x}^{\prime}=\bm{f}(\bm{x};\bm{A})+\int_{0}^{t}\bm{k}(t-s,\bm{x}(s))\operatorname{d\\!}s,$
(9)

where $\bm{f}$ is given in (7), and
$\bm{k}(t,\bm{x})\mathrel{\mathop{\ordinarycolon}}=\mathcal{P}\mathcal{L}e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{x}$.

Note that, (9) is _not_ an approximation—it is an _exact_ representation of
the $\bm{x}$ part of the original problem (6). The equation (9) can be
interpreted as a _mean-field_ equation, where the two terms on the right hand
side are called the _streaming term_ (corresponding to the mean-field
dynamics) and _memory term_ , respectively. The mean-field dynamics provide
the _main drift_ of the evolution, and the memory term in a convolution form
is for vital _adjustment_. This inspires us to approximate the memory term as
a time convolution on $\bm{x}$, which naturally yields a delay differential
equation reduced a continuous-time neural network, as shown in the next
subsection.

###  2.3 Memory approximation and delay differential equation

To compute the evolution (9) of $\bm{x}$, we consider an approximation of the
Mori-Zwanzig memory term by a neural network $\bm{\varepsilon}$ with time
convolution of $\bm{x}$ as follows,

$\int_{0}^{t}\bm{k}(t-s,\bm{x}(s))\operatorname{d\\!}s\approx\bm{\varepsilon}(\bm{x}(t),\bm{h}(t);\bm{\eta})\quad\mbox{where}\quad\bm{h}(t)=\int_{0}^{t}\bm{K}(t-s;\bm{w})\bm{x}(s)\operatorname{d\\!}s.$
(10)

In (10), $\bm{K}(\cdot;\bm{w})$ is a convolutional operator with parameter
$\bm{w}$, and $\bm{\varepsilon}(\bm{x},\bm{h};\bm{\eta})$ is a deep neural
network with $(\bm{x},\bm{h})$ as input and $\bm{\eta}$ as parameter. Both
$\bm{w}$ and $\bm{\eta}$ are to be trained by the cascade data $\mathcal{D}$.
Hence, (9) reduces to the _delay differential equation_ which involves a time
integral $\bm{h}(t)$ of past $\bm{x}$:

$\bm{x}^{\prime}=\tilde{\bm{f}}(\bm{x},\bm{h};\bm{\theta})\mathrel{\mathop{\ordinarycolon}}=\bm{f}(\bm{x};\bm{A})+\bm{\varepsilon}(\bm{x},\bm{h};\bm{\eta}).$
(11)

The initial condition of (11) with source set $\mathcal{S}$ is given by

$\bm{x}(0)={\bm{\chi}}_{\mathcal{S}},\quad\bm{h}(0)=\bm{0},\quad\mbox{and}\quad\bm{x}(t)=\bm{h}(t)=\bm{0},\quad\forall\,t<0.$
(12)

We call the system (11) with initial (12) the _neural mean-field_ (NMF)
dynamics.

The delay differential equation (11) is equivalent to a coupled system of
$(\bm{x},\bm{h})$ which is shown in the following theorem, whose proof is
provided in Section A.3.

######  Proposition 2.1.

The delay differential equation (11) is equivalent to the following coupled
system of $(\bm{x},\bm{h})$:

$\displaystyle\bm{x}^{\prime}(t)$
$\displaystyle=\tilde{\bm{f}}(\bm{x}(t),\bm{h}(t);\bm{A},\bm{\eta})=\bm{f}(\bm{x}(t);\bm{A})+\bm{\varepsilon}(\bm{x}(t),\bm{h}(t);\bm{\eta})$
(13a) $\displaystyle\bm{h}^{\prime}(t)$
$\displaystyle=\int_{0}^{t}\bm{K}(t-s;\bm{w})\tilde{\bm{f}}(\bm{x}(s),\bm{h}(s);\bm{A},\bm{\eta})\operatorname{d\\!}s$
(13b)

with initial condition (12). In particular, if
$\bm{K}(t;\bm{w})=\sum_{l=1}^{L}\bm{B}_{l}e^{-\bm{C}_{l}t}$ for some
$L\in\mathbb{N}$ with
$\bm{w}=\\{(\bm{B}_{l},\bm{C}_{l})_{l}\mathrel{\mathop{\ordinarycolon}}\bm{B}_{l}\bm{C}_{l}=\bm{C}_{l}\bm{B}_{l},\,\forall\,l\in[L]\\}$,
then (13) can be solved by a non-delay system of $(\bm{x},\bm{h})$ with (13a)
and $\bm{h}^{\prime}=\sum_{l=1}^{L}(\bm{B}_{l}\bm{x}-\bm{C}_{l}\bm{h})$.

In the remainder of this paper, we only consider the linear kernel
$\bm{K}(t;\bm{w})=\bm{B}e^{-\bm{C}t}$ where $\bm{B}$ and $\bm{C}$ commute for
simplicity. As shown in Proposition 2.1, NMF (11) reduces to a non-deday ODE
system of $(\bm{x},\bm{h})$ with (13a) and
$\bm{h}^{\prime}=\bm{B}\bm{x}-\bm{C}\bm{h}$. Solving such a system for the
optimal parameter $\bm{\theta}=(\bm{A},\bm{w},\bm{\eta})$ has been cast as the
so-called neural ODE (NODE) in [5]. In the following subsection, we establish
a direction connection between mathematical optimal control and NODE, and
provide a rigorous proof that NODE exactly evaluates the gradient of the
target payoff function (likelihood function in our case) during optimization
from the optimal control point of view. Compared to [5], our proof is based on
calculus of variation which is more mathematically rigorous. Moreover, we show
how to incorporate the running payoff (or loss) function at scattered
observation times through a rigorous derivation, as needed in continuous-time
NMF training.

Note that, once the optimal $\bm{\theta}$ is obtained, we can extract $\bm{A}$
for network inference. Moreover, we can compute $\bm{x}(t)$ for all $t$ using
(11) with any given source set $\mathcal{S}$, which solves the influence
estimation problem. Therefore, the network inference and influence estimation
problems can be tackled simultaneously by the parameter training of NMF.

###  2.4 Parameter training and optimal control

To obtain explicit form of NMF (11) for influence estimation and network
inference, we need to know the network parameters
$\bm{\theta}=(\bm{A},\bm{\eta},\bm{w})$. Let $[0,T]$ be the time horizon of
the cascade data
$\mathcal{D}=\\{\mathcal{C}_{k}=(\mathcal{S}_{k},\bm{\tau}_{k})\mathrel{\mathop{\ordinarycolon}}k\in[K]\\}$,
i.e., all cascades are recorded up to time $T$. Given any particular
$\mathcal{C}=(\mathcal{S},\bm{\tau})\in\mathcal{D}$ where
${\bm{\tau}}=\\{t_{i}\in[0,T]\cup\\{\infty\\}\mathrel{\mathop{\ordinarycolon}}i\in[n]\\}$,
it suffices to derive the negative log-likelihood function of the infection
probabilities $\bm{x}$ given (11) with parameter $\bm{\theta}$ for this
cascade $\mathcal{C}$. The total negative log-likelihood function is thus the
sum of such function over all the $K$ cascades in $\mathcal{D}$.

Recall from Section 2.1 that $\bm{X}(t)$ is the jump stochastic process
describing the infection state of the nodes and $\bm{x}(t)$ is the infection
probabilities. Therefore, $\bm{x}^{\prime}(t)$ is essentially the (non-
conditional) intensity of $\bm{X}(t)$. In other words, $\bm{X}(t)$ is
identical to a non-homogeneous Poisson process with intensity function
$\bm{x}^{\prime}(t)$ for $t$ almost everywhere. Due to the relation between
the intensity function and the likelihood function of a point process [23],
the negative log-likelihood function of $\bm{x}^{\prime}(t)$ given the cascade
$\mathcal{C}=(\mathcal{S},\bm{\tau})$ can be easily obtained, and it is also
the “loss function” $\ell$ we need to minimize:

$\ell(\bm{x};\mathcal{C})=\sum_{i=1}^{n}\mathinner{\Bigl{(}-\log
x_{i}^{\prime}(t_{i})+{x}_{i}(T)\Bigr{)}}=\int_{0}^{T}r(\bm{x}(t),\bm{\theta})\operatorname{d\\!}t+\bm{1}^{\top}\bm{x}(T),$
(14)

where the running loss function is defined by

$r(\bm{x}(t),\bm{\theta})=\sum_{i=1}^{n}-\delta(t-t_{i})\log
x_{i}^{\prime}(t)=\sum_{i=1}^{n}-\delta(t-t_{i})\log(\tilde{\bm{f}}(\bm{x},\bm{h};\bm{\theta}))_{i},$
(15)

and $\delta(\cdot)$ is the Dirac delta function. The running loss takes into
account the changes of $\bm{x}(t)$ at intermediate times during $[0,T]$.

We can further add regularization or incorporate prior information to (14). In
particular, if $\mathcal{E}$ is given, we know that $\bm{A}$ must be supported
on $\mathcal{E}$, which serves as the constraint of $\bm{A}$. If we know the
network has low density (sparse edges), then we can enforce a sparsity
regularization such as $\|\bm{A}\|_{1}$ (the $l_{1}$ norm of the vectorized
$\bm{A}\in\mathbb{R}^{n^{2}}$). In general, $\bm{A}$ can be interpreted as the
convolution to be learned from a graph convolution network (GCN)[29, 59]. The
support and magnitude of $\bm{A}$ implies the network structure and strength
of interaction between pairs of nodes, respectively. We will provide more
details of our choice of regularization and its parameter setting in Section
4.

To summarize, the optimal parameter $\bm{\theta}$ of (11) can be obtained by
minimizing the loss function in (14) for the given cascade $\mathcal{C}$:

$\displaystyle\min_{\bm{\theta}}\quad$
$\displaystyle\ell(\bm{\theta};\mathcal{C})\mathrel{\mathop{\ordinarycolon}}=\int_{0}^{T}r(\bm{x}(t),\bm{\theta})\operatorname{d\\!}t+\bm{1}^{\top}\bm{x}(T),$
(16a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\bm{m}^{\prime}(t)=\bm{g}(\bm{m}(t);\bm{\theta}),\quad\bm{m}(0)=[{\bm{\chi}}_{\mathcal{S}_{k}};\bm{0}],\quad
t\in[0,T],$ (16b)

where $\bm{m}(t)=[\bm{x}(t);\bm{h}(t)]\in\mathbb{R}^{2n}$ and

$\bm{g}(\bm{m};\bm{\theta})=\begin{pmatrix}\bm{A}\bm{x}-\operatorname{diag}(\bm{x})\bm{A}\bm{x}+\bm{\varepsilon}(\bm{x},\bm{h};\bm{\eta})\\\
\bm{B}\bm{x}-\bm{C}\bm{h}\end{pmatrix}.$ (17)

This is the parameter training problem given one cascade $\mathcal{C}$, and
can be trivially extended to the case $\mathcal{D}$ which consists of $K$
cascades. In what follows, we drop the notation $\mathcal{C}$ for conciseness.

In (17), $\bm{g}(\bm{m};\bm{\theta})$ is the NMF dynamics derived in (13) with
parameter $\bm{\theta}=(\bm{A},\bm{w},\bm{\eta})$ and
$\bm{w}=(\bm{B},\bm{C})$. In particular, $\bm{\eta}$ stands for the network
parameters of the deep neural network $\bm{\varepsilon}$ that wraps the
augmented state $\bm{m}$ to approximate the MZ memory term (9). It is worth
noting that, the so-called control variable $\bm{\theta}$ is constant and
time-invariant in NODE [5] as well as in NMF. Therefore, it is considerably
easier to handle than that in classical optimal control $\bm{\theta}$.
Specifically, we can develop an algorithm for solving $\bm{\theta}$ in (16)
which is easy to implement. Moreover, we can derive rigorous proof of the
relation between the gradient of the loss function and the solution of the
augmented dynamics.

As we can see, to find the optimal $\bm{\theta}$ of (16), the key is to
compute $\nabla_{\bm{\theta}}\ell(\bm{\theta})$ for any $\bm{\theta}$. To this
end, we recall that the _Hamiltonian_ function associated with the control
problem (16) is

$H(\bm{m}(t),\bm{p}(t);\bm{\theta})=\bm{p}(t)\cdot\bm{g}(\bm{m}(t);\bm{\theta})+r(\bm{m}(t),\bm{\theta}),$
(18)

where $\bm{p}(t)\in\mathbb{R}^{2n}$ is the co-state variable (also known as
the adjoint variable) associated with $\bm{m}(t)$. Here, $\bm{p}(t)$ plays the
role of Lagrange multiplier for the ODE constraint (16b). The standard optimal
control theory states that the co-state $\bm{p}(t)$ follows the ODE backward
in time as follows:

$\begin{cases}\bm{p}^{\prime}(t)=-\nabla_{\bm{m}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{p}(t)-\nabla_{\bm{m}}r(\bm{m}(t),\bm{\theta}),&\quad
T\geq t\geq 0,\\\ \bm{p}(T)=[\bm{1};\bm{0}].\end{cases}$ (19)

The terminal condition $\bm{p}(T)=[\bm{1};\bm{0}]$ has this simple form
because the “terminal loss” in (16) is given by
$[\bm{1};\bm{0}]\cdot\bm{m}(T)=\bm{1}\cdot\bm{x}(T)$.

Now we show that $\nabla_{\bm{\theta}}\ell$ can be obtained by solving the ODE
(16b) forward in time and an augmented ODE backward in time. To this end, we
need the following theorem, whose proof is given in Appendix A.4.

######  Theorem 3.

The gradient $\nabla_{\bm{\theta}}\ell(\bm{\theta})$ of the loss function
$\ell$ defined in (16) for any parameter $\bm{\theta}$ and cascade data
$\mathcal{C}$ is given by

$\nabla_{\bm{\theta}}\ell(\bm{\theta})=\int_{0}^{T}\mathinner{\Bigl{(}\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{p}(t)+\nabla_{\bm{\theta}}r(\bm{m}(t),\bm{\theta})\Bigr{)}}\operatorname{d\\!}t.$
(20)

Moreover, if $\bm{m}^{*}$ is the solution of (16b) using the optimal solution
$\bm{\theta}^{*}$ to (16), and $\bm{p}^{*}$ is the co-state determined by (19)
with $\bm{m}^{*}$ and $\bm{\theta}^{*}$, then
$\int_{0}^{T}\nabla_{\bm{\theta}}H(\bm{m}^{*}(t),\bm{p}^{*}(t);\bm{\theta})\operatorname{d\\!}t=\bm{0}$.

The formula (20) in Theorem 3 suggests that we can compute
$\nabla_{\bm{\theta}}\ell$ by tracking an auxiliary variable $\bm{q}$ that
follows the backward differential equation and terminal condition:

$\begin{cases}\bm{q}^{\prime}(t)=-\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t),\bm{\theta})^{\top}\bm{p}(t)-\nabla_{\bm{\theta}}r(\bm{\theta},\bm{m}(t)),&\quad
T\geq t\geq 0,\\\ \bm{q}(T)=\bm{0}.\end{cases}$ (21)

Then (20) implies that
$\nabla_{\bm{\theta}}\ell(\bm{\theta})=\bm{q}(T)-\int_{0}^{T}\bm{q}^{\prime}(t)\operatorname{d\\!}t=\bm{q}(0)$.

Before closing this section, we need to clarify one implementation issue with
the running loss $r$. Suppose that the infection times in the cascade
$\mathcal{C}$ can be sorted
$0<t^{(1)}<t^{(2)}<\cdots<t^{(m)}<t^{(m+1)}\mathrel{\mathop{\ordinarycolon}}=T$.
That is, there are $m$ infections (excluding the infections at the source
nodes) during the cascade $\mathcal{C}$. (Note that any two infection times
coincide with probability 0 since the point process is simple.) For notation
simplicity, suppose that at time $t^{(i)}$, the new infected node is $i$.

Then the integral of the running loss reduces to

$\displaystyle\int_{0}^{T}\nabla_{\bm{\theta}}r(\bm{\theta},\bm{m})\operatorname{d\\!}t$
$\displaystyle=\sum_{i=0}^{m}\nabla_{\bm{\theta}}\left(-\log\bm{g}_{i}(\bm{m}(t^{(i)});\bm{\theta})\right),$
(22)

where $\bm{g}_{i}(\bm{m}(t),\bm{\theta})$ is the $i$th component of
$\bm{g}(\bm{m}(t),\bm{\theta})$. Hence, we need to compute $\bm{q}(0)$ by
solving the ODE (21) backward in each time interval as

$\displaystyle\bm{q}(t^{(i-1)})$
$\displaystyle=\bm{q}(t^{(i)})-\int_{t^{(i)}}^{t^{(i-1)}}\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t),\bm{\theta})^{\top}\bm{p}(t)\operatorname{d\\!}t-\nabla_{\bm{\theta}}\log\bm{g}_{i-1}(\bm{m}(t^{(i-1)});\bm{\theta}).$
(23)

Similarly, we have

$\displaystyle\bm{p}(t^{(i-1)})$
$\displaystyle=\bm{p}(t^{(i)})-\int_{t^{(i)}}^{t^{(i-1)}}\nabla_{\bm{m}}\bm{g}(\bm{m}(t),\bm{\theta})^{\top}\bm{p}(t)\operatorname{d\\!}t-\nabla_{\bm{m}}\log\bm{g}_{i-1}(\bm{m}(t^{(i-1)});\bm{\theta}).$
(24)

The ODE of $\bm{m}(t)$ remains the same as in (16b) since it does not involve
the running loss $r$.

To summarize, in order to compute $\nabla_{\bm{\theta}}\ell(\bm{\theta})$ for
any given $\bm{\theta}$, we need to first solve the ODE (16b) of $\bm{m}(t)$
forward in time from $0$ to $T$; Then we need to solve the ODE system (16b),
(19), and (21) of $(\bm{m}(t),\bm{p}(t),\bm{q}(t))$ backward in time from $T$
to $0$. In particular, we need to solve the backward ODE such that the last
term (24) and (23) are added for $\bm{p}(t)$ and $\bm{q}(t)$ in each time
interval $(t^{(i-1)},t^{(i)}]$. Finally, we obtain
$\nabla_{\bm{\theta}}\ell(\bm{\theta})=\bm{q}(0)$. The complete training
process is summarized in Algorithm 1, where mini-batches of cascades are used
to compute the stochastic gradient in searching the (local) minimizer
$\bm{\theta}$. We did not include the gradient of the regularization of
$\bm{\theta}$, but its computation is standard and can be easily added to
$\nabla_{\bm{\theta}}\ell(\bm{\theta})$.

Algorithm 1 Neural mean-field (NMF) dynamics

1: Input:
$\mathcal{D}=\\{\mathcal{C}_{k}=(\mathcal{S}_{k},\bm{\tau}_{k})\mathrel{\mathop{\ordinarycolon}}k\in[K]\\}$.

2: Initialization: Network architecture $\bm{g}(\cdot;\bm{\theta})$ and
parameter $\bm{\theta}=(\bm{A},\bm{\eta},\bm{w})$.

3: for $k=1,\dots,\text{MaxIterations}$ do

4: Sample a mini-batch of cascades $\hat{\mathcal{D}}\subset\mathcal{D}$.

5: Compute $\bm{m}(t)$ in (16b) forward in time for each
$\mathcal{C}\in\hat{\mathcal{D}}$. (Forward pass)

6: Compute
$\sum_{\mathcal{C}\in\hat{\mathcal{D}}}\nabla_{\bm{\theta}}\ell(\bm{\theta};\mathcal{C})$
using the BackwardMode below. (Backward pass)

7: Update parameter $\bm{\theta}$ using ADAM with stochastic gradient
$\sum_{\mathcal{C}\in\hat{\mathcal{D}}}\nabla_{\bm{\theta}}\ell(\bm{\theta};\mathcal{C})$.

8: end for

9: Output: Network parameter $\bm{\theta}$.  BackwardMode

10: Input: Cascade $\mathcal{C}=(\mathcal{S},\bm{\tau})$ with
$\bm{\tau}\mathrel{\mathop{\ordinarycolon}}0=t^{(0)}<t^{(1)}<\cdots<t^{(m+1)}=T$
and $\bm{m}(T)$.

11: Terminal augmented state:
$[\bm{m}(T);\bm{p}(T);\bm{q}(T)]=[\bm{m}(T);[\bm{1};\bm{0}];\bm{0}]$.

12: for $i=m+1,\dots,1$ do

13: Solve the ODE below backward in time $(t^{(i-1)},t^{(i)}]$:
$\begin{pmatrix}\bm{m}^{\prime}(t)\\\ \bm{p}^{\prime}(t)\\\
\bm{q}^{\prime}(t)\end{pmatrix}=\begin{pmatrix}\bm{g}(\bm{m}(t);\bm{\theta})\\\
-\nabla_{\bm{m}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{p}(t)\\\
-\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{p}(t)\end{pmatrix}$ with
terminal condition $[\bm{m}(t^{(i)});\bm{p}(t^{(i)});\bm{q}(t^{(i)})]$.

14:
$\bm{p}(t^{(i-1)})\leftarrow\bm{p}(t^{(i-1)})-\nabla_{\bm{m}}\log\bm{g}_{i-1}(\bm{m}(t^{(i-1)});\bm{\theta})$.

15:
$\bm{q}(t^{(i-1)})\leftarrow\bm{q}(t^{(i-1)})-\nabla_{\bm{\theta}}\log\bm{g}_{i-1}(\bm{m}(t^{(i-1)});\bm{\theta})$.

16: end for

17: Output:
$\nabla_{\bm{\theta}}\ell(\bm{\theta};\mathcal{C})\leftarrow\bm{q}(0)$.

##  3 Influence Maximization with Learned NMF

In this section, we show how the proposed NMF can be used to tackle an
important but very challenging problem known as the influence maximization.
Suppose we have trained an NMF with parameters $\bm{\theta}$ in Algorithm 1,
such that we can estimate $\bm{x}(t)$ for any $t\in[0,T]$ and any given source
node set ${\bm{\chi}}_{\mathcal{S}}$. Then the goal of influence maximization
is to identify ${\bm{\chi}}_{\mathcal{S}}\in\\{0,1\\}^{n}$ such that its
influence at the prescribed time $T$ (or any other prescribed $t\in(0,T)$) is
maximized. Namely, our goal is to solve the following optimization problem

$\max_{{\bm{\chi}}_{\mathcal{S}}}\
\sigma(T;{\bm{\chi}}_{\mathcal{S}})\mathrel{\mathop{\ordinarycolon}}=\bm{1}_{n}^{\top}\bm{x}(T;{\bm{\chi}}_{\mathcal{S}}),\quad\mbox{s.t.}\quad{\bm{\chi}}_{\mathcal{S}}\in\\{0,1\\}^{n},\quad\bm{1}_{n}^{\top}{\bm{\chi}}_{\mathcal{S}}=n_{0},$
(25)

where $n_{0}\in\mathbb{N}$ is the given budget size. Note that
$\bm{x}(T;{\bm{\chi}}_{\mathcal{S}})$ is the first $n$ components of
$\bm{m}(T)$ computed by forward NMF dynamics with initial value
$\bm{m}(0)=[{\bm{\chi}}_{\mathcal{S}};\bm{0}]$. However, (25) is an NP-hard
combinatorial optimization problem [18], we propose to relax the binary-valued
decision vector ${\bm{\chi}}_{\mathcal{S}}$ to $\bm{u}\in[0,1]^{n}$ in the
continuous hypercube $[0,1]^{n}$ as

$\min_{\bm{u}\in\mathcal{U}}\
L(\bm{u})\mathrel{\mathop{\ordinarycolon}}=\mathcal{R}(\bm{u})-\bm{1}_{n}^{\top}\bm{x}(T;\bm{u}),\quad\mbox{where}\quad\mathcal{U}\mathrel{\mathop{\ordinarycolon}}=\\{\bm{u}\in[0,1]^{n}\mathrel{\mathop{\ordinarycolon}}\bm{1}_{n}^{\top}\bm{u}=n_{0}\\},$
(26)

and $\mathcal{R}(\bm{u})$ is a regularizer that encourages all components of
$\bm{u}$ to take values close to either $0$ or $1$. In our experiments, we
simply set $\mathcal{R}(\bm{u})=\sum_{i=1}^{n}u_{i}(1-u_{i})$. Then we employ
the projected gradient descent (PGD) method to solve (26):

$\bm{u}_{l+1}=\Pi_{\mathcal{U}}(\bm{u}_{l}-\gamma_{l}\nabla_{\bm{u}}L(\bm{u}_{l}))\mathrel{\mathop{\ordinarycolon}}=\operatorname*{\mathrm{arg\,min}}_{\bm{u}\in\mathcal{U}}\,\|\bm{u}-(\bm{u}_{l}-\gamma_{l}\nabla_{\bm{u}}L(\bm{u}_{l}))\|^{2},$
(27)

where $l$ is the iteration counter of PGD, $\tau_{l}>0$ is the step size, and
$\Pi_{\mathcal{U}}$ denotes the orthogonal projection onto $\mathcal{U}$. If
$\nabla_{\bm{u}}L(\bm{u}_{l})$ is known, then (27) is a standard quadratic
program (QP) and can be solved efficiently by off-the-shelf solvers.
Therefore, the only remaining question is to compute
$\nabla_{\bm{u}}L(\bm{u})$ for any given $\bm{u}$. The following theorem
states that this quantity can be computed very efficiently using the proposed
NMF dynamics. The proof is provided in Appendix A.5.

######  Theorem 4.

Let $[\bm{m}(t);\bm{s}(t)]$ be the solution of the augmented NMF system:

$\begin{pmatrix}\bm{m}^{\prime}(t)\\\
\bm{s}^{\prime}(t)\end{pmatrix}=\begin{pmatrix}\bm{g}(\bm{m}(t);\bm{\theta})\\\
\nabla_{\bm{x}}\bm{g}_{\bm{x}}(\bm{m};\bm{\theta})^{\top}\bm{s}(t)\end{pmatrix}$
(28)

with initial value $[\bm{m}(0);\bm{s}(0)]=[[\bm{u};\bm{0}];\bm{1}]$ forward in
time $[0,T]$, where $\bm{g}_{\bm{x}}$ is the first $n$ components of $\bm{g}$.
Then $\nabla_{\bm{u}}L(\bm{u})=\nabla_{\bm{u}}\mathcal{R}(\bm{u})-\bm{s}(T)$.

Theorem 4 implies that $\nabla_{\bm{u}}L(\bm{u})$ can be easily computed by
solving NMF augmented by an auxiliary variable $\bm{s}(t)$ forward in time
$[0,T]$ as in (28). Note that the computation complexity of (28) is linear in
the network size $n$ and standard numerical ODE integrators can quickly solve
the ODE to high accuracy. We summarize the steps for solving (26) in Algorithm
2. Note that the output $\bm{u}$ may not be binary, and thus we can set the
largest $n_{0}$ components of $\bm{u}$ to $1$ and the rest to $0$ as the final
source set selection.

Algorithm 2 Influence maximization via neural mean-field dynamics (NMF-InfMax)

1: Input: Trained NMF with $\bm{g}(\cdot;\bm{\theta})$ from Algorithm 1,
budget $n_{0}\in\\{1,\dots,n-1\\}$

2: Initialization: $\bm{u}\in\mathcal{U}$.

3: for $l=1,\dots,\text{MaxIterations}$ do

4: Solve $[\bm{m}(T),\bm{s}(T)]$ from (28) forward in time with initial
$[[\bm{u};\bm{0}];\bm{1}]$. (Forward pass)

5: Set $\hat{\bm{u}}\leftarrow\bm{u}-\gamma\nabla_{\bm{u}}L(\bm{u})$ where
$\nabla_{\bm{u}}L(\bm{u})=\bm{s}(T)$.

6: Solve a QP:
$\bm{u}\leftarrow\operatorname*{\mathrm{arg\,min}}_{\bm{u}\in\mathcal{U}}\,\|\bm{u}-\hat{\bm{u}}\|^{2}$.

7: end for

8: Output: Source set selection $\bm{u}$.

##  4 Numerical Experiments

###  4.1 Implementation details

In our NMF implementation, the neural mean field dynamic
$\bm{g}(\cdot;\bm{\theta})$ derived from Proposition 2.1 is learned as
Algorithm 1, where $\bm{\varepsilon}(\bm{x},\bm{h};\bm{\eta})$ is a three-
layer fully connected network. Specifically, the input layer size of
$\bm{\varepsilon}$ is $2n$, and both of the hidden and output layer sizes are
$n$. We use Exponential Linear Unit (ELU) as the activation function. The
output is truncated into $[0,1]$. We use the $\ell_{0}$-norm regularization
approximated by log-sum introduced by [50]. The NMF networks are trained and
tested in PyTorch [49] by Adam optimizer [28] with default parameters
(lr=0.001, $\beta_{1}$=0.9, $\beta_{2}$=0.999, $\epsilon$=1e-8) on a Linux
workstation with Intel i9 8-Core Turbo 5GHz CPU, 64GB of memory, and an Nvidia
RTX 2080Ti GPU. InfluLearner and NetRate are both trained by the Matlab code
published by the original authors. All experiments are performed on the same
machine. Given ground truth node infection probability $\bm{x}^{*}$, the Mean
Absolute Error (MAE) of influence (Inf) and infection probability (Prob) of
the estimated $\bm{x}$ are defined by $|\bm{1}\cdot(\bm{x}(t)-\bm{x}(t)^{*})|$
and $\|\bm{x}(t)-\bm{x}(t)^{*}\|_{1}/n$ for every $t$, respectively. For all
the experiments related on the influence estimation, we also use the scaled
influence MAE $|\bm{1}\cdot(\bm{x}(t)-\bm{x}(t)^{*})|/n$ as an evaluation
metric.

###  4.2 Infection probability and influence function estimation

We first apply the proposed NMF to synthetic diffusion networks where ground
truth node infection probabilities are available for quantitative evaluation.

#### Networks

We use three types of network models [31] to generate these synthetic
networks: hierarchical (Hier) network [8], core-periphery (Core) network [33]
and Random (Rand) network with parameter matrices [0.9,0.1;0.1,0.9],
[0.9,0.5;0.5,0.3], and [0.5,0.5;0.5,0.5], respectively.

For each of these three types of networks, we randomly generate 5 networks of
$(n,d)=(128,4)$ and another 5 networks of $(n,d)=(1024,4)$, where $n$ is the
total number of nodes on the network and $d$ is the average out-degree per
node.

#### Diffusion models and parameters

We simulate the diffusion on these networks such that the infection time are
modeled by exponential distribution (Exp), Rayleigh distribution (Ray), and
general Weibull distribution (Wbl).

Note that our theoretical results in this work are based on diffusion models
using exponential distribution, however, we still conduct experiments on other
distributions to test the performance of NMF empirically. In particular, we
draw the parameters $\alpha_{ji}$ from Unif[0.1,1] to simulate the
heterogeneous interactions between nodes for exponential and Reyleigh
distributions. We generate both of the shape and scale parameters of Weibull
distribution from Unif[1,10] randomly.

#### Training and testing data

We randomly generate 900 source node sets of size varying between 1 and 10,
and simulate 10 diffusion cascades for each source set for training. Thus the
training data consists of $K$=9,000 cascades, all of which are truncated into
time window $[0,T]$ with $T=20$.

We generate 100 additional source sets in a similar way, and then split them
as 50%-validation and 50%-test with ground truth of infection probability and
influence estimated by simulating 10,000 cascades for each source set. This
setting on validation and test data will be used for all the experiments
related to influence estimation and all networks and cascades are generated
using the SNAP package [34].

#### Algorithm and parameter settings

In the training of NMF, the batch size of cascade data is set to 300 and the
number of epochs is 50. The coefficients of the regularization term on
$\bm{A}$ and weight decay in Adam optimizer are set to (0.01,1) and (0.001,0)
for network of size 128 and 1024, respectively. We use Runge-Kutta 4th order
(rk4) method with 40 time steps to solve the ODEs numerically.

#### Comparison algorithm

For comparison, we use InfluLearner [12], which is a state-of-the-art method
that can estimate individual node infection probability directly from cascade
data in the CIC setting as our method.

InfluLearner draws a set of random binary features from certain distribution
for each node $j$ indicating the reachabilities of $j$ by other nodes, and
then uses a convex combination of random basis function to parameterize the
conditional infection probability of the node given a source set over these
binary vectors. To estimate the reachability distribution, InfluLearner
calculates the mean frequency of node $j$ being influenced by a source node
$s$, average over all cascades in the training dataset with the source $s$.

In our test, we set the number of random features to 200 as suggested in [12].

It is worth noting that InfluLearner requires additionally the source identity
for each infection to estimate the coverage functions. That is, InfluLearner
also needs to know the original source node in the source set for each and
every new infection occurred in the cascade in the training data. This
additional information is provided in our simulated data in favor of
InfluLearner. However, it is often unavailable in real-world applications such
as epidemic spreads. The proposed NMF method does not have such restriction.

Moreover, to quantify estimation error, we compute the MAE of node infection
probability and influence at $t_{l}=l$ for $l=1,\dots,20$, and average each
over the 50 test source sets. Since InfluLearner needs to learn the coverage
function for a prescribed time $t$, we have to run it for each of the 20 time
points one by one. In contrast, the proposed NMF is more advantageous since it
can directly estimate the entire evolution of infection probabilities during
$[0,T]$, which is more computationally efficient.

#### Comparison results

We show the numerical results of InfluLearner and NMF for influence estimation
on the three aforementioned synthetic diffusion networks (i.e., Hier, Core,
and Rand) in Figure 2. For each of these three networks, we simulate three
types of diffusion times (i.e., Exp, Ray, and Wbl). Therefore, we have 9
network/diffusion combinations in total. For each of these 9 combinations, we
show the scaled influence MAE (top) and probability MAE (bottom) of
InfluLearner and NMF on networks of size 128 and 1024 as explained above. In
each plot of Figure 2, we show the mean (center line) and standard deviation
(shade) averaged over 5 instances.

As we can observe in Figure 2, the error of NMF is much smaller than that of
InfluLearner for almost all times, except at some early stages and on
Hierarchical network with Weibull distribution.

This demonstrates that NMF is a much more accurate method in influence
estimation.

(a) Core + Exp

(b) Core + Ray

(c) Core + Wbl

(d) Rand + Exp

(e) Rand + Ray

(f) Rand + Wbl

(g) Hier + Exp

(h) Hier + Ray

(i) Hier + Wbl

Figure 2: MAE of scaled influence (top) and node infection probability
(bottom) by InfluLearner [12] and NMF on each of the 9 different combinations
of Core-periphery (Core), Random (Rand) and Hierarchical (Hier) networks, and
exponential (Exp), Rayleigh (Ray) and Weibull (Wbl) diffusion models. Mean
(centerline) and standard deviation (shade) over 50 test source sets are
shown. Each network has two configurations of $(n,d)$: $(128,4)$ and
$(1024,4)$, where $n$ is the number of nodes in the diffusion network, and $d$
is the average out-degree per node.

###  4.3 Network structure inference

In addition to influence estimation, the proposed NMF can also learn the
network structure and the transmission rate matrix $\bm{A}$ as a byproduct
during training. In this test, we examine the quality of the learned $\bm{A}$.
We set the recovered adjacency matrix $\mathcal{E}$ to the binary indicator
matrix $\bm{A}^{\top}\geq\epsilon$. More precisely, once we learned $\bm{A}$
in NMF training, we set the edge $\mathcal{E}$ as $\mathcal{E}_{ij}=1$ if
$\alpha_{ij}=(\bm{A})_{ji}\geq 0.01$ and $0$ otherwise. We set the threshold
$\epsilon=0.01$ because all the transmission rates are between $[0.01,1]$.

#### Evaluation criteria

To evaluate the quality of $\mathcal{E}$ and $\bm{A}$, we use four metrics:
precision (Prc), recall (Rcl), accuracy (Acc), and correlation (Cor), defined
as follows,

$\displaystyle\text{Prc}(\mathcal{E},\mathcal{E}^{*})$
$\displaystyle=\textstyle\frac{|\mathcal{E}\cap\mathcal{E}^{*}|}{|\mathcal{E}^{*}|},\
\
\qquad\qquad\text{Rcl}(\mathcal{E},\mathcal{E}^{*})=\textstyle\frac{|\mathcal{E}\cap\mathcal{E}^{*}|}{|\mathcal{E}|},$
$\displaystyle\text{Acc}(\mathcal{E},\mathcal{E}^{*})$
$\displaystyle=1-\textstyle\frac{|\mathcal{E}-\mathcal{E}^{*}|}{|\mathcal{E}|+|\mathcal{E}^{*}|},\qquad\text{Cor}(A,A^{*})=\textstyle\frac{|\mathrm{tr}(A^{\top}A^{*})|}{\|A\|_{F}\|A^{*}\|_{F}},$

where $|\mathcal{E}|$ counts the number of nonzero entries in $\mathcal{E}$,
and the $\mathcal{E}^{*}$ and $\bm{A}^{*}$ are the ground truths,
respectively. In Cor, $\|A\|_{F}^{2}=\mathrm{tr}(A^{\top}A)$ is the Frobenius
norm of the matrix $A$. Prc is the ratio of edges in $\mathcal{E}^{*}$ that
are recovered in $\mathcal{E}$. Rcl is the ratio of correctly recovered edges
in $\mathcal{E}$. Acc indicates the ratio of the number of common edges shared
by $\mathcal{E}$ and $\mathcal{E}^{*}$ against the total number of edges in
them. Cor measures similarity between $A$ and $A^{*}$ by taking their values
into consideration. All metrics are bounded between $[0,1]$, and higher value
indicates better accuracy.

#### Comparison algorithm

For comparison purpose, we also applied NetRate [16], a state-of-the-art
algorithm that uncovers the network structure and transmission rates from
cascade data. It is worth noting that NetRate requires the knowledge of the
specific diffusion model (e.g., Exp, Ray, or Wbl), so that the likelihood
function can be explicitly expressed. Moreover, NetRate can only estimate
$\bm{A}$ of diffusion networks, but not the influence. In contrast, NMF
tackles both network inference and influence extimation simultaneously. In
terms of computation efficiency, we observed that the implementation of
NetRate provided in [16] runs very slowly for large networks. Therefore, we
only perform comparisons on networks of size $n=128$ in this experiment.

#### Comparison results

We compared the estimated $\mathcal{E}$ and $\bm{A}$ using NetRate and NMF
using the four criteria mentioned above in Table 1 for three types of networks
(Random, Hierarchical, and Core-periphery) and two diffusion models
(Exponential and Rayleigh). In all of these tests, NMF consistently
outperforms NetRate in all accuracy metrics.

Table 1: Performance of network structure inference using NetRate [16] and the proposed NMF on Random, Hierarchical, and Core-periphery networks consisting of 128 nodes and 512 edges with Exponential and Rayleigh as diffusion distribution on edges. Quality of the learned edge set $\mathcal{E}$ and distribution parameter $\bm{A}$ are measured by precision (Prc), recall (Rcl), accuracy (Acc), and correlation (Cor). Larger value indicates higher accuracy. Diffusion | Network | Method | Prc | Rcl | Acc | Cor
---|---|---|---|---|---|---
Exponential | Random |  NetRate | 0.457 | 0.821 | 0.515 | 0.438
NMF | 0.459 | 0.997 | 0.622 | 0.910
Hierarchical |  NetRate | 0.395 | 0.748 | 0.515 | 0.739
NMF | 0.595 | 0.997 | 0.745 | 0.928
Core-periphery |  NetRate | 0.277 | 0.611 | 0.264 | 0.264
NMF | 0.292 | 0.997 | 0.450 | 0.839
Rayleigh | Random |  NetRate | 0.481 | 0.399 | 0.434 | 0.465
NMF | 0.883 | 0.905 | 0.894 | 0.909
Hierarchical |  NetRate | 0.659 | 0.429 | 0.519 | 0.464
NMF | 0.889 | 0.936 | 0.911 | 0.913
Core-periphery |  NetRate | 0.150 | 0.220 | 0.178 | 0.143
NMF | 0.649 | 0.820 | 0.724 | 0.820

We also draw $\bm{A}$ inferred by NetRate and NMF for a visual comparison in
Figure 3. In Figure 3, we show the ground truth $\bm{A}^{*}$ (left), the
matrix $\bm{A}$ inferred by NetRate (middle), and $\bm{A}$ learned by NMF
(right). The values of $\alpha_{ij}$ are indicated by the color—the darker the
red is, the higher the value of $\alpha_{ij}$—and the white pixels represent
where $\alpha_{ij}$ is zero. As we can see, $\bm{A}$ learned by NMF is much
more faithful to $\bm{A}^{*}$ than that by NetRate. This result shows that NMF
is very versatile and robust in learning network structure from cascade data.

(a) True

(b)  NetRate

(c) NMF

Figure 3: Ground truth $\bm{A}^{*}$ (left) and $\bm{A}$ inferred by NetRate
(middle) and NMF (right) in same color scale using cascades from a
Hierarchical network consisting of 128 nodes and 512 edges with exponential
diffusion model. Darker pixel indicates larger value of an entry of $\bm{A}$.

Since NetRate code [16] was implemented in MATLAB and is executed on CPU in
our experiment, the computation times of NetRate and NMF cannot be directly
compared. However, we notice that NetRate takes approximately 10+ hours on
average to infer each network structure $\bm{A}$ in Table 1, whereas NMF only
requires about 300 seconds on average to return both more accurate $\bm{A}$
and an influence estimation mechanism.

(a) Probability MAE

(b) Influence MAE

(c) Train time vs $d$

(d) Train time vs $n$

Figure 4: (a)–(b) MAE of infection probability and influence obtained by
InfluLearner [12] and NMF on Hierarchical networks of size $n=128$ and
increasing $d$ from 4 to 6. (c) Training time (in seconds) of NMF versus
density (average out-degree per node) $d$. (d) Training time (in seconds)
versus network size $n$.

(a) Varying training set size

(b) Influence vs $n_{0}$

Figure 5: (a) Influence generated the source sets selected by NMF-InfMax
trained using increasing number of cascades on Hierarchical networks with
1,024 nodes and 4,096 edges. (b) Influence generated by the source sets
selected by IMINFECTOR and NMF-InfMax on the MemeTracker dataset at $T=10$
hours.

###  4.4 Scalability to network size and density

In this test, we will demonstrate the robustness of NMF in influence
estimation when the network size $n$ and density $d$ vary. Recall that $d$
stands for the average out-degree per node. The larger and/or denser the
network is, the more challenging the estimation becomes. In all the
experiments, we use training data consisting of 9,000 cascades generated from
Hierarchical network and exponential diffusion model and set the batch size to
300.

#### Network size

Recall that we have showed in Figure 2 that NMF consistently outperforms than
InfluLearner when the network size is set to 128 and 1024. To show the
scability of NMF, we further test NMF on increasing network size $n$ from 128
to 2048 (with density $d=4$). To test the training time of NMF, we terminate
the computation when the average MAE of infection probability on validation
data over 20 timepoints $t_{\ell}=\ell$ ($\ell=1,2,\dots,20$) is below 0.07.
Euler method with 40 steps is employed as the ODE solver and the learning rate
of the Adam optimizer is set to 0.0001 for network with 2048 nodes. The
training time of NMF is shown in Figure 4(d), which demonstrate that NMF is
scalable for large network size $n$.

#### Network density

We also test the performance of NMF for varying network density $d$. We
compare the infection probability and influence MAE of InfluLearner and NMF
for varying edge density $d$ set to 4, 5, and 6 on a Hierarchical network on
exponential diffusion model with 128 nodes. Figure 4(a) and Figure 4(b) show
that the MAE of infection probability and influence estimation obtained by
InfluLearner and NMF. These two plots show that NMF is very robust when the
density of the network increases by consistently generating estimates of low
MAE. Figure 4(c) shows the training time of NMF versus network density $d$
while $n=128$ is fixed. In this plot, the computation time is recorded when
the training MAE at time $t_{h}$ is below 0.04, where $t_{h}$ is the time when
on average half of the nodes on the network are infected as indicated by the
ground truth. Here, rk4 method with 40 steps is employed as the ODE solver.
Similarly, Figure 4(d) shows the training time versus network size $n$ while
$d=4$ is fixed. From Figures 4(c) and 4(d), we can see that the computational
cost of NMF grows approximately quadratic in density $d$ and linear in size
$n$.

###  4.5 Influence maximization

This part of the experiment is dedicated to performance evaluation in
influence maximization. Specifically, we use the trained NMF to find the
optimal source set with limited budget for maximal influence by following
Algorithm 2 which is referred to as NMF-InfMax.

#### Comparison algorithms

For comparison purpose, we also test the following methods for influence
maximization.

  * •

IMINFECTOR [46]: IMINFECTOR represents the cascade data into two datasets
consisting of seed-cascade length pairs and seed-influenced node pairs to
approximate the influence spread and infection probability of each node by a
regression model and a probability classifier, respectively. The outputs are
used to reduce the number of candidate seeds and reformulate the computation
of the influence spread in a greedy solution to influence maximization. Like
our method, IMINFECTOR only uses cascade data as inputs, with embedding size
50 and sampling percentage 120, trained for 50 epochs with a learning rate of
0.1. The reduction percentage $P$ is set to 100 to keep full information of
cascades.

  * •

IMM [54]: IMM is a reverse reachable (RR) sketch based method which applies
the standard greedy algorithm for maximum coverage to derive a budget size
node set that covers a large number of RR sets sampled from the given network.
We consider the case when IMM return $(1-1/e-\varepsilon)$-approximate
solution with $\epsilon=0.1$ and parameter $\ell=1$, following the experiments
in [54].

  * •

InfluMax[19, 21]: InfluMax speed up the greedy influence maximization
algorithm by exploiting submodularity. We incorporate it with the influence
estimation algorithm ConTinEst[13]. For ContinEst, we draw 10,000 random
samples, each of which has 5 random labels for each node.

Since IMM and InfluMax both require the knowledge of the transmission matrix
$\bm{A}$, we apply NetRate method to learn $\bm{A}$ from cascade data first,
then feed $\bm{A}$ to these two methods. We remark that NetRate and IMINFECTOR
are both in favor of training data consisting of cascades with the source sets
of size 1 (i.e., only one source node). In contrary, NMF-InfMax does not have
this restriction and thus is more flexible. However, for comparison purpose,
we only feed cascade data with single source node to all methods in this
experiment.

#### Experiment setting

We again use three types of Kronecker graph models: Hierarchical (Hier), Core-
periphery (Core) and Random (Rand) networks, and simulate the diffusion
processes using exponential distribution with transmission rates randomly
sampled from Unif[0.1,1]. For each type of network model, we generate two
networks of size $n=1,024$ with $d=2$ and $d=4$, respectively. We sample 100
source nodes, and for each source node we simulate 10 cascades. Hence we have
a total of $K$=1,000 cascades for training. To train NMF, we set the batch
size to 300 and the number of epochs to 50. The coefficients of the
regularization term on $\bm{A}$ is set to 0.001, and the rk4 method with 40
time steps is employed as the ODE solver. To train NMF-InfMax, we set the step
size to constant 0.01, and terminate PGD if either the iteration number
reaches 500 or the computed influence does not change for 10 consecutive
iterations. In Figure 5(a), we show the accuracy of NMF-InfMax when the number
of training cascades increases from 1,000 to 5,000 for each fixed source set
of size from 1 to 10. As we can observe, the accuracy increases significantly
when the number of cascades grows from 1,000 to 2,000 but then improvements
become insignificant. This suggests that 2,000 cascades is necessary to obtain
more accurate influence maximization results for NMF-InfMax. However, due to
the limited scalibilty of NetRate which performs extremely slowly when the
number of cascades is over 1,000 and average out-degree is over $4$. We also
tested IMINFECTOR with larger training data set, but unlike our method, the
accuracy of IMINFECTOR does not improve over 1,000. Hence we still only feed
1,000 cascades to all the compared methods despite that this choice is only in
favor of three existing methods.

It is also important to note that both InfluMax and IMM require the knowledge
of diffusion model given by the shape and scale parameters of edges for the
computation of NetRate and their own. Thus, they are more vulnerable to model
mis-specification. In this test, we assume they know the ground truth
diffusion model, so they can attain their highest accuracy. However, it is
worth noting that the network inference by NetRate can be very expensive
computationally. For example, it took NetRate up to 160 hours to infer the
network structure from 1,000 cascades of a core-periphery network of $n=1,024$
and $d=4$ in Figure 6(c). In contrast, the computational cost of IMINFECTOR is
very low, but IMINFECTOR is more restrictive on data because it requires that
the training cascades contain the nodes to be selected. This may not be
feasible in practice. Moreover, the influence maximization results obtained by
IMINFECTOR also appear to be worse than that by NMF, as shown below.

#### Comparison results

The influence maximization results obtained by the aforementioned algorithms
and NMF are shown in Figure 6. As we can see, NMF-InfMax consistently returns
more influential source sets with all budget $n_{0}$ for all varying network
structure, density, and budget.

(a) Hierarchical

(b) Random

(c) Core-periphery

Figure 6: Influence of the source sets selected by the compared methods on
three different types of networks: (a) Hierarchical, (b) Random, and (c) Core-
periphery, with exponential diffusion model at $T=10$ and varying source sizes
$n_{0}$ from 1 to 10. Each network consists of 1024 nodes and 2048 edges (top)
or 4096 edges(bottoms).

#### Real data

We extract diffusion cascades from the MemeTracker dataset [30] which includes
300 million blog posts and articles collected from 5,000 active media sites
between March 2011 and February 2012. Following [12], we select the group of
cascades with the keyword ”apple and jobs” and then split them as 60%-train
and 40%-validation for the influence maximization models. As the diffusion
model of real-world cascade data is unknown, we only test IMINFECTOR and NMF-
InfMax. We follow the setting in [12] to compute the influence of any selected
source set: we uniformly sample one cascade from the data for each node in the
set and take the union of all sampled cascades as the set of infected nodes.
We repeat this process for 1,000 times and take the average as the true
influence of the selected set. Figure 5(b) shows the result of influence
maximization results. In Figure 5(b), we set $T=10$ and the source size
$n_{0}=10,20,\cdots,60$, and plot the influence of the source sets selected by
IMINFECTOR and NMF-InfMax. As we can see, NMF-InfMax consistently selected
more influential combination of nodes that generate greater influence than
those selected by IMINFECTOR do.

##  5 Related Work

In this section, we conduct a comprehensive review of the literature related
to the present work. There are several topics involved in the proposed method,
namely, influence estimation, network inference, and influence maximization,
which just emerged within the past decade. These topics are considered
independently, and the methods developed are mostly heuristic or sample-
demanding. In what follows, we discuss these topics and their related work in
order.

###  5.1 Influence estimation

Sampling-based influence estimation methods have been considered for discrete-
time and continuous-time diffusion models. Discrete-time models assume node
infections only occur at discrete time points. Under this setting, the
independent cascade (IC) and linear threshold (LT) models are considered and
the propagation spread of a source set $S$ is simply estimated by the expected
reachable set size of $S$ taken over the randomness of the influence
propagation process in [26]. To improve the efficiency of Monte Carlo
simulations used in influence estimation, a method with provable performance
guarantee is developed which iterates over a sequence of guesses on the true
influence until the verifier accepts in [39]. In [39], the verifier estimates
the influence on multiple sampled graphs using a standard Riemann sum of the
influence function, and accepts if this value is close to the guesses. In [3],
the reverse reachable (RR) sets of nodes are adopted which proved the expected
spread equals $n$ times the fraction of sampled RR sets covered by the source
set. The sample size is controlled by a given threshold [3], a pre-calculated
parameter [55], or some stop conditions [54] to achieve a balance between
efficiency and accuracy. Instead of using the full network structure as the
methods above, sketch-based approaches only characterize propagation instances
for influence computation, such as the method in [10], which considers per-
node summary structures defined by the bottom-$k$ min-bash [9] sketch of the
combined reachability set. In contrast to discrete-time models, continuous-
time diffusion models allow arbitrary event occurrence times and hence are
more accurate in modeling real-world diffusion processes. In continuous-time
independent cascade (CIC) models, influence estimation can be reformulated as
the problem of finding the least label list which contains information about
the distance to the smallest reachable labels from the source [13, 21].
Compared to methods using a fixed number of samples, a more scalable
approximation scheme with a built-in block is developed to minimize the number
of samples needed for the desired accuracy [45]. Inspired by [54], algorithms
proposed in [3, 54, 55] can be extended from the IC model to other discrete-
time models and CIC models by generalizing the definition of RR sets. In [25],
a neural mean-field dynamics approach is proposed, which employs the Mori-
Zwanzig (MZ) formalism to derive the node infection probabilities in discrete-
time setting. The influence function can also be approximated by solving a
jump stochastic differential equation [60] or a deterministic differential
equation that governs the evolution of the influence counter [7].

The aforementioned methods require knowledge of cascade traces [10] or the
diffusion networks, such as node connectivity and node-to-node infection
rates, as well as various assumptions on the diffusion of interests. However,
such knowledge about the diffusion networks may not be available in practice,
and the assumptions on the propagation or data formation are often
application-specific and do not hold in most other problems. InfluLearner [12]
is a state-of-the-art method that does not require knowledge of the underlying
diffusion network. InfluLearner estimates the influence directly from cascades
data in the CIC models by learning the influence function with a
parameterization of the coverage functions using random basis functions.
However, the estimation of random basis function suggested by [12] requires
knowledge of the original source node for every infection, which can be
difficult or impossible to be tracked in real-world applications, such as
epidemic spreads.

In recent years, deep learning techniques have been employed to improve the
scalability of influence estimation on large networks. In particular,
convolutional neural networks (CNNs) and attention mechanism are incorporated
with both network structures and user specific features to learn users’ latent
feature representation in [51]. By piping represented cascade graphs through a
gated recurrent unit (GRU), the future incremental influence of a cascade can
be predicted [36]. RNNs and CNNs are also applied to capture the temporal
relationships on the user-generated contents networks (e.g., views, likes,
comments, reposts) and extract more powerful features in [61]. In methods
based on graph structures, graph neural networks (GNNs) and graph convolution
networks (GCNs) are widely applied. In particular, two coupled GNNs are used
to capture the interplay between node activation states and the influence
spread [4], while GCNs integrated with teleport probability from the domain of
page rank in [35] enhanced the performance of method in [51]. However, these
methods depend critically on the structure or content features of cascades
which is not available in many real-world applications.

###  5.2 Network structure inference

Inference of diffusion network structure is an important problem closely
related to influence estimation. In particular, if the network structure and
infections rates are unknown, one often needs to first infer such information
from a training dataset of sampled cascades, each of which tracks a series of
infection times and locations on the network. Existing methods have been
proposed to infer network connectivity [17, 20, 38, 14] and also the infection
rates between nodes [43, 16, 18]. Submodular optimization is applied to infer
network connectivity [17, 20, 38] by considering the most probable [17] or all
[20, 38] directed trees supported by each cascade. One of the early works that
incorporate spatio-temporal factors into network inference is introduced in
[38]. Utilizing convex optimization, transmission functions [14], the prior
probability [43], and the transmission rate [16] over edges are inferred from
cascades. In addition to static networks, the infection rates are considered
but also in the unobserved dynamic network changing over time [18]. Besides
cascades, other features of dynamical processes on networks have been used to
infer the diffusion network structures. To avoid using predefined transmission
models, the statistical difference of the infection time intervals between
nodes in the same cascade versus those not in any cascade was considered in
[52]. A given time series of the epidemic prevalence, i.e., the average
fraction of infected nodes was applied to discover the underlying network. The
recurrent cascading behavior is also explained by integrating a feature vector
describing the additional features [57]. A graph signal processing (GSP)
approach is developed to infer graph structure from dynamics on networks [41,
11].

###  5.3 Influence maximization

Influence maximization is an important but very challenging problem in real-
world applications of diffusion networks, such as commercial advertising and
epidemic controls. Influence maximization is shown to be an NP-hard problem
under most of diffusion models [37] (e.g., LT, IC, CIC). It was first
formulated in [26] as a combinatorial optimization problem. Under certain
assumptions, the influence function $\sigma(\cdot)$ is a non-negative monotone
submodular function, and a standard greedy method [26, 21] can be applied to
obtain provable sub-optimal solution. Specifically, the greedy method starts
from an empty set $\mathcal{S}$ and gradually add one node $i$ that maximizes
the marginal gain $\sigma(\mathcal{S}\cup\\{i\\})-\sigma(\mathcal{S})$ to
$\mathcal{S}$. Note that this requires repeatedly evaluation of influences
$\sigma(\mathcal{S})$ which affects the result of influence maximization
significantly.

Instead of searching all the nodes in the greedy iterations, a GCN is trained
by a probabilistic greedy mechanism, such that it selects a node with
probability proportional to its marginal gain to identify noise and predict
the node quality for the propagation spread [40]. The computations on the
reward for adding the node to set $\mathcal{S}$ is performed in another
Q-learning network. The importance of nodes can also be measured by exploiting
the submodularity [22, 32] or only considering one-hop and two-hoop spread
benefit measures on nodes in [24]. The influence maximization problem is also
modeled as the maximum coverage problem of selecting the budget number of
nodes to cover the maximum number of sampled RR sets in [3, 54, 55]. For
instances without the information of network structure, the influence
relationships between nodes are representation learned from cascade date
initiated by a single node to derive a greedy solution in [46, 47].

##  6 Conclusion

We propose a novel framework using neural mean-field dynamics for inference
and estimation on diffusion networks. Our new framework is derived from the
Mori-Zwanzig formalism to obtain exact evolution of node infection
probabilities. The Mori-Zwanzig memory can be approximated by convolutions,
which renders the system as a delay differential equation for highly
interpretable parameterization. Directly using information diffusion cascade
data, our framework outperforms many state-of-the-art methods in network
structure inference and influence estimation. Our framework can also
effectively tackle influence maximization on networks, which is known to be a
challenging NP-hard problem. Extensive numerical experiments were conducted to
show the promising accuracy and efficiency of the proposed framework on both
synthetic and real-world data sets. We expect that the proposed framework can
be applied to many other optimization and control problems arising from
diffusion network applications, such as optimal campaigning, propagation
control, and source identification, which will also be investigated in our
future work.

##  Appendix A Proofs

###  A.1 Proof of Theorem 1

###### Proof.

Let $\lambda_{i}^{*}(t)$ be the conditional intensity of node $i$ at time $t$,
i.e., $\mathbb{E}[\operatorname{d\\!}X_{i}(t)|$
$\mathcal{H}(t)]=\lambda_{i}^{*}(t)\operatorname{d\\!}t$. In the standard
diffusion model, the conditional intensity $\lambda_{i}^{*}(t)$ of a healthy
node $i$ (i.e., $X_{i}(t)=0$) is determined by the total infection rate of its
infected neighbors $j$ (i.e., $X_{j}(t)=1$). That is,

$\lambda_{i}^{*}(t)=\sum_{j}\alpha_{ji}X_{j}(t)(1-X_{i}(t)).$ (29)

By taking expectation $\mathbb{E}_{\mathcal{H}(t)}[\cdot]$ on both sides of
(29), we obtain

$\displaystyle\lambda_{i}(t)\mathrel{\mathop{\ordinarycolon}}=\ $
$\displaystyle\mathbb{E}_{\mathcal{H}(t)}[\lambda_{i}^{*}(t)]=\mathbb{E}_{\mathcal{H}(t)}\mathinner{\Bigl{[}\alpha_{ji}X_{j}(t)(1-X_{i}(t))\big{|}\mathcal{H}(t)\Bigr{]}}$
$\displaystyle=\ $
$\displaystyle\sum_{j}\alpha_{ji}(x_{j}-x_{ij})=\sum_{j}\alpha_{ji}(x_{j}-y_{ij}-e_{ij}).$
(30)

On the other hand, there is

$\lambda_{i}(t)\operatorname{d\\!}t=\mathbb{E}_{\mathcal{H}(t)}[\lambda_{i}^{*}(t)]\operatorname{d\\!}t=\mathbb{E}_{\mathcal{H}(t)}[\operatorname{d\\!}X_{i}(t)|\mathcal{H}(t)]=\operatorname{d\\!}\mathbb{E}_{\mathcal{H}(t)}[X_{i}(t)|\mathcal{H}(t)]=\operatorname{d\\!}x_{i}.$
(31)

Combining (30) and (31) yields

$\displaystyle x_{i}^{\prime}$
$\displaystyle=\frac{\operatorname{d\\!}x_{i}(t)}{\operatorname{d\\!}t}=\sum_{j}\alpha_{ji}(x_{j}-y_{ij}-e_{ij})=(\bm{A}\bm{x})_{i}-(\operatorname{diag}(\bm{x})\bm{A}\bm{x})_{i}-\sum_{j}\alpha_{ji}e_{ij}$

for every $i\in[n]$, which verifies the $\bm{x}$ part of (6). Similarly, we
can obtain

$\displaystyle x_{I}^{\prime}$ $\displaystyle=\sum_{i\in I}\sum_{j\notin
I}\alpha_{ji}(x_{I}-x_{I\cup\\{j\\}})=\sum_{i\in I}\sum_{j\notin
I}\alpha_{ji}(y_{I}+e_{I}-y_{I\cup\\{j\\}}-e_{I\cup\\{j\\}}).$ (32)

Moreover, by taking derivative on both sides of $x_{I}(t)=y_{I}(t)+e_{I}(t)$,
we obtain

$\displaystyle x_{I}^{\prime}=\sum_{i\in
I}y_{I\setminus\\{i\\}}x_{i}^{\prime}+e_{I}^{\prime}=\sum_{i\in
I}y_{I\setminus\\{i\\}}\sum_{j\neq
i}\alpha_{ji}(x_{j}-x_{i}x_{j}-e_{ij})+e_{I}^{\prime}.$ (33)

Combining (32) and (33) yields the $\bm{e}$ part of (6).

It is clear that $\bm{x}_{0}={\bm{\chi}}_{\mathcal{S}}$. For every $I$, at
time $t=0$, there is $x_{I}(0)=\prod_{i\in I}X_{i}(0)=1$ if
$I\subset\mathcal{S}$ and $0$ otherwise; and the same for $y_{I}(0)$. Hence
$e_{I}(0)=x_{I}(0)-y_{I}(0)=0$ for all $I$. Hence
$\bm{z}_{0}=[\bm{x}_{0};\bm{e}_{0}]=[{\bm{\chi}}_{\mathcal{S}};\bm{0}]$, which
verifies the initial condition of (6). ∎

###  A.2 Proof of Theorem 2

###### Proof.

Consider the system (6) over a finite time horizon $[0,T]$, which evolves on a
smooth manifold $r\subset\mathbb{R}^{N}$. For any real-valued phase
(observable) space function
$g\mathrel{\mathop{\ordinarycolon}}r\to\mathbb{R}$, the nonlinear system (6)
is equivalent to the linear partial differential equation, known as the
Liouville equation:

$\begin{cases}\partial_{t}u(t,\bm{z})=\mathcal{L}[u](t,\bm{z}),\\\
u(0,\bm{z})=g(\bm{z}),\end{cases}$ (34)

where the Liouville operator
$\mathcal{L}[u]\mathrel{\mathop{\ordinarycolon}}=\bar{\bm{f}}(\bm{z})\cdot\nabla_{\bm{z}}u$.
The equivalency is in the sense that the solution of (34) satisfies
$u(t,\bm{z}_{0})=g(\bm{z}(t;\bm{z}_{0}))$, where $\bm{z}(t;\bm{z}_{0})$ is the
solution to (6) with initial value $\bm{z}_{0}$.

Denote $e^{t\mathcal{L}}$ the Koopman operator associated with $\mathcal{L}$
such that $e^{t\mathcal{L}}g(\bm{z}_{0})=g(\bm{z}(t))$ where $\bm{z}(t)$ is
the solution of (6). Then $e^{t\mathcal{L}}$ satisfies the semi-group
property, i.e.,

$e^{t\mathcal{L}}g(z)=g(e^{t\mathcal{L}}z)$ (35)

for all $g$. On the right hand side of (35), $\bm{z}$ can be interpreted as
$\bm{z}=\bm{\iota}(\bm{z})=[\iota_{1}(\bm{z}),\dots,\iota_{N}(\bm{z})]$ where
$\iota_{j}(\bm{z})=z_{j}$ for all $j$.

Now consider the projection operator $\mathcal{P}$ as the truncation such that
$\mathcal{P}g(\bm{z})=\mathcal{P}g(\bm{x},\bm{e})=g(\bm{x},0)$ for any
$\bm{z}=(\bm{x},\bm{e})$, and its orthogonal complement as
$\mathcal{Q}=I-\mathcal{P}$ where $I$ is the identity operator. Note that
$\bm{z}^{\prime}(t)=\frac{\operatorname{d\\!}\bm{z}(t)}{\operatorname{d\\!}t}=\frac{\partial}{\partial
t}e^{t\mathcal{L}}\bm{z}_{0}$, and
$\bar{\bm{f}}(\bm{z}(t))=e^{t\mathcal{L}}\bm{f}(\bm{z}_{0})=e^{t\mathcal{L}}\mathcal{L}\bm{z}_{0}$
since $\mathcal{L}\iota_{j}(\bm{z})=\bm{f}_{j}(\bm{z})$ for all $\bm{z}$ and
$j$. Therefore (6) implies that

$\frac{\partial}{\partial
t}e^{t\mathcal{L}}\bm{z}_{0}=e^{t\mathcal{L}}\mathcal{L}\bm{z}_{0}=e^{t\mathcal{L}}\mathcal{P}\mathcal{L}\bm{z}_{0}+e^{t\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}.$
(36)

Note that the first term on the right hand side of (36) is

$e^{t\mathcal{L}}\mathcal{P}\mathcal{L}\bm{z}_{0}=\mathcal{P}\mathcal{L}e^{t\mathcal{L}}\bm{z}_{0}=\mathcal{P}\mathcal{L}\bm{z}(t).$
(37)

For the second term in (36), we recall that the well-known Dyson’s identity
for the Koopman operator $\mathcal{L}$ is given by

$e^{t\mathcal{L}}=e^{t\mathcal{Q}\mathcal{L}}+\int_{0}^{t}e^{s\mathcal{L}}\mathcal{P}\mathcal{L}e^{(t-s)\mathcal{Q}\mathcal{L}}\operatorname{d\\!}s.$
(38)

Applying (38) to $\mathcal{Q}\mathcal{L}\bm{z}_{0}$ yields

$\displaystyle e^{t\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}$
$\displaystyle=e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}+\int_{0}^{t}e^{s\mathcal{L}}\mathcal{P}\mathcal{L}e^{(t-s)\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}\operatorname{d\\!}s$
$\displaystyle=e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}+\int_{0}^{t}\mathcal{P}\mathcal{L}e^{(t-s)\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}e^{s\mathcal{L}}\bm{z}_{0}\operatorname{d\\!}s$
(39)
$\displaystyle=e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}+\int_{0}^{t}\mathcal{P}\mathcal{L}e^{(t-s)\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}(s)\operatorname{d\\!}s.$

Substituting (37) and (39) into (36), we obtain

$\frac{\partial}{\partial
t}e^{t\mathcal{L}}\bm{z}_{0}=\mathcal{P}\mathcal{L}\bm{z}(t)+e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}_{0}+\int_{0}^{t}\mathcal{P}\mathcal{L}e^{(t-s)\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}(s)\operatorname{d\\!}s,$
(40)

where we used the fact that
$e^{t\mathcal{L}}\mathcal{P}\mathcal{L}\bm{z}_{0}=\mathcal{P}\mathcal{L}e^{t\mathcal{L}}\bm{z}_{0}=\mathcal{P}\mathcal{L}\bm{z}(t)$.
Denote
$\bm{\phi}(t,\bm{z})\mathrel{\mathop{\ordinarycolon}}=e^{t\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}$,
then we simplify (40) into

$\frac{\partial}{\partial
t}e^{t\mathcal{L}}\bm{z}_{0}=\mathcal{P}\mathcal{L}\bm{z}(t)+\bm{\phi}(t,\bm{z}_{0})+\int_{0}^{t}\bm{k}(t-s,\bm{z}(s))\operatorname{d\\!}s,$
(41)

where
$\bm{k}(t,\bm{z})\mathrel{\mathop{\ordinarycolon}}=\mathcal{P}\mathcal{L}\bm{\phi}(t,\bm{z})=\mathcal{P}\mathcal{L}e^{t\mathcal{L}}\mathcal{Q}\mathcal{L}\bm{z}$.

Now consider the evolution of $\bm{\phi}(t,\bm{z})$, which is given by

$\partial_{t}\bm{\phi}(t,\bm{z}_{0})=\mathcal{Q}\mathcal{L}\bm{\phi}(t,\bm{z}_{0}),$
(42)

with initial condition
$\bm{\phi}(0,\bm{z}_{0})=\mathcal{Q}\mathcal{L}\bm{z}_{0}=\mathcal{L}\bm{z}_{0}-\mathcal{P}\mathcal{L}\bm{z}_{0}=\bar{\bm{f}}(\bm{x}_{0},\bm{e}_{0})-\bar{\bm{f}}(\bm{x}_{0},\bm{0})=\bm{0}$
since $\bm{e}_{0}=\bm{0}$. Applying $\mathcal{P}$ on both sides of (42) yields

$\partial_{t}\mathcal{P}\bm{\phi}(t,\bm{z}_{0})=\mathcal{P}\mathcal{Q}\mathcal{L}\bm{\phi}(t,\bm{z}_{0})=\bm{0},$

with initial $\mathcal{P}\bm{\phi}(0,\bm{z}_{0})=\bm{0}$. This implies that
$\mathcal{P}\bm{\phi}(t,\bm{z}_{0})=\bm{0}$ for all $t$. Hence, applying
$\mathcal{P}$ to both sides of (40) yields

$\frac{\partial}{\partial t}\mathcal{P}\bm{z}(t)=\frac{\partial}{\partial
t}\mathcal{P}e^{t\mathcal{L}}\bm{z}_{0}=\mathcal{P}\mathcal{L}\bm{z}(t)+\int_{0}^{t}\mathcal{P}\bm{k}(t-s,\bm{z}(s))\operatorname{d\\!}s.$
(43)

Restricting to the first $n$ components, $\mathcal{P}\bm{z}(t)$ reduces to
$\bm{x}(t)$ and $\mathcal{P}\bm{k}(t-s,\bm{z}(s))$ reduces to
$\bm{k}(t-s,\bm{x}(s))$. Recalling that
$\mathcal{P}\mathcal{L}\bm{z}(t)=\mathcal{P}\bar{\bm{f}}(\bm{z}(t))=\bar{\bm{f}}(\bm{x}(t),\bm{0})=\bm{f}(\bm{x}(t))$
completes the proof. ∎

###  A.3 Proof of Proposition 2.1

###### Proof.

From the definition of $\bm{h}(t)$ in (44), we obtain

$\bm{h}(t)=\int_{0}^{t}\bm{K}(t-s;\bm{w})\bm{x}(s)\operatorname{d\\!}s=\int_{-\infty}^{t}\bm{K}(t-s;\bm{w})\bm{x}(s)\operatorname{d\\!}s=\int_{0}^{\infty}\bm{K}(s;\bm{w})\bm{x}(t-s)\operatorname{d\\!}s$
(44)

where we used the fact that $\bm{x}(t)=0$ for $t<0$. Taking derivative on both
sides of (44) yields

$\displaystyle\bm{h}^{\prime}(t)$
$\displaystyle=\int_{0}^{\infty}\bm{K}(s;\bm{w})\bm{x}^{\prime}(t-s)\operatorname{d\\!}s=\int_{0}^{\infty}\bm{K}(s;\bm{w})\tilde{\bm{f}}(\bm{x}(t-s),\bm{h}(t-s);\bm{A},\bm{\eta})\operatorname{d\\!}s$
$\displaystyle=\int_{-\infty}^{t}\bm{K}(t-s;\bm{w})\tilde{\bm{f}}(\bm{x}(s),\bm{h}(s);\bm{A},\bm{\eta})\operatorname{d\\!}s=\int_{0}^{t}\bm{K}(t-s;\bm{w})\tilde{\bm{f}}(\bm{x}(s),\bm{h}(s);\bm{A},\bm{\eta})\operatorname{d\\!}s$

where we used the fact that
$\bm{x}^{\prime}(t)=\tilde{\bm{f}}(\bm{x}(t),\bm{h}(t);\bm{A},\bm{\eta})=0$
for $t<0$ in the last equality.

If $\bm{K}(t;\bm{w})=\sum_{l}\bm{B}_{l}e^{-\bm{C}_{l}t}$, then we can take
derivative of (44) and readily deduce that
$\bm{h}^{\prime}=\sum_{l=1}^{L}(\bm{B}_{l}\bm{x}-\bm{C}_{l}\bm{h})$. ∎

###  A.4 Proof of Theorem 3

###### Proof.

Let $\bm{\zeta}\in\mathbb{R}^{m}$ and $\varepsilon\geq 0$ be arbitrary.
Consider the variation of any control $\bm{\theta}$ given by
$\bm{\theta}_{\varepsilon}\mathrel{\mathop{\ordinarycolon}}=\bm{\theta}+\varepsilon\bm{\zeta}$
and denote $\bm{m}_{\varepsilon}(t)$ the state process following (16b) with
$\bm{\theta}_{\varepsilon}$. Then we have

$\bm{m}_{\varepsilon}(t)=\bm{m}(t)+\varepsilon\bm{y}(t)+o(\varepsilon),\qquad
0\leq t\leq T,$

where the first-order perturbation $\bm{y}(t)$ satisfies

$\begin{cases}\bm{y}^{\prime}(t)=\nabla_{\bm{m}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{y}(t)+\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{\zeta},\qquad
0\leq t\leq T,\\\ \bm{y}(0)=\bm{0}.\end{cases}$

Therefore, the directional derivative of $\ell$ defined in (16a) at
$\bm{\theta}$ along the direction $\bm{\zeta}$ is

$\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}\varepsilon}\ell(\bm{\theta}_{\varepsilon})\Big{|}_{\varepsilon=0}$
$\displaystyle=\int_{0}^{T}\mathinner{\Bigl{(}\nabla_{\bm{m}}r(\bm{m}_{t},\bm{\theta})\bm{y}(t)+\nabla_{\bm{\theta}}r(\bm{m}(t),\bm{\theta})\bm{\zeta}\Bigr{)}}\operatorname{d\\!}t+\bm{p}(T)\bm{y}(T).$
(45)

On the other hand, we have

$\displaystyle(\bm{p}\cdot\bm{y})^{\prime}=\bm{p}^{\prime}\cdot\bm{y}+\bm{p}\cdot\bm{y}^{\prime}$
$\displaystyle=-\mathinner{\Bigl{(}\nabla_{\bm{m}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{p}(t)+\nabla_{\bm{m}}r(\bm{m}(t),\bm{\theta})\Bigr{)}}\cdot\bm{y}$
$\displaystyle\quad+\bm{p}\cdot\mathinner{\Bigl{(}\nabla_{\bm{m}}\bm{g}(\bm{m}(t);\bm{\theta})^{\top}\bm{y}(t)+\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})^{\top}\bm{\zeta}\Bigr{)}}$
$\displaystyle=-\nabla_{\bm{m}}r(\bm{m}(t),\bm{\theta})^{\top}\bm{y}(t)+\bm{p}(t)^{\top}\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{\zeta}.$

Since $\bm{y}(0)=\bm{0}$, we know

$\displaystyle\bm{p}(T)\cdot\bm{y}(T)$
$\displaystyle=\int_{0}^{T}\mathinner{\Bigl{(}-\nabla_{\bm{m}}r(\bm{\theta},\bm{m}(t))^{\top}\bm{y}(t)+\bm{p}(t)^{\top}\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})\bm{\zeta}\Bigr{)}}\operatorname{d\\!}t.$
(46)

Substituting (46) into (45) yields

$\frac{\operatorname{d\\!}}{\operatorname{d\\!}\varepsilon}\ell(\bm{\theta}_{\varepsilon})\Big{|}_{\varepsilon=0}=\mathinner{\Bigl{\\{}\int_{0}^{T}\mathinner{\left(\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})^{\top}\bm{p}(t)+\nabla_{\bm{\theta}}r(\bm{\theta},\bm{m}(t))\right)}\operatorname{d\\!}t\Bigr{\\}}}\cdot\bm{\zeta}.$

As $\bm{\zeta}$ is arbitrary, we know that the gradient
$\nabla_{\bm{\theta}}\ell(\bm{\theta})$ is as claimed in (20).

Note that the integrand in (20) is

$\mathinner{\left(\nabla_{\bm{\theta}}r(\bm{\theta},\bm{m}(t))+\nabla_{\bm{\theta}}\bm{g}(\bm{m}(t);\bm{\theta})^{\top}\bm{p}(t)\right)}=\nabla_{\bm{\theta}}H(\bm{m}(t),\bm{p}(t);\bm{\theta}).$

Hence, at the optimal $\bm{\theta}^{*}$ of $\ell$, we have

$\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}\varepsilon}\ell(\bm{\theta}^{*}_{\varepsilon})\Big{|}_{\varepsilon=0}=\mathinner{\Bigl{(}\int_{0}^{T}\nabla_{\bm{\theta}}H(\bm{m}^{*}(t),\bm{p}^{*}(t);\bm{\theta}^{*})\operatorname{d\\!}t\Bigr{)}}\cdot\bm{\zeta}\geq
0,$ (47)

for all $\bm{\zeta}\in\mathbb{R}^{m}$, from which we readily deduce the
identity regarding $H$ at $\bm{\theta}^{*}$. ∎

###  A.5 Proof of Theorem 4

###### Proof.

Let $\bm{v}\in\mathbb{R}^{n}$ be arbitrary and consider the variation
$\bm{u}_{\epsilon}\mathrel{\mathop{\ordinarycolon}}=\bm{u}+\epsilon\bm{v}+o(\epsilon)$
of $\bm{u}$ with $\epsilon>0$. Let $\bm{x}_{\epsilon}(t)$ be the $\bm{x}$-part
of the solution $\bm{m}_{\epsilon}(t)$ to (16b) with initial
$[\bm{u}_{\epsilon};\bm{0}]$. Suppose
$\bm{x}_{\epsilon}(t)=\bm{x}(t)+\epsilon\bm{w}(t)+o(\epsilon)$ for all
$t\in[0,T]$ as $\epsilon\to 0$, then $\bm{w}(t)$ solves

$\begin{cases}\bm{w}^{\prime}(t)=\nabla_{\bm{x}}\bm{g}_{\bm{x}}(\bm{m}(t);\bm{\theta})\bm{w}(t),\quad
0\leq t\leq T,\\\ \bm{w}(0)=\bm{v}.\end{cases}$ (48)

Note that (48) is a linear ODE of $\bm{w}$ and thus has an analytic solution
as follows:

$\bm{w}(T)=e^{\int_{0}^{T}\nabla_{\bm{x}}\bm{g}_{\bm{x}}(\bm{m}(t);\bm{\theta})\operatorname{d\\!}t}\bm{v}.$

Next, we compute the directional derivative of $L$ defined in (26) at $\bm{u}$
along direction $\bm{v}$:

$\displaystyle\frac{\operatorname{d\\!}}{\operatorname{d\\!}\epsilon}L(\bm{u}_{\epsilon})\Big{|}_{\epsilon=0}$
$\displaystyle=\frac{\operatorname{d\\!}}{\operatorname{d\\!}\epsilon}\mathinner{\Bigl{(}\mathcal{R}(\bm{u}_{\epsilon})-\bm{1}\cdot\bm{x}_{\epsilon}(t)\Bigr{)}}\Big{|}_{\epsilon=0}$
$\displaystyle=\nabla_{\bm{u}}\mathcal{R}(\bm{u})\cdot\bm{v}-\bm{1}\cdot\bm{w}(T)$
$\displaystyle=\mathinner{\Bigl{(}\nabla_{\bm{u}}\mathcal{R}(\bm{u})-e^{\int_{0}^{T}\nabla_{\bm{x}}\bm{g}_{\bm{x}}(\bm{x}(t);\bm{\theta})^{\top}\operatorname{d\\!}t}\bm{1}\Bigr{)}}\cdot\bm{v}.$

As $\bm{v}$ is arbitrary, we know the gradient $\nabla_{\bm{u}}L(\bm{u})$ is

$\displaystyle\nabla_{\bm{u}}L(\bm{u})=\nabla_{\bm{u}}\mathcal{R}(\bm{u})-e^{\int_{0}^{T}\nabla_{\bm{x}}\bm{g}_{\bm{x}}(\bm{x}(t);\bm{\theta})^{\top}\operatorname{d\\!}t}\bm{1}.$
(49)

It is clear that the second term on the right hand side of (49) is $\bm{s}(T)$
solved from

$\displaystyle\begin{cases}\bm{s}^{\prime}(t)=\nabla_{\bm{x}}\bm{g}_{\bm{x}}(\bm{x}(t);\bm{\theta})^{\top}\bm{s}(t),\quad
0\leq t\leq T,\\\ \bm{s}(0)=\bm{1}.\end{cases}$ (50)

This completes the proof. ∎

## References

  * [1] Á. Bodó, G. Y. Katona, and P. L. Simon.  SIS epidemic propagation on hypergraphs.  Bulletin of mathematical biology, 78(4):713–735, 2016.
  * [2] M. Boguná and R. Pastor-Satorras.  Epidemic spreading in correlated complex networks.  Physical Review E, 66(4):047104, 2002.
  * [3] C. Borgs, M. Brautbar, J. Chayes, and B. Lucier.  Maximizing social influence in nearly optimal time.  Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, Dec 2013.
  * [4] Q. Cao, H. Shen, J. Gao, B. Wei, and X. Cheng.  Popularity prediction on social platforms with coupled graph neural networks.  In Proceedings of the 13th International Conference on Web Search and Data Mining, WSDM ’20, pages 70–78, New York, NY, USA, 2020. Association for Computing Machinery.
  * [5] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud.  Neural ordinary differential equations.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 6571–6583. Curran Associates, Inc., 2018.
  * [6] A. J. Chorin, O. H. Hald, and R. Kupferman.  Optimal prediction and the mori–zwanzig representation of irreversible processes.  Proceedings of the National Academy of Sciences, 97(7):2968–2973, 2000.
  * [7] S.-N. Chow, X. Ye, H. Zha, and H. Zhou.  Influence prediction for continuous-time information propagation on networks.  Networks and Heterogenous Media, 13(4):567–583, 2018.
  * [8] A. Clauset, C. Moore, and M. E. J. Newman.  Hierarchical structure and the prediction of missing links in networks.  Nature, 453(7191):98–101, May 2008.
  * [9] E. Cohen.  Size-estimation framework with applications to transitive closure and reachability.  Journal of Computer and System Sciences, 55(3):441–453, 1997.
  * [10] E. Cohen, D. Delling, T. Pajor, and R. F. Werneck.  Sketch-based influence maximization and computation: Scaling up with guarantees.  In Proceedings of the 23rd ACM International Conference on Conference on Information and Knowledge Management, CIKM ’14, pages 629–638, New York, NY, USA, 2014. Association for Computing Machinery.
  * [11] X. Dong, D. Thanou, M. Rabbat, and P. Frossard.  Learning graphs from data: A signal representation perspective.  IEEE Signal Processing Magazine, 36(3):44–63, 2019.
  * [12] N. Du, Y. Liang, M.-F. Balcan, and L. Song.  Influence function learning in information diffusion networks.  In Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32, ICML’14, pages II–2016–II–2024. JMLR.org, 2014.
  * [13] N. Du, L. Song, M. Gomez-Rodriguez, and H. Zha.  Scalable influence estimation in continuous-time diffusion networks.  In Advances in Neural Information Processing Systems, pages 3147–3155, 2013.
  * [14] N. Du, L. Song, M. Yuan, and A. J. Smola.  Learning networks of heterogeneous influence.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 25, pages 2780–2788. Curran Associates, Inc., 2012.
  * [15] M. Farajtabar, X. Ye, S. Harati, L. Song, and H. Zha.  Multistage campaigning in social networks.  In Advances in Neural Information Processing Systems, pages 4718–4726, 2016.
  * [16] M. Gomez-Rodriguez, D. Balduzzi, and B. Schölkopf.  Uncovering the temporal dynamics of diffusion networks.  In Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11, pages 561–568, Madison, WI, USA, 2011. Omnipress.
  * [17] M. Gomez-Rodriguez, J. Leskovec, and A. Krause.  Inferring networks of diffusion and influence.  ACM Transactions on Knowledge Discovery from Data (TKDD), 5(4):21, 2012.
  * [18] M. Gomez-Rodriguez, J. Leskovec, and B. Schölkopf.  Structure and dynamics of information pathways in online media.  CoRR, abs/1212.1464, 2012.
  * [19] M. Gomez Rodriguez and B. Schölkopf.  Influence maximization in continuous time diffusion networks.  In 29th International Conference on Machine Learning (ICML 2012), pages 1–8. International Machine Learning Society, 2012.
  * [20] M. Gomez-Rodriguez and B. Schölkopf.  Submodular inference of diffusion networks from multiple trees.  In ICML, 2012.
  * [21] M. Gomez-Rodriguez, L. Song, N. Du, H. Zha, and B. Schölkopf.  Influence estimation and maximization in continuous-time diffusion networks.  ACM Transactions on Information Systems (TOIS), 34(2):1–33, 2016\.
  * [22] A. Goyal, W. Lu, and L. V. Lakshmanan.  Celf++: Optimizing the greedy algorithm for influence maximization in social networks.  In Proceedings of the 20th International Conference Companion on World Wide Web, WWW ’11, pages 47–48, New York, NY, USA, 2011. Association for Computing Machinery.
  * [23] J. Gulddahl Rasmussen.  Lecture Notes: Temporal Point Processes and the Conditional Intensity Function.  arXiv e-prints, page arXiv:1806.00221, June 2018.
  * [24] Q. He, X. Wang, Z. Lei, M. Huang, Y. Cai, and L. Ma.  Tifim: A two-stage iterative framework for influence maximization in social networks.  Applied Mathematics and Computation, 354:338–352, 2019.
  * [25] S. He, H. Zha, and X. Ye.  Network diffusions via neural mean-field dynamics.  In Advances in Neural Information Processing Systems 33, 2020.
  * [26] D. Kempe, J. Kleinberg, and É. Tardos.  Maximizing the spread of influence through a social network.  In Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 137–146. ACM, 2003.
  * [27] D. Kempe, J. Kleinberg, and É. Tardos.  Influential nodes in a diffusion model for social networks.  In Automata, languages and programming, pages 1127–1138. Springer, 2005.
  * [28] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  International Conference on Learning Representations, Dec 2014.
  * [29] T. N. Kipf and M. Welling.  Semi-supervised classification with graph convolutional networks.  In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net, 2017.
  * [30] J. Leskovec, L. Backstrom, and J. Kleinberg.  Meme-tracking and the dynamics of the news cycle.  In Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’09, pages 497–506, New York, NY, USA, 2009. Association for Computing Machinery.
  * [31] J. Leskovec, D. Chakrabarti, J. Kleinberg, C. Faloutsos, and Z. Ghahramani.  Kronecker graphs: An approach to modeling networks.  The Journal of Machine Learning Research, 11:985–1042, 2010.
  * [32] J. Leskovec, A. Krause, C. Guestrin, C. Faloutsos, J. VanBriesen, and N. Glance.  Cost-effective outbreak detection in networks.  In Proceedings of the 13th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’07, pages 420–429, New York, NY, USA, 2007. Association for Computing Machinery.
  * [33] J. Leskovec, K. J. Lang, A. Dasgupta, and M. W. Mahoney.  Statistical properties of community structure in large social and information networks.  In Proceedings of the 17th International Conference on World Wide Web, WWW ’08, pages 695–704, New York, NY, USA, 2008. Association for Computing Machinery.
  * [34] J. Leskovec and R. Sosič.  Snap: A general-purpose network analysis and graph-mining library.  ACM Transactions on Intelligent Systems and Technology (TIST), 8(1):1, 2016.
  * [35] C. K. Leung, A. Cuzzocrea, J. J. Mai, D. Deng, and F. Jiang.  Personalized deepinf: Enhanced social influence prediction with deep learning and transfer learning.  In 2019 IEEE International Conference on Big Data (Big Data), pages 2871–2880, 2019.
  * [36] C. Li, J. Ma, X. Guo, and Q. Mei.  Deepcas: An end-to-end predictor of information cascades.  In Proceedings of the 26th international conference on World Wide Web, pages 577–586, 2017.
  * [37] Y. Li, J. Fan, Y. Wang, and K.-L. Tan.  Influence maximization on social graphs: A survey.  IEEE Transactions on Knowledge and Data Engineering, 30(10):1852–1872, 2018.
  * [38] Y. Liang, Z. Jiang, and Y. Zheng.  Inferring traffic cascading patterns.  In Proceedings of the 25th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems, SIGSPATIAL ’17, New York, NY, USA, 2017. Association for Computing Machinery.
  * [39] B. Lucier, J. Oren, and Y. Singer.  Influence at scale: Distributed computation of complex contagion in networks.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’15, pages 735–744, New York, NY, USA, 2015. Association for Computing Machinery.
  * [40] S. Manchanda, A. MITTAL, A. Dhawan, S. Medya, S. Ranu, and A. Singh.  Gcomb: Learning budget-constrained combinatorial algorithms over billion-sized graphs.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 20000–20011. Curran Associates, Inc., 2020.
  * [41] G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro.  Connecting the dots: Identifying network structure via graph signal processing.  IEEE Signal Processing Magazine, 36(3):16–43, May 2019.
  * [42] J. C. Miller and I. Z. Kiss.  Epidemic spread in networks: Existing methods and current challenges.  Mathematical modelling of natural phenomena, 9(2):4, 2014.
  * [43] S. A. Myers and J. Leskovec.  On the convexity of latent social network inference.  In Proceedings of the 23rd International Conference on Neural Information Processing Systems - Volume 2, NIPS’10, pages 1741–1749, Red Hook, NY, USA, 2010. Curran Associates Inc.
  * [44] M. Newman.  Networks: an introduction.  Oxford University Press, 2010.
  * [45] H. T. Nguyen, T. P. Nguyen, T. N. Vu, and T. N. Dinh.  Outward influence and cascade size estimation in billion-scale networks.  Proc. ACM Meas. Anal. Comput. Syst., 1(1), June 2017.
  * [46] G. Panagopoulos, F. Malliaros, and M. Vazirgiannis.  Multi-task learning for influence estimation and maximization.  IEEE Transactions on Knowledge and Data Engineering, pages 1–1, 2020.
  * [47] G. Panagopoulos, F. D. Malliaros, and M. Vazirgianis.  Influence maximization using influence and susceptibility embeddings.  Proceedings of the International AAAI Conference on Web and Social Media, 14(1):511–521, May 2020.
  * [48] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani.  Epidemic processes in complex networks.  Reviews of modern physics, 87(3):925, 2015.
  * [49] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [50] C. Qiao, Y. Shi, Y.-X. Diao, V. D. Calhoun, and Y.-P. Wang.  Log-sum enhanced sparse deep neural network.  Neurocomputing, 407:206–220, 2020.
  * [51] J. Qiu, J. Tang, H. Ma, Y. Dong, K. Wang, and J. Tang.  Deepinf: Social influence prediction with deep learning.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 2110–2119, 2018.
  * [52] Y. Rong, Q. Zhu, and H. Cheng.  A model-free approach to infer the diffusion network from event cascade.  In Proceedings of the 25th ACM International on Conference on Information and Knowledge Management, CIKM ’16, pages 1653–1662, New York, NY, USA, 2016. Association for Computing Machinery.
  * [53] F. D. Sahneh and C. Scoglio.  Epidemic spread in human networks.  In Decision and Control and European Control Conference (CDC-ECC), 2011 50th IEEE Conference on, pages 3008–3013. IEEE, 2011.
  * [54] Y. Tang, Y. Shi, and X. Xiao.  Influence maximization in near-linear time: A martingale approach.  In Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data, pages 1539–1554. ACM, 2015.
  * [55] Y. Tang, X. Xiao, and Y. Shi.  Influence maximization: Near-optimal time complexity meets practical efficiency.  In Proceedings of the 2014 ACM SIGMOD International Conference on Management of Data, SIGMOD ’14, pages 75–86, New York, NY, USA, 2014. Association for Computing Machinery.
  * [56] M. Vergeer, L. Hermans, and S. Sams.  Online social networks and micro-blogging in political campaigning the exploration of a new campaign tool and a new campaign style.  Party Politics, 19(3):477–501, 2013.
  * [57] L. Wang, S. Ermon, and J. E. Hopcroft.  Feature-enhanced probabilistic models for diffusion network inference.  In Machine Learning and Knowledge Discovery in Databases, pages 499–514. Springer, 2012.
  * [58] J. Wortman.  Viral marketing and the diffusion of trends on social networks.  Technical Reports (CIS), page 880, 2008.
  * [59] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and P. S. Yu.  A comprehensive survey on graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–21, 2020.
  * [60] Y. Zang, G. Bao, X. Ye, H. Zha, and H. Zhou.  A jump stochastic differential equation approach for influence prediction on heterogenous networks.  Communications in Mathematical Sciences, 18(8):2341–2359, 2020\.
  * [61] Y. Zhu, J. Xie, and Z. Chen.  Predicting the popularity of micro-videos with multimodal variational encoder-decoder framework.  arXiv preprint arXiv:2003.12724, 2020.