Diffusion Network Inference for Cross-layer Cascades
Authors: Siyu Huang
Yubai Yuan
Abdul Basit Adeel
Venue: NeurIPS 2025 poster ()
Keywords: diffusion networks
low-rank optimization
network mixture model
social networks
variational inference
Primary area: machine_learning_for_sciences

Abstract: A cascade over a network refers to the diffusion process where behavior changes occurring in one part of an interconnected population lead to a series of sequential changes throughout the entire population. In recent years, there has been a surge in interest and efforts to understand and model cascade mechanisms since they motivate many significant research topics across different disciplines. The propagation structure of cascades is governed by underlying diffusion networks that are often hidden. Inferring diffusion networks thus enables interventions in cascading process to maximize information propagation and provides insights into the Granger causality of interaction mechanisms among individuals. In this project, we propose a novel double network mixture model for inferring latent diffusion network in presence of strong cascade heterogeneity. The new model represents cascade pathways as a distributional mixture over diffusion networks that capture different cascading patterns at the population level. We develop a data-driven optimization method to infer diffusion networks using only visible temporal cascade records, avoiding the need to model complex and heterogeneous individual states. Both statistical and computational guarantees are established for the proposed method. We apply the proposed model to analyze research topic cascades in social sciences across U.S. universities and uncover the latent research topic diffusion network among top U.S. social science programs.

Diffusion Network Inference for Cross-layer Cascades
Siyu Huang* Department of Statistics
Yubai Yuan" Department Statistics
The Pennsylvania State University sph6006@psu.ed edu
The Pennsylvania State University
yvy5509Qpsu.ed edu
Abdul Basit Adeel
Department of Sociology and Criminology The Pennsylvania State University axa63720psu. edu
Abstract
A cascade over network refers to the diffusion process where behavior changes occurring in one part of an interconnected population lead to series of sequen- tial changes throughout the entire population. In recent years, there has been surge in interest and efforts to understand and model cascade mechanisms since they motivate many significant research topics across different disciplines. The propagation structure of cascades is governed by underlying diffusion networks that often hidden. Inferring diffusion networks thus enables interventions in cascading process to maximize information propagation and provides insights into the Granger causality of interaction mechanisms among individuals. In this project, we propose a novel double network mixture model for inferring latent diffusion network in presence of strong cascade heterogeneity. The new model represents cascade pathways as a distributional mixture over diffusion networks that capture different cascading patterns the population level. We develop a data-driven optimization method to infer diffusion networks using only visible temporal cascade records, avoiding the need model complex and heterogeneous individual states. Both statistical and computational guarantees are established for the proposed method. We apply the proposed model to analyze research topic cascades in social sciences across U.S. universities and uncover the latent research topic diffusion network among top U.S. social science programs.
1
Introduction
Cascades over network refer to the diffusion processes where behavior changes in part of an interconnected population lead to series of sequential changes throughout the entire population. In recent years, there are surging interests and efforts to understand and model the cascade mechanism since it motivates many significant research topics in different areas, including social influence 12 13 20], information propagation via social media [32 [] diffusion policy and social norms viral marketing 25 9], and contagion of infectious diseases 21 35].
4
33
understand the diffusion networks that govern cascade propagation patterns. However, diffusion networks are often hidden and need to be inferred from observed cascading behaviors. For example, in the case of infectious diseases, we can observe when an individual is infected but need to impute the missing information on who infects this individual.
One fundamental problem
*Co-first authors. Corresponding author.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

More importantly real-world cascading behaviors often exhibit strong heterogeneity and are jointly governed by different diffusion patterns. For example, in information cascade and epidemiology. cascades can diffuse among population via different transmission channels building on various social relations, and thus lead to heterogeneous propagation speeds and scales [24 30 Furthermore, cascading heterogeneity originates from the variability of individuals' statuses in engaging the cascade
[39 [T. In social media like X and Instagram, the spreading pattern and the speed of messages
heavily depend on users' activity status. An active user will respond more instantly to interesting messages and accelerate the information cascade compared to inactive users [TO]. As for volatility cascades in financial markets, structure heterogeneity in volatillity diffusion depends different time horizons of the agents the market [41]. In these scenarios, individual statuses determine the transmission channels engaging in the cascades, and changes of individual statuses can also change the downstream cascade diffusion patterns. To conclude from these examples, diffusion patterns can exhibit combinatorial complexity population grows. We present example of cascade diffusion when transmission channels of individuals vary in Figure (
(a) Observed network
(b) Latent diffusion network
Figure cascade diffuses through nodes {A, B, D,C.D,E, F} on two-layer network. Yellow nodes: activated nodes; gray nodes: inactivated nodes. {B, D} are activated via the observed network and (E, C, F} are activated via the latent network.
In this paper, we propose a novel double network mixture model to infer multiple diffusion networks simultaneously from heterogeneous cascade data. The proposed model introduces distributional mixture of diffusion networks to capture the heterogeneous cascading patterns, where diffusion networks provide complementary connection information. The main advantage the distributional mixture is to avoid modeling the complex individual status changes. Specifically, the proposed model can describe the diffusion process over multi-layer networks where cascades propagate across different layers alternatively. Compared to existing methods, the proposed method can uncover latent diffusion networks even when the number of diffusion patterns is exponential to the number of nodes. Furthermore, the parameter estimation in model can be solved by sequence of convex optimization problems, which leads to both statistical and computational guarantees for our diffusion network estimation.
2 Related works
Various directed probabilistic graphical models have been developed to infer diffusion networks from observed cascade samples 11
g 22 18]. Generally, these models treat infection time continuous random variable and construct the likelihood of cascade samples based on the target diffusion network under local Markov assumption. To capture heterogeneous diffusion patterns, several multi-pattern cascade models have been proposed [37 40) where cascade samples are adaptively clustered into groups and each group corresponds to a distinct diffusion network. Among these methods, ConNIe (2), NetRate [T], and MMRate (77 are popular representatives. Specifically, ConNIe employs a maximum likelihood formulation via convex programming, incor- porating an 11-type penalty to promote sparsity the inferred network. Building upon ConNle, NetRate explicitly represents diffusion as continuous-time probabilistic process. characterized by edge-specific transmission rates governing edge-wise diffusion probabilities. MMRate further extends this framework by accommodating multiple distinct diffusion patterns, assuming multiple heterogeneous latent networks, with each cascade diffusing via one network according to a certain probability.
19
31

---

3 Methodology
3.1 Continuous-time cascade on single network
Consider a network with N nodes where each node has two infection conditions in a single cascade: infected (activated) and uninfected (inactivated), and an infected node will always remain infected. A cascade a N-dimensional temporal record t = (4,,22 ... AN), where I is the infection time of the i-th node. Instead of infinite time horizon, length T, i.e., Vi, 4 E [40, ln + TJ. where lo = minj si v {4, is the infection time of source node. We denote the infection time of the nodes not infected in the observation window to + T. Without loss of generality, we assume to 0 in this paper.
observe cascade within finite time window of
A cascade difffuses nodewisely over edges of a diffusion network. The continuous-time model formulates the cascade transmission from an infected node j to another infected node with survival analysis models. Specifically, given node j being infected at time ty. let f(a ty, AyE) denote the likelihood of node i being infected by node j at time t., where t; > t,. The transmission rate λÿá represents how fast cascades diffuse from node node i. Accordingly, the cumulative probability S
function is F(ts Eg.
f(s t,, Aj:)ds We consider the hazard rate function defined as
5(4 (44, S(t, t, Aj) exp (-f," H(t (t,, Ag)) 1 F(t, t,, Ags) S(t, Ay,)
1(4 Ay,
(-( H(()t))dt)
where S(4 L,, 23s) = - 1-F F(4 lj Azs) is the survival function that indicates the probability of node not being infected by node j until time 4-- Typical parametric forms of hazard rate functions are Exp model: H(t łj, Aye) = Ay; Pow model: H (t, A,,) = - Ray model: H4: (jx Ajs) Xys(4; (j). Given that Ayz controls the likelihood and speed of transmission between node to node i, the global cascading pattern over network can be captured by matrix = (Aig) e R N where rows represent sender nodes and columns represent receiver nodes. Note that A can be asymmetric, i.e., Aij * Xjjt when node i and j different in the capacity of launching transmission.
The diffusion process is typically modeled by independent cascade models [22 where any infected nodes can infect a node independently and a node stays infected once another node infects Therefore, one formulate the likelihood of node being infected by potential parent nodes (j bi:tz < t.) at time t, as
Pt)t;
Σ f(t. bg Azs) juytt
II Ai
S(ta te Aus)
(I)
Σ f4
H(t:
ty Age)
II
S(t, ARi)
On the other hand, if node is not activated by any parent nodes in the observation window, then the corresponding likelihood for {i: 14 7) is
(T; P(t. (ty < <t))
II S(T Ag)) Jit, <T
(2)
Since infections of each node are conditionally independent given the corresponding parent nodes, we can decompose the likelihood of a cascade t as the multiplication of series of conditional probabilities:
N A(L:) I P41 A.) 1 < (T)((A))))))))))) 1(4 T)P:(T;A..) 1=1
Therefore, the full likelihood of independent cascade samples (1())) is [T_, P( θ). The advantage of the above model is that each cascade sample induces a directed acyclic graph, whose local Markov property allows the decomposition of likelihood and thus reduction in computational complexity of inferring the diffusion network O
3

---

3.2 Double network mixture model
In many applications such as social media, there exist multiple networks among an interconnected population that reflect different types of relations. Cascade diffusion pattern can change over different networks, and cascade can proceed on different networks alternatively and simultaneously due to the inter-layer interactions. The diffusion process on multi-layer networks can greatly increase the degree freedom possible diffusion pathways, which leads to strong heterogeneity in cascade observations. To model the multi-network cascade behavior, we propose the double network mixture model. Consider two diffusion networks O and v among the same population with N units, where Oj and Vi ij are the transmission rate from node i to node corresponding to two relations. We introduce double diffusion indicators 5 E {0, 1] such that
z iid Bern (T:), 1, c.,,c......
For cascade c, node activated via diffusion pathway on network O 2 1, and 2 0f via network V. Denote π { [0, 11 where is the probability of node i engaging cascade via network O Then the diffusion pathway of the cascade e can be represented as
E RR NxN = ( + (1- "
O + (1 ZN v NN c= 1,..,,,, C. (3)
Therefore, -r column-wiss mixture of @ and V, which can vary for different cascades. Different to conventional mixture model, the proposed method allows Zí to vary across different cascade samples and nodes. Therefore, the proposed model (3) can generate up to 2 different types of diffusion patterns, which grows exponentially as network size increases. Therefore, the proposed method has the principled heterogeneity modeling for the diffusion patterns. The likelihood of the proposed double mixture model can be also explicitly formulated. Following Pı and Po in [ we have the probability for node i being infected in the c-th cascade given specific network
P( Z5: U1 15 < (1) θr)
(4); {.)] 27
[Pr(4 V V.) V.)
(4)
and following 2 the probability for node not being infected is
P(T Z :t < -
(T;
(5)
N.C Denote Z = (2 c= the joint distribution of diffusion indicators is
P(Z) II P(Z))
2 (1 m
(6)
c=1
1=1
Then we have the joint distribution of cascade samples and diffusion indicators Z as
11 P Z:
II C-1
I iT
P. ( 4 [11-[1]]]] PT U.) 1-)
(7)
cml
I juza
(T 25 ( (o 12
(8)
where Ω (TT, ① Wt denote model parameters. And the marginal distribution of cascade samples is
1| P P( 10
11
Il SEE
T. -1-
(9)
T.Po + +(1++++1110111111))) - (I
Another advantage of our model is that the posterior distribution of diffusion indicator {2 can be calculated in an explicit form
ift <
# : = -1
T.
(10)
if
P.

---

3.3 Model identification
In this subsection, we establish the identifiability of the double network mixture model. Denote R == (j (1...... | O ji > 0} and R² == {j E {1. [1...... .N L js > 0} as the sets of nodes that can directly reach node i on O and V, and tR. [o, T]jR UR := [t, | R. u R as the infectious times of these nodes. We have the following identifiability result for model in (9.
Proposition 1. For each node i, i= 1,..., ,N, assume that 1) |[Ox|l1 + | > 0 and there exists j + i such that ju * Lú, and 2) survival function satisfies S(t, tg. Asa) = exphty;h(t, - t,)] for some differentiable function h(.)). Then, the parameters (T;; O v.:)) associated with node in 9 are identifiable, i.e., if
(iP O.) + (1 (L- 7; 3 Pt; L.): or T. (t; ,) + -1-18 (1- Pu v.) = n. Po (t: + +1 n. (o (.:)
for any t2 and bR.. then n. = n. O. . and L V.
Assumption 2) can be satisfied by popular risk models including Exp model, Pow model, Ray model, and other additive risk models of information propagation, such kernel hazard functions [10 and
feature-enhanced hazard functions [36]. Proposition I shows that the network preferences { are identifiable and diffusion networks + and are column-wiss identifiable. However, similar to the labeling non-identifiabiiit issue in finite mixture model [23], O and y may still not be globally identifiable without structural constraints due to column permutation. Specifically, the data distribution does not change if we swap O and L. in O and & with other columns fixed.
Lpeer-specific network structure constraint Motivated by world applications, one can interpret @ as the diffusion pathways over an observed social network A E (0, N)NN Therefore, we can impose support constraints θ
O > 0f A.; 1; ,, = 0 A.j 0.
Due to the sparse nature of real-world social networks, the support constraint also implicitly imposes sparsity constraint on + On the other hand, one can interpret V as the diffusion pathways via latent social relations of individuals that are not captured by the social network A. The magnitude of L ij reflects the social distance between individual and in terms of their latent factors. It has been found that social distance typically has or can be approximated by low-rank structure << 266 22]. since gigh-dimensional social factors can always be embeded into a low dimensional latent space that preserves social distances (34]. Therefore, we impose low-rank structure on V as
rank V r, 1 << N.
Imposing the above support constraint and low-rank structure allows 2 and V capture comple- mentary diffusion patterns driven by different types of relations. In addition, the structure constraints solve the above column-wise permutation issue [6]. Specifically, we introduce the matrix sub- w space A1(O) = {Ne N support N) = A). We also perform SVD on V = UVVT
where U, V E R and r is the rank of L. Then, we define another matrix subspace as A2( ( = (UX++YV
X,Y e R-Nkk) J. We have the following result:
Proposition 2. Given that the assumptions in Proposition I hold, then O and & are identifiable if: (11)
max (M)
|N1]2 x
max NEA2 NeAi(
||N < 1,
where || ||2 and -11zo denote matrix operation norm and largest element in magnitude.
The first term in (33 controls the rank of θ given fixed sparsity level, where a larger value indicates a lower rank. The second term controls the sparsity level of V given a fixed rank, where larger value indicates a lower sparsity level. Intuitively, networks @ and V can be globally identified given that they are well-separated in terms of either rank or sparsity.
3.4 Model estimation
Combining the distribution of cascade samples in Section 3.2 and the network structure constraints in Section 3.3, we estimate the model parameters f via constrained likelihood maximization as
1 C
> log P(t()): (2) s.t. + (I A)= 0, rank ( < r,
arg max Ω- ( LO.VVR )
C
cml
5

---

where I is N by--N matrix with all elements being 1. However, both the likelihood function and the rank regularization are difficult to directly optimize. Therefore, we maximize the evidence lower bound of the log-likelihood function and replace the low-rank penalty with its convex relaxation as nuclear norm || Il.. The optimization problem can thus be reformulated as follow
arg max 100mmax222 (z) l C (()).2
Σ log P(t 3 [Ba(z) log q( [Z) s.t. O A) 0,
& P.
where ZZ( is the expectation of Z over distribution ( The above objective function can be optimized via EM algorithm by iteratively updating ( and Q. Specifically, with the estimated parameters (2(*) from the s-th step:
N 11 P(Z
Gestep:q(Z:9<4
c=1i=1
M-step: arg max
Σ C
1
E P log P(t Z: s.t. @ A) 0, | <.p
Q
In E-step, the posterior distribution of network indicators #े = P(Z ; (())) can
be explicitly updated via (10). In M-step, the objective function Q
2
2o_1 E F
[lo log P(t(o), (() can be decomposed Q ( Q Ω )( +
Q2( (2()))++ + Qa( | (2())) The arguments θ, v and π can thus be updated parallelly in M-step
M.1 : (++1) arg max Qi(e)
2 (+)
A) s.t. 1 A) -0
(12)
M.2 (0++1)_ arg max Q2(4V ( Ω (+)) s.t. 1 &p
(13)
V
M.3 π (++1) arg max Q3 (T (2()))
(14)
When the structural network A is not observed, one can impose nnrr penalty to pursue sparsity structure in θ. Accordingly, the constraint in optimization (12) is replaced by |0]| & s', where > 0 is the sparsity tuning parameter. The main advantage of the proposed relaxation that the M-step becomes series of convex optimization problems. Specifically, denote parameter spaces (s := O E [o, 31 N XN
O A) 0} with s = A o. Ex()) := ( [0, [2] XN | p). and + := [ E [e, c JM where B1, B2 are nonnegative constants and O < e < 0.5.
We have the following result:
Theorem 3.1. The parameter spaces (s) (p),,, and Ex are convex sets for any 0, p
and Q3(T (2()) is concave on π. If the hazard function H (²,A) satisfies (H (2) = 0 for 0x
t t t', then @i(@ ()(e)) and Q2 L (2(=)) concave θ and L respectively. Furthermore, if for any node i, the probabilities of being source node P(v) > 0 for e R. where R denotes the set of
nodes from which is reachable via directed path, then E, Q ( [(())]] and E, Ee(Q2(V) ( are also strictly concave in terms of O and N respectively.
Theorem 3.1 guarantees that the optimization in M-step has a unique and optimal solution. Combining the theorem with the convergence guarantee of EM algorithm for convex ancillary function Q(
( 38, the above likelihood maximizer Ω is guaranteed to converge to the true Ω. addition, the convexity assumption on hazard function can be satisfied by popular risk models including Exp model, Pow model, Ray model, and other additive risk models. We summarize and provide details for the above optimization process in the Appendix.
4 Numerical experiments on synthetic cascading data
We investigate the performance of the proposed method recovering the diffusion networks based on synthetic cascading data. The performance of global transmission rates estimation on a network O is measured by normalized mean absolute error (MAE) MAE(@) Σx 1 +0y1/O, The

---

performance of network probability estimation is measured as MAE(T)) = - In addition, we investigate the performance recovering the structure of diffusion network. Specifically, given a diffusion network ④ we consider three topology estimation metrics accuracy, precision, and recall as Ace( ()) = E.. [ .y) I(0,,)l/(2 I( + b.. I(),))) Pre(O) = ش,, (I)),)) IO, Σi,; I(@iz), and Recall () = Σ,, (T),)) I(O,))) E., I(O,)) where I(a) - 1 if
> 0 and I(a) = 0 otherwise. In the following numerical experiments, we fix the size of diffusion networks at N 200.
4.1 Benchmark comparison under different network topologies
We compare the proposed method with baseline methods including NetRate 31 MMRate 37 and
ConNle [29] on diffusion network recovery. We generate diffusion networks and vi to mimic different structure of real-world networks. Specifically, we consider the diffusion network @ and its support as random network, network with community structure, and scale-free network. On the other band, we fix the latent diffusion network to be both low-rank (rank 5) and sparse (edge density 0.05). Given and V,, we generate C = 2, 000 independent cascade samples based on double mixture model with Exp transmission model and observation window length being T 10.
Table shows the proposed method outperforms both NetRate and MMRate by achieving lower MAE of transmission rate estimation on diffusion network @ under three network topologies and lower MAE of corresponding latent diffusion network V MAE comparison does not include ConNle since it only estimates network topology. In addition, we compare the proposed method with baseline methods in terms of network topology recovery on both & and joint network Ô U Vi for a fair comparison since benchmark methods NetRate and ConNIe do not distinguish the different diffusion networks by design. Table D shows the proposed method achieves higher accuracy on topology recovery via both latent network v and joint network θ V than NetRate, MMRate, and ConNIc under different underlying structures in O It also shows that under different settings, our method significantly outperforms all other methods in precision while remains higher recall than MMRate and ConNle, and achieves similar recall to NetRate. Since only MMRate differentiates diffusion networks and estimates probabilities of cascading diffusing over specific network, compare MAB( T from the proposed method with that from MMRate. Table [ shows that our method also outperforms MMRate in estimation of network selection probability under different network topologies.
Table 1: Diffusion network estimations from different methods under three θ topologies (standard deviation in parenthesis, best performance highlighted in blue).
MrEE
MdE.
MAE-
Arce
Pres
Recy
Arceu
Pree_*
Recevs
Rand
0.315/0.311)
0.5440.0007
(1047(0,000)
0.842/0.
0.757
0 (06)
0.889
0.82900.006)
0.955
OurAlg
Com
0.264:0.005)
0.507(0.001
0.02
0.836:0.005)
0.767(0.007)
0.912/0.004)
0.888/0.001)
0.832/0.006
0.954C0-003)
Scale
0.312/0.006)
0.545(0.006)
0.050(0.002)
0.833/0.005
0.75000.007)
0.937V0.005
0.883(0.004)
0.826:0.007)
0.949X0.003)
Rand
0.865/0.001)
2.9980.006)
0.172/0.002)
0.09(00.001)
0.90460.005
0.116/0.022) 0.207(0.003)
0.842/0.004)
NaRate
Com
0.851c0.000;
3.0030.014)
0.18000 002)
0.09%0.0013
0975500004)
0.216(0.002)
0.12200.011
0.955(00015)
Scale
0.855/0.000)
2.8970.018)
0.17940.002)
0.09%0.0011
0.929(0.003)
0.220/0.002)
0.125/0.001)
0.922/0.013)
Rand
0:637/0.007)
1.900000000
0,888(0 133)
0,201 023)
0.13000.015)
063350.023)
0.265/0 028)
0,164(0.019)
0.6050.019
MMRase
Cam
0.62000.000)
(665600036)
177)
0.219(0.017)
0.13000.012)
0.705/0.008)
0.28100.020)
0.17500015)
0.75500000)
Scale
0.665(0.006)
1.47((00227)
0.928(0.028)
0,211(0.019)
0.126(0.013)
0.65000.009)
0.278(0.023)
0.173:0.018)
0.710000007
Rand
0.519(0:006)
0.397(0.805)
0.722000100
0.550/0.006) 0.638/0006))
0.80088008)
ConNle
Com
0.519(0.005)
0.395/0.005)
0.748/0.00 007)
0.63/000055
0.532/0.007>
0.79600000)
Scale
0.523/0 006)
0.40200.005)
0.748/0.009)
0.646(0.006)
0.543/0.006)
0.79700006)
4.2 Network recovery under different transmission models
In this subsection, we investigate the performance of our diffusion network estimation when cascade samples are generated from popular transmission models including Exp, Pow, and Ray models, respectively. Table 2 illustrates the transmission rates recovery and network topology recovery on w under different transmission models and cascade sample sizes C. As the sample size C' increases. both the parameter estimations (MAE) and topology recovery metrics (Acc, Pre, Rec) improve under different transmission models. In addition, the degree of improvement decreases as more cascade samples become available. This pattern indicates the consistency of the proposed diffusion network estimators, and the convergence of the proposed EM-type optimization. Notice that the diffusion network recovery based on cascade samples generated from Pow transmission model is better than
7

---

Exp and Ray model. This is because the Pow model introduces an additional parameter as time lag lower bound, which lowers the variation of activate time and overestimation transmission rates.
Table 2: Diffusion network estimations from the proposed method under different cascade models and sample sizes.
Mr&Ep
MAE.E.
Arce
Pree
Rice
C- 500
0.8710.019) 0.320)0.19)
0.564/0.010) 0.695(0.06)
0.431(0.011)
0.817(00110)
C= 1000
0.23510.013)
0.66700.077)
0.091:0.005
0.745/0.009)
0.65300.01))
0.862/0.007)
1500
0.216/0.012)
0.61:00000))
0.08500.014)
0.11200006)
0.905/000000) 0.746000008)
0.18700.11)
0.562/0.007)
0.062/0.004)
0.85700.044)
0.794(0.006)
0.9300.005)
0.28710.018)
L05310056)
0.15210.008)
0.605/0.118)
0.52548:024
0,714(0013)
Ray
C= 1000
1500
0.179/0.011)
0.9000 0.211000155
0.62910.91))
0.79500.0175 0.134)0.04)
0.11200019) 0030)
0.702/0.026)
605000026)
0.915/0.007)
06660/811 7)
C- 2000
0. (1188)00112)
0.60500081)
0.13710.004)
0.807(0.113
D. 0.709000119
0.937((0011)
C-500 500
0433100012) 0.171(0.008)
0114200110)
0.862000115)
0.832/0.024)
C= 1000
0. 01129000088
0.32800007)
0.956/0.004) 0.13100.00)
0.932/0.007)
0.962/0.004)
Pow
C= 1500
0111
29000 004)
0.969(0.007) 0.12000008)
0.94600.044)
0.99%(0.007)
2=0000
0,106/0.006)
0.21100004)
0.1290.011)
0.973/0.003
0,9500.005)
0,997(0.001
5 Real cascading data analysis
In this section, we study the cascading patterns of research topics in sociology by discovering the diffusion networks among US universities. Geographic proximity is known to facilitate idea exchanges through collaborations and citations among colleagues within the same institution or nearby locations (3, which can be due to dependence shared research resources and the need for coordination (28). It is also known that research diffusion can happen via the latent network of scholars, which is also sometimes called the "invisible colleges" 7 to highlight the role of informal networks and latent ties. Our goal is to infer the latent research topic diffusion network and compare its difference with the geographic network in terms of topic cascading.
Data preparation and preprocess We select universities in both the list of 1965 ASA Guide to Graduate Departments of Sociology and the list of 2022 US News Best Sociology Programs in America. In addition, we exclude the universities whose sociology programs were established after 1965. Based the above conditions, we finalize N = 104 universities, which considered the target population of our study. We then create a geographical network E (0, 1) N among the selected universities, where Asj- = 1 if university and j are located in the same state.
(a) Geographical network
(b) Latent diffusion network
Figure 2: Cascade transmission rates over geographical network C and latent diffusion network V.
To construct research topics, we use Elsevier's Scopus API to compile a dataset of 29,725 unique articles from 23 top generalist sociology journals, including information authors, their affilia- tions, article titles, and abstracts. Multiple keywords are extracted from the titles and abstracts of articles. After these preprocesses, fields in sociology. For each research topic c, we construct the corresponding cascade samples as (19, 15...., ,t%) where if := E - fgg with ig heing the publication date of the first article that involves topic c, and E being the publication date of the first article that involves topic and is published by any scholar affiliated with university i. If university i never publishes any article that involves topic c, we set is to be year 2022.
obtain 3,033 unique research topics that cover major research
Diffusion networks inference We simultaneously estimate the geographic diffusion network and the latent diffusion network y based on the proposed model. Figure (2) illustrates the inferred diffu-
8

---

sion networks ② and L respectively. The size of nodes represents the out-degree of corresponding universities on networks. We see that the geographic diffusion network has a strong local community structure, especially among universities in California, the Great Lakes, and the northeast coast. Figure (2) also shows that the latent diffusion ner- work demonstrates a decentralized structure and contains many connections between the east and
Disributiinn of Transmissoon Rates Geographic vs. Networks
west coasts. This captures the national collabo- ration and academic mobility. We also compare estimated pairwise transmission rates on diffusion networks θ and & in Figure (3). In summary, the latent diffusion network is denser than the geo- graphic network, and the magnitude of transmis- sion rates on the geographic network is larger and has more variation compared to the latent diffusion network.
Courts
Pairwise Transmission Rate
To provide an interpretation of the inferred geo- graphic diffusion network @ and latent diffusion network L, we connect universities' positions on the network with universities' U.S. News Univer- sity rankings the sociology program. We choose these rankings since they are systematic and popular summary of academic factors. Figure (4a) illustrates the association between program ranking and the node-wing wise betweenness centrality latent diffusion network. The Betweenness centrality of a node i is defined as Σ مجزز (i) where jı is the total number of shortest paths from node to k. and (i)( is the number of those paths that pass through i. Based on the figure. the universities'
Figure 3: Comparison of distributions of pairwise transmission rates from geographic network O, and latent network ij-
University Rank vs, Betweenness Centrality
University Rank vs. PageiRan Centrality
0.25
Betwoornesss Contrality
0.20
Pagokank Contrality
University Ranik
0.15
0110-
0.05
9.00
000000
University Rank
University Rank
100
(a) Betweenness centrality vs Ranking on T
(b) Pagerank centrality Ranking on T
Figure 4: Node-wise centrality vs Ranking on latent network V.
betweenness centralities are relatively uniform and do not significantly decrease as ranking increases. This suggests that there does not exist any strong community or centralized topology in the latent network. When universities are grouped according to their ranking (top 20, middle 60, low 20 etc.), each group has different high-betweenness universities, serving as the bridges in idea diffusion. We also investigate the relation between university rankings and Pagerank centrality. Pagerank centrality measures the influence of a node on a network based on how many influential nodes it connects to. Figure 4bb) shows that the universities' Pagerank centralities decrease their rankings increase, which suggests that the idea exchanges among high-ranked universities are more frequent and faster than those of other universities. This pattern appears because high-ranked universities have a higher level of research activities. University prestige is also considered as an important proxy for the quality of ideas, which increases the likelihood of research ideas from higher-ranked universities to diffuse.
Scalability analysis
The implicit convex nature of our objective function allows the optimization process to be as efficient as gradient descent, even though the optimization of our method is based on EM-type update. The E-step in our algorithm has explicit form in (I0). Furthermore, have proved in Theorem

---

3.1 that the objective function Q( in the M-step is strictly concave and in each EM-iteration, we only perform gradient ascent once instead of maximizing Q(22) until convergence. Under these conditions, established proofs in 2 show that the above EM-update is equivalent to gradient decent, and our optimization thus enjoys gradient descent's s geometrie convergence rate.
Computational complexity of the proposed method is K O(N2 C) + K O(N*), where N is network size, C' is sample size, and K is the number of iterations executed. In practice, we require Cs N for reasonable estimation. The first term thus dominates and the computational complexity is approximately the same as that of gradient descent. We can replace the standard SVD operation in the M-step by randomized SVD or Lanczos algorithm. Then, we can further reduce the computational complexity of the second term from ((N3) (r or (( + I) where r < N is the rank of network V NN is the size of network, and is constant usually between 10 and 20.
We also numerically investigate both time and estimation performance of the proposed method in recovering diffusion networks hased on synthetic cascading data in large network settings. We adopt the same criteria used in Section 4 to evaluate estimation performance and use average time per iteration to evaluate time performance. We generate diffusion networks O and y of different sizes in similar way to Section 4.2, fixing Vi to be both low-rank (rank 5) and sparse (edge density 0.01). Given O and Y we generate C' = 50, 000 independent cascade samples based on double mixture model with Exp transmission model and observation window length being T 10.
Table 3: Execution time of different methods under different network sizes (unit: second per iteration).
N 500
N 1000
N 2000
NV 4000
OurAle
2.395(0.001)
8. 8.201(0.019)
33 33.442(0.031)
132, $2.955500040)
NetRate
1,390(0.001)
4.929(0.011)
20.039(0.030)
79.992(0.040)
Table shows that the computation time of the proposed method per iteration is proportional the squared network size N This aligns with the earlier theoretical analysis of computational complexity. In addition, for all network sizes N 500, 1000, 2000, 4000, the ratio between computational time per iteration of the proposed method and NetRate i3] approximately constant 1.65. Since NetRate (31] is well-known for its scalability, this result demonstrates the computational efficiency and scalability of the proposed method.
Table illustrates the estimation performance of the proposed method for large networks. For networks of all sizes N = 500, 1000, 2000, the proposed method outperforms NetRate BI on both transmission rate estimation and network topology recovery. For networks of sizes N = 500, 1000, the proposed method achieves estimation performances comparable to those of smaller networks in Section 4. Additionally, when sample size is fixed, the performance of the proposed method decreases, indicating the need for more samples on larger networks to achieve accurate estimations.
Table 4: Diffusion network estimations from the proposed method under different network sizes.
MdE
MAE.
Acco
Pres
Recy
N 500
0.398(0.001)
0,711(0.11))
0.89440.0033
0.820(0.005)
0.9833(0000)
OurAlg
1000
0.411(0.005)
0.7970.007)
0.779(0.081)
0.666/0.001)
0.939(0.003)
N= 2000
0.457(0.003)
0.814(0.008)
0.571(0.002)
0.523(0.002)
0.630X0.002)
N 500
0.743(0.04)
1.053(0038))
0.078(0.008)
0.818c0.007 0.841(0.002)
NetRate
N 1000
0.755(0.010)
(.94(00019)
0.08600.001)
0.045(0.001)
0.91400.007)
= 2000
0.759(0.003)
1.096/0.002)
0.053(0.000)
0.04450.000)
0.913(0.000)
7 Conclusion
In this paper, we propose novel double network mixture model for heterogeneous cascading process on multi-layer networks. Our method can identify the latent diffusion network complementary to the observed network. Due its convex formulation, our method has both statistical and computational guarantee in terms of estimating diffusion networks. A major future work is to generalize the mixture graph model system with more than two-layer networks and derive the model identification conditions. Extending the proposed method to inference on time-varying networks with time- dependent parameters (O(t,, V (t), (t) or nonparametric transmission models are also interesting directions for future works.
10

---

References
[1] Eytan Adar and Lada A Adamic. Tracking information epidemics blogspace. In The 2005 IEEE/WIC/ACM International Conference on Web Intelligence (WI005) pages 207-214. IEEE, 2005.
[2] Sivaraman Balakrishnan, Martin J. Wainwright, and Bin Yu. Statistical guarantees for the EM algorithm: From population to sample-based analysis. The Annals of Statistics, 45(1):77 120, 2017.
[3] Pierre-Alexandre Balland and David Rigby. The geography of complex knowledge. Economic geography, 93(1):1-23, 2017.
[4] Sushil Bikhchandani, David Hirshleifer, and Ivo Welch. Learning from the behavior of others: Conformity, fads, and informational cascades. Journal of economic perspectives, 12(3):151-170, 1998.
[5] Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal on optimization, 20(4):1956- 1982, 2010.
[6] Venkat Chandrasekaran, Sujay Sanghavi, Pablo A Parrilo, and Alan S Willsky, Rank-sparsity incoherence for matrix decomposition. SIAM Journal on Optimization, 21(2):772-566, 2011.
[7] Diana Crane. Invisible Colleges; Diffusion of Knowledge in Scientific Communities. University of Chicago Press, Chicago,, 1972.
[8] Hadi Daneshmand, Manuel Gomez-Rodriguez,, Le Song, and Bernhard Schoelkopf. Estimating diffusion network structures: Recovery conditions, sample complexity & soft-thresholding algorithm. In International conference on machine learning, pages 793-801. PMLR, 2014.
[9] Nan Du, Song, Manuel Gomez Rodriguez, and Hongyuan Zha. Scalable influence estimation in continuous-tima diffusion networks. Advances in neural information processing systems, 26, 2013.
[10] Nan Du, Le Song, Ming Yuan, and Alex Smola. Learning networks of heterogeneous influence. Advances in neural information processing systems, 25, 2012.
[11] Wenjing Duan, Bin Gu, and Andrew B Whinston. Informational cascades and software adoption on the internet: an empirical investigation. MIS quarterly, pages 23 23-48, 2009.
[12] John RP French. The bases of social power. Studies in social power/Universitty of Michigan Press, 1959.
[13] Noah E Friedkin and Eugene C Johnsen. Social influence network theory: A sociological examination of small group dynamics, volume 33. Cambridge University Press, 2011.
[14] Jerome Friedman, Trevor Hastie, Holger Höfling, and Robert Tibshirani. Pathwise coordinate optimization. The annals of applied statistics, 1(2):302-332, 2007.
[15] P. E. Gill, W. Murray, and M. A. Saunders. Snopt: An sqp algorithm for large-scale constrained optimization. SIAM Review, 47:99-131, 2005.
[16] Philip E. Gill, Walter Murray, Michael A. Saunders, and Elizabeth Wong. User's Guide for SNOPT Version 7.7: Software for Large-Scale Nonlinear Programming. Department of Mathematics, University of California, San Diego, La Jolla, CA 92093-0112, USA, March 2018. CCoM Technical Report 18-1.
[17] Manuel Gomez-Rodriguez, Jure Leskovec, and Andreas Krause. Inferring networks of diffusion and influence. ACM Transactions on Knowledge Discovery from Data (TKDD), 5(4):1-37, 2012.
[18] Manuel Gomez-Rodriguez, Jure Leskovec, and Bernhard Schölkopf. Modeling information propagation with survival theory. In International conference on machine learning, pages 666-674. PMLR, 2013.

---

[19] Manuel Gomez Rodriguez, Jure Leskovec, and Bernhard Schölkopf. Structure and dynamics of information pathways in online media. In Proceedings of the sixth ACM international conference on Web search and data mining, pages 23-32, 2013.
[20] Márton Karsai, Gerardo Iñiguez, Riivo Kikas, Kimmo Kaski, and János Kertész. Local cascades induced global contagion: How heterogeneous thresholds, exogenous effects, and unconcerned behaviour govern online adoption spreading. Scientific reports, 6(1):27178 2016.
[21] Matt J Keeling and Ken TD Eames. Networks and epidemic models. Journal of the royal society interface, 2(4):295-307, 2005.
[22] David Kempe, Jon Kleinberg, and Êva Tardos. Maximizing the spread of influence through social network. In Proceedings of the ninth ACM SIGKDD international conference Knowledge discovery and data mining, pages 137-146, 2003.
[23] Daeyoung Kim and Bruce Lindsay. Empirical identifiability finite mixture models. Annals of the Institute of Statistical Mathematics, 67:745-772, 2015.
[24] Alden S Klovdahl, John J Potterat, Donald E Woodhouse, John B Muth, Stephen Q Muth, and William W Darrow. Social networks and infectious disease: The colorado springs study. Social science & medicine, 38(1):79 88, 1994.
[25]
Jure Leskovec, Lada A Adamic, and Bernardo A Huberman. The dynamics of viral marketing. ACM Transactions on the Web (TWEB), 1(I):5-es, 2007.
[26] David Liben-Nowell and Jon Kleinberg. The link-prediction problem for social networks. journal of the association for information science and technology (2007). Google Scholar Google Scholar Digital Library Digital Library, 2007.
[27] Aditya Krishna Menon and Charles Elkan. Link prediction via matrix factorization. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2011, Athens, Greece, September 5-9, 2011, Proceedings, Part II 22, pages 437-452. Springer, 2011.
[28] Sarah Morrison-Smith and Jaime Ruiz. Challenges and barriers in virtual teams: a literature review. SN Applied Sciences, 2(6): 1096, 2020.
[29] Seth A Myers and Jure Leskovec. On the convexity of latent social network inference. Advances in Neural Information Processing Systems, 2, 2010.
[30] Jeongha Oh, Anjana Susarla, and Yong Tan. Examining the diffusion of user-generated content in online social networks. Available at SSRN 1 1182631, 2008
[31] Manuel Gomez Rodriguez, Jure Leskovee, David Balduzzi, and Bernhard Schölkopf. Un- covering the structure and temporal dynamics of information propagation. Network Science, 2(1):26-65, 2014.
[32] Daniel M Romero, Brendan Meeder, and Jon Kleinberg. Differences in the mechanics of information diffusion across topics: idioms, political hashtags, and complex contagion on twitter. In Proceedings of the 20th international conference on World wide web, pages 695-704, 2011.
[33] Charles R Shipan and Craig Volden. The mechanisms of policy diffusion. American journal of political science, 52(4):840-857 2008
[34] Madeleine Udell and Alex Townsend. Why are big data matrices approximately low rank?
SIAM Journal on Mathematics of Data Science, 1(1):144-160, 2019.
[35] Jacco Wallinga and Peter Teunis. Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures. American Journal of epidemiology, 16006):50--516, 2004.
[36] Liaoruo Wang, Stefano Ermon, and John Hoperoft. Feature-enhanced probabilistic models for diffusion network inference. In Machine Leurning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2012, Bristol, UK, September 24-28, 2012. Proceedings
Part II 23, pages 499-514. Springer, 2012.
12

---

[37] Senzhang Wang, Xia Hu, Philip S Yu, and Zhoujun Li. Mmrate: Inferring multi-aspect diffusion networks with multi-patten cascades. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 1246-1255, 2014.
[38] Zhaoran Wang, Quanquan Gu, Yang Ning, and Han Liu. High dimensional expectation- maximization algorithm: Statistical optimization and asymptotic normality. arXiv preprint arXiv:1412.8729, 2014.
[39] Duncan J Watts. A simple model of global cascades random networks. Proceedings the National Academy of Sciences, 99(9):5766-57)) 2002.
[40] Ming Yu, Varun Gupta, and Mladen Kolar. Estimation of a low-rank topic-based model for information cascades. Journal of Machine Learning Research, 21(71):1-47, 2020.
[41] Gilles Zumbach and Paul Lynch. Heterogeneous volatility cascade in financial markets. Physica
A: Statistical Mechanics and its Applications, 298(3-44:551-555 2001.
13

---

A
Additional details on methodology
A.1
Detailed expression of objectives in M-step
@ ) Q2 v
(2(()), and Q3( (2(+)) can be formulated in detail as:
C
Q )
Σ
log (1 + + Σ
(T;;
C
jusT
02(V/19
Σ
Σ
(1 - +5)
Σ
(1-
(1- - 15)
log (o(T:
C
c=1
jdgz1
N
Qs(T 2
Σ
#s log ग. + (1 +1
#;) log(1 - n: 1
C
where PI, Pu follow the definitions in Section 3.1.
A.2 Optimization algorithm and discussions
We summarize the optimization of the reformulated problem in Section 3,4 via EM algorithm as Algorithm 1.
Algorithm 1 First-order projected EM algorithm
Require: initialization ()) (e (0((), x T()), observed network A, low-rank penalty μ4, learn- ing rate X, and stopping criterion e. while (5+(++1))
()))) |l| y (++ (()(()) (2()) μ y (+)| c do.
E-step: update af P(Z t": ())) via its posterior distribution based on (2(s) for i
1, 1... N, c I.... C.
M-step: decompose Q ( Q1(@) 2 + Q2 v (2())) Q2(T (+())
M.1: Update O via Q1 ( (()):
( max(O (+) + -
( O A,
M.2: Update V via ( ((())):
M.2.1:
aQ2
(=)
M.2.2: perform SVD decomposition on L = UAVT
M.2.3: ((+11) <- max - Xx) V,0)
M.3: Update via | (a+1) 0 1-= 1. N
end while
In Algorithm 1, we utilize projected gradient ascent and proximal gradient ascent to update diffusion networks @ and in M. and M.2 of M-step, where the latter can be implemented via singular value soft-thresholding operation [14 5] For computational efficiency and stability, we utilize the first-order EM algorithm 2 such that the ELBO increased via one-step gradient ascend instead of maximized in M-step. In addition, both the gradients :Q1
@ and eQ2 have closed forms and can be V efficiently calculated. The closed form gradients are provided in A.3.
Likelihood-based parameter tuning The low-rank penalty μ in the above Algorithm can be selected in a data-adapted manner. Specifically, we can first randomly separate the total cascade samples into training subset Ctrain and validation subset Crais and estimate model parameters Gy on Ctrain given a specific μ. Then we calculate the logelikelihood of validation samples as [Caat] ExeCra log P(t(r), ②t,), and select such that
Σ CEConl
log ((t/c),..))
= arg max WEG
|Ceatl
14

---

where G is grid of candidates for penalty values.
Sparsity pursuit θ When the observed network A is not sparse enough the sample size of cascade is relatively small, one can further add / regularizer on θ in the objective function Q1(e (2())), Accordingly, the matrix is updated via an additional coft-thresholding operator in the M.1 step of Algorithm 1. Previous results shows that adopting the 4 regularization can improve the performance and efficiency to recover the diffusion network structure [8 31..
A.3 Closed form gradients in Algorithm 1
Denote (5 15s = S(t 5. Ox) (5 15s = S(t Ej. WR). (H² = H(t Ej. Oj.), and (H² )5 (t; Ey. v ji), where S(- .))) and H( ..)) are survival and hazard functions, respectively.
In addition, we introduce two cascade sample indicator variables (c) and ( (c) for each node pair (j. 2) such that
[1, if < and, 17 & T,
(c)jr
(0, otherwise
and
[1. , TT. and, t2 & T,
1 (c)
(0, otherwise
8Q1 and 2Q₂ over each element in @ and N respectively, can be explicitly formulated Then,
8Q1
Σ (c), T(Ag= =1) (50) j
*
(s
(c), x I(Agr = 2 <IC
#
H
11
20
8Q2 ( kj
Σ
(c)
A 015 já
+ (12)
(c)
A
0
H
)
jr
Vi
H )
Uy jz
The above gradients can be calculated by matrix operations after translating two indicator functions into matrices.
B Additional details on experimental setup
B.1 Benchmark comparison under different network topologies
In this experiment, we generate diffusion networks with different topologies. Specifically, we consider the diffusion network @ and its support A as random network, network with community structure, and scale-free network, For random network, we generate A via Erdos-Renyi model with edge generation probability being 0.01. For the community structure, that contains four equal-sized communities. The generation probability of within-community edge is 0.05 and that of wetween-communine
generate A with stochastic block model
edge is 0.01. For scale-fie free network, Barabási-AlBes Albert model where we set the number of edges to attach from a new node to existing nodes to be 1. After generating support A, set its diagonal be 0 and assign each non-zero edge weight that follows Unif(L 5) to construct θ.
generate A via
We generate the latent diffuusion network as & = V V. where 11 U2 E R We fix the proportion of non-zero edges in V, 12 0.1 and sample weights from fni 2) for these non-zero edges. The generated V has an edge density 0.05 and transmission rates between and 8 on non-zero edges.
In addition, we impute edges on A to ensure that each row and column of o has at least one non-zero element, i.e., the combined network is connected. Given O and V we generate C' = 2, 000 independent cascade samples based on the double mixture model with Exp transmission model and observation window length of T 10.
15

---

B.2 Network recovery under different transmission models
In this experiment, we generate the latent diffusion network v with low-rank structure using the same procedure as in B.1.
Given L and its support B, we then generate O's support A by foreing one overlap with B per column, i.. {(i,j) A, = 1]. and one non-overlap with B per column, i.c. {(i.j) A. = B, 0). For Exp and Pow transmission models, we create θ by sampling transmission rates from f2, 5) for non-zero edges of A. We slightly modify the generation process for Ray model by decreasing the transmission rates on non-zero edges of O and w. We instead sample weights from U0., 0.8) for non-zero edges in 1 12 and weights from f 2) for those in O This modification aims to force significantly large differences in infection time that contribute to large differences in probability of infection among potential parents.
We also add edges in θ to avoid all-zero row or column in Q V to ensure diffusion network connectivity. The resulting @ and L have edge density at 0.01 and 0.05, respectively. Based on O and v, we generate different numbers of cascade samples at C 500, 1000, 1500, 2000 with fix observation window length at T = 10. In the Pow transmission model, we select delay parameter 8
B.3 Network recovery under large network settings
In this experiment, for each network size N 500, 1000, 2000, 4000, we generate the latent diffusion network & following the same procedure as in B.2 but fixing its rank at 5 and edge density at 0.01. Given V,, we also follow the same procedure in B.2 generate O. The resulting has an edge density of 0.001. Based on θ and V we generate C' 50000 independent cascade samples with fixed observation window length of 10.
To cope with memory limits when the network size N and the sample size C are large, in each outer iteration, stream the data in batches of size B: we process one batch time to compute the contributions (likelihood terms, gradients) required by both the E-step and the M-step, accumulate these quantities across all [M/B] batches, and then carry out the parameter update using the ag- gregated totals-thus matching the effect of a full-batch EM update while keeping memory usage bounded. When more memory is available, one can increase the batch size B accordingly for better parallelization on GPU and thus shorter execution time per iteration.
C
Additional numerical results
C.1 Network recovery under different support overlap
In this subsection, we investigate the performance of our proposed method under different levels of overlap between θ and L defined as overlap( O Mt Σi, I(A.y x Bij)/ E, I(A.j B,y), where A, B are the support of θ Y, respectively.
We generate = 1 V following similar procedure in B. except that we sample weights of non- zero edges in 1 and V. from f 2, 1.5) and increase its rank to 30 to better control its overlap degree with @ The edge density of the generated V 0.025 and the transmission rates on its non-zero edges range from 0.05 to 2.
Given v and its support B. we generate different supports A following the same procedure in B.2 while changing their overlaps and non-overlaps per column with B so that overlap( θ \r) takes three levels about 0.1, 0.3, 0.5. Then we create θ by sampling transmission rates from (ni1 0.2) o (1.9, 2.0) for the overlap {{i,j) A.j Bij = 1}, and sampling from (1, for non-overlap support [(. A. = 1, Bi3 03.
We also add edges in with transmission rates from (1, to avoid all-zero row or column in θ J & to ensure diffusion network connectivity. Based on O and L generate different numbers of cascade samples from Exp transmission model, and fix the observation window length at T 10.
Table 5 shows that the diffusion network recovery performance of the proposed method improves as sample sizes increases. Additionally, as overlap level between and V increases, the transmission rates estimations MAE on both θ y and decreases while estimation of V support improves.
16

---

This pattern uncovers trade-off between network differentiation and topology recovery in network mixture model. As the overlap level increases. the joint diffusion complexity on O and & decreases due to the shared diffusion pathways and support information. On the other hand, it becomes more difficult and requires more samples to differentiate @ and & on their overlap component.
Table 5:
Diffusion network estimations from the proposed method under different degrees of overlap- ping between diffusion networks.
overlup O, v
M MALe
MAL
MAE,
Acoe
Proe
Roop
C= 500
0.325(0.014)
0.314(0.013) 0,510(00118)
0.690(0.007)
0.553(0.009)
0.919(0.005)
C= 1000
0.414(0.016) 0.305(
0.300(0.016)
0.793(0.008)
0.974(0.003) (.669(0.010)
0.1
C= 1500
0.307(0.016)
0.391(0.018)
0.297(0.016)
0.711 0.827(0.013)
8)
0.989(0.002)
C= 2000
0.277(0.012)
0.360(0.012)
0.288(0.015)
0.853(0.004)
0.748(0.007)
0.994(0.002)
C= 500
0.577(0.017) 0.393(0.012)
0.764(0.008) 0.435(0.016)
0.647(0.010)
0.933(0.008)
C= 1000
0.388(0.018)
0.519(0.016)
0.45((0.015)
0,852(0.006)
0.756(0.010)
0.977(0.003)
0.3
C= 1500
0.387(0.013)
0A46(00012)) 0.507(.018)
0.887(0.008)
0.804(0.013)
0.989(0.002)
C= 2000
0.366(0.017)
0.489(0.022)
0.452(0.018)
0.908(0.003)
0.835(0.004)
0.994(0.002)
C 500
0,458(0.014)
0.5330.017) 0.618(0.014)
0.810(0.017)
0.712(0.023)
0.940K0.000 006)
C 1000
0.438(0.018)
0.555(0.013)
0.555(0.011)
0.884(0.034)
0.805)00053)
0.975(0.006)
0.5
C = 1500
0.442(0.018)
0,919(0.027) 0.560(0.011) 0.556(0.026)
0,986(0.004) (.863(0
C 2000
0.421(0.019)
0.526(0.019)
0.557(0.011)
0.935(0.030)
0.886(0.052)
0.992(0.004)
C.2 Performance under different time window length
In this subsection, investigate the network recovery performance under different observation window lengths, which ranges in T' {1, 2, 5, 10).
We generate A and T following the same procedure in C.1 and fix their overlap level at about 0.30. We generate C' 1500 cascade samples from transmission models Exp. Pow, Ray.
We illustrate the performance in Figure 5 [] In general, both transmission rate estimations and network topology recovery improve as the observation window becomes T longer since effective diffusion information increases within each cascade sample. Specifically, by Figure the proposed method achieves better MAE of θ and V under Ray (green line with shaded one standard deviation error range) and Pow (yellow line) cascade model. MAE criterion under Exp model is not strictly monotone, but its overall trend matches that of Pow and Ray. Figure 6 shows that accuracy and precision of y also increases as T increases and then remains stable. The results in this subsection suggest to set large T in practice, which is also consistent with our theoretical analysis.
1.0
1.5
1.0
0.6
0.6
0.6
02
0.2
0.2
10
(a) M AAE
(b) M AE
(c) MAE
Figure 5: Parameter estimation under different time window lengths.
C.3 Network recovery without support information
In this subsection, we investigate the performance of the proposed method when support A of diffusion network O is unobserved.
We consider diffusion networks θ and V of size = 100, generated by similar scheme as C. Since the size of network is reduced, rank of V,, V2, where v v V is accordingly reduced to 15. The generated ② and & have edge densities of 0.02 and 0.05, receptively. Non-zero edges of and V have ranges 0.1 to 2 and 0.05 to 2.2, respectively. The level of overlap between the two networks is
17

---

0.8
0.5
0.6
0.4
0.4
02
0.7
10
(a) Accy
(b) Prev
(c) Recally
Figure 6: Topology recovery of latent network estimation under different time window lengths.
fixed at 0.15. We generate cascade samples of different sizes C = 2000, 3000, 4000, 5000 and from all three different transmission models Exp, Pow, and Ray.
Table shows that when support of is not observed, performance of the proposed method on transmission rate estimation still increases number of samples increases for all three cascade generation models. Benchmark values Table 6 are obtained from the proposed method with and and show the relationship between precision and recall
information the support of @. Figure of the estimated θ and V when their entries below some varying thresholds are truncated to O. As sample size increases, the performance of the proposed method network topology recovery increases (shown by the larger area under P-R curves). Although the proposed method without support information of @ cannot achieve as good performance as when the information available. its performance on both transmission rate estimation and network topology recovery still exhibits converging pattern the benchmark performance as sample size increases.
Table 6: Parameter estimation from the proposed method when support of θ unknown.
M
Acceo
M AE
Aone
C 2000
0.601700.228)
0.628(0.013)
0.672(0.027)
0.803(0001.5
C= 3000
0.533(0.038)
0.702(0.016)
0.604(0.030)
0.883(0.007)
C= 4000
0.447(0.107)
0.742(0.053)
0.520(0.089)
0.902(0.005) (000))
5000
0.327(0.113)
0.787(0.047)
0.414(0.095)
0.904(0.024)
Benchmark
0.088(0.005)
0.976(0.003)
0.210(0.008)
0.926(0.006)
C 2000
(665(0.044)
0.699(0.039) 0.573(0.021)
0.601
C= 3000
0.553(0.112)
0.680(0053)
0.592(0.093)
0.699(0.040)
Ray
C= 4000
0.292(0.054)
0.824(0.02))
0.376(0.042)
0.811(0.016)
5000
0.232(0.022)
0.808)
0.326(0.017)
0.837(0.010)
Benchmark
0.067(0.002)
1.000(00000)
0.177(0.002)
0.941(0.023)
2000
0.604(0.004)
0.732(0.026)
0.537(0.051)
0.855(0.010)
3000
0.447(0.051)
0.786(0.043)
0.23(00033))
0.871(0.008)
Pow
C= 4000
0.331(0.062)
0.822(0.025)
0.344(0.048)
0.879(0.010)
C 5000
0.346(0.047)
(0.02)))22
0.341(0.034)
0.86900010)
Benchmark
0.976(0.003) (0.00)
0.175(0.008)
0.916(0.002)
Pertkon
2000
5000
0.4
0/6
000
0.2
10
0.0
0.4
05
Recal
Recwi
Recel
(a) Precision-recaa -ecall, Exp
(b) Precision-recall, Ray
(c) Precision- Pow
Figure 7: Precision-recall plot of +
18

---

1.0
1.0
601)
2000
200
3003
3000
3000
Recall
Recall
Recall
(a) Precision- Exp
(b) Precision-reeall, Ray
(c) Precision-reeeee recall, Pow
Figure 8: Precision-recall plot of yp
D Additional details on method implementation
To give proper initializations to proposed method, estimation result A from NetRate B1. For experiments in Section 4.1, 4.2, C.I, and C.2, we initialize O, = A c A and Vo AQ (1 A), where A the support of ④
We employ new initialization method in C.3 since information on A, the support of θ is missing. We initialize To the truncated SVD of A at rank We then initialize Go the top 7% entries of the residual A o Both r and are tunable hyperparameters, where In the optimization process of C.3, instead of forcing gradient updates of O taking non-zero values only on A, we introduce an additional I regularizer on θ in the loss function to promote its sparsity and update O via additional I, soft-threshldding operator.
r & N and 0 < y < 100.
We also provide short discussion this section on the general strategy to choose the best hyperpa- rameters aex Xu, μ. pe. pv), where À denotes learning rate, i denotes low-rank penalty, p denotes 4 penalty, and the subscript denotes the network associated with.
In general, network size, sample size, or edge density increases, smaller λe, X are required for stable and good performance of the algorithm. In terms of different cascade transmission models, the optimal learning rates for Pow model have the largest magnitude, while those for Ray model have the smallest. Level of overlap, length of observation window, or the rank & does not significantly affect learning rate choices.
For fixed (θ, Vr, π), there exist best threshold A for singular value soft-thresholding operation and best thresholds Xe pe. X pv for 4 soft-threshee loft-thresholding operation. Thus, for fixed (O, V, π). we need adjust penalties inversely proportional learning rates for the best performance of the proposed algorithm. ( penalties should be increased as transmission rates on θ y increase, while should be increased as rank of increases.
E
Proofs for theoretical results
Proposition 1
Proof. To prove the column-wise identification between O.; and & only need to consider node i's direct parents in two networks as R := (j E (1..... .N} O); > 0) and R := {j (1,..., ,N] 3 0} the sets of nodes that directly reach node i on + and v. Denote the
activation times of the direct parents of node as tx. [0, (ि (t, | e R u R
n. (1; O.;) + (1 n.) P1 V.
#. Pa (1: +.) + (1 m. Pr (t: (..), then we have
n
H(ty:
S(t, 0
m)) Σ H(ty:
S(ty: 1e
(15)
(1 म. Σ H(t;: + |[ S(t, 0,,, (1 n.) Σ
H(t; T ()) I
S(ty V j)).
(16)
where the summation and multiplication are over j E R. Ryy with p :- |R R. Given that survival function satisfies S(tr;A))) exp(Ay:h(t, - } for some non-negative function h(.),, then
19

---

the infection likelihood and hazard function are
f(tziAgs) = 23 exp(=>),h(t1 - -t)) H(tzidjj) Aysht (t:-t))
Denote wi (1t) = 7,,' - 11 Stta: Rt), uz(tz) -1- - (1) II St2k: M.).
W1 (()) = T.,'4 - - S(4; kn.) and m2(()) = (1- -2,))1(44 - 11) S(1ks xx)) Then
15 can be re-written
Ewi(t)) + Σ 12(tt) 2 - Σ 11 (t,)) - Σ 12 (2) = 0.
(17)
which should hold for any LR.. Consider a 4p 4p matrix F with each row be- ing (w)).)..... (((p)) (w2())).... w2((), (i())..... ..((p), (2((1), (eal1),...,,(((())) w2((p)] where each row takes different values (R. Based on assumption that there exist o i such that 0 * V then up to permutation between + and L.,, the columns corresponding to different ij will be linear independent. Without loss of generality, we only need to investigate the 4p x 4 submatrix F with each row being [w> ((f)),/u((;)) W1 (tz), (t)())). We prove that WI(ty) = (z w2(ty) (ty(t)) and θj Gjt, Ue Ue by contradiction argument. For two different values ł and ty., we have
w (t,)) (2((t) 11, S(t, ²xx) I1, S(t, vR
(18)
w1 (t,) w2(t,) 11, O₁r) [I, S(t, ,(
Σ(0, - (n))(( - (1) -h(1, - -4)) = 0,
(19)
and we can choose tj # tjJ for those j* such that ,,:; # ,.:.- Then the above equation implies ,²+; L j* which cause contradiction. Therefore, the rank of F is larger than 1.
If there exists b such that x wi (tz) + b x w2(t,) w (t,) hold for any tR.- Then we lett t, for all j R Ryy except j*. Then we have
a x axTh t,+) exp(h(t; - ty+) ( + +b (1 m.) ty#) explh(t; - tf+ )
(20)
-7,h'(t, - (20) exp (h(t: (,) :} aal exp(((tt(() - (0) - ,.+))} } + b exp(h(t; - ( (V,= }} 1,
(21)
(22)
where a' = a x 1/ = b- Notice that at least one of θ ,5²; and ,, + , is not zero and the monotonicity of exp(Xt) in terms of I, the above exponential function equation cannot hold as long as h44;- t) is not constant over t. Therefore, the rank of F can not be 2. Using the same argument above, we can show the rank of 1 can not be 3 as well.
*
If the rank of T is 4. i.e., F is full rank, then I is also full rank by applying above argument on each; j R R.y except j Based on (7). have r[O.] V θ. V
JJ 0, which leads Therefore, wi (l,) wi (4,), w2 w2((1)) (2,) and θj 0js; V V for all j. Finally, by the definition of
to θ V
v 0 and contradicts to assumption 10+ | 11/1 >
wn (1z) and bildew (1z), we can derive म, मं. The statement based on Pt: can be similarly derived.
Proposition 2
Proof. Given Proposition 1, we have the column-wise identification of O.; and V.,,, and only the identification of @ and V up to the column exchange between θ and V. Notice that Proposition 1 guarantee the identification of C = O + V. Then if
max NEA, (AA1))
||
max NEA, (E),
| < 1,
(23)
we have (e) C A₂ V
Ø based on Proposition in G Notice that E Ai( and I E Ag (L), then ② and v are identifiable conditioning on C' is identifiable.
O
20

---

Theorem 3.1
Proof. We first show (s) Ex (p), and Er convex sets for any > 0, p 0 by checking that the convex combination of any two elements in (+), (v), or _ still lies in the set.
For any 1 O E (e(s) and 0 & À & 1,
0 f X e jj + (1 - () < A3 + (1 2) = 31
(24)
(YY + (1- x)=))
A) 1|o, (1 +]] (1 2) ( A]]= =0
E is thus a convex set. Note that >> 0 avoids the trivial case where A zero matrix.
is easy to argue that E is convex set using similar argument to 24
For any L D (p(p) and 0 & X E 1,
| + (1 A) 2211 < X 1 11+ -1
A) 2
(25)
By 25 and similar argument to 124p E (p) is a convex set.
Next, we show Qa(T) ((())) is concave over л. Recall the formulation of Q3(T (2()) in A.1
C
Q2(T T (2()))
C Σ [
log(T) + (1 *²) log(1 π]]
When Ex, both log(T)) and log(1 π) concave
From the above formulation,
QT(T (2()) is concave over by lincarity of concavity.
Without loss of generality, we show Q1 ( is concave on O under the assumption that hazard function Ht t,, A) satisfies àt²(t t', A)/aa22 for t >
Note that this assumption immediately implies the concavity of H(I ()) and the log concavity of S( E A) over X.
Recall the definition of Pr and Py in Section 3.1 we have
log Pr(1 ..)) log( Σ
(t f4 Oj:)
II
S(t; th. Ok:)
(26)
jutyet
log (T: (++)
log( II J4,T
S(T) ty. Oj,)
(27)
Rewrite (26 and 27 in the following form,
log (t O.,) log( Σ Htt tz. ,,))) Σ log S(t, tx, @w))
log Pu O.) = log S(T) ty. O
Then by concavity and monotonicity of log function, composition rule for concavity, and linearity of concavity, both log Py (44 @.,) and log Ptr +..) are concave over O.,.
Concavity of Q ( ( over O immediately follows from composition rule for concavity and the proof for (T(T (T (+))) Concavity of Q2(+) (p(e) directly follows this proof by changing θ
By the above proof, the linear combination Σj A log Py (T; (T; O.y) where nे depends on (2)s instead of + concave over O Therefore, E, (E jug=T 4 log Po (T: O.,)) is concave in terms
of O. In addition, notice that (e. is (N2 N) x (N² N) diagonal block matrix (o(
e
(Q.00) S 0
as 00
0 j * I, which leads to blocks
(N-1)x(N-1)'
1. 1.
N. And the ith block corresponds to the second derivative of
# (1g57) log Pt(t; O.:) + ((ę-)) log (T)
21

---

Based on the above arguments, we only need to show that the strictly concavity for each block in (. ) i.e., Ee _ (1; of log (t; O.)). We then set to be larger enough
@ such that Co D [(fEr)) N and the index of cascade where is infected as 1, Co.
and
Co (1fTr oi log Pt(t, =Ee(
Σ Es 0
Σ
i log (t(t;:
Then we denote Q = c Co o log (r (.;) and follow the argument of Lemma 10 in prove Q is positive definite. Following their notations and a O can write
g to
@ = D(a) + C x(a)[x(aa)'
where D(a) t -. a));a) and a) is a diagonal matrix with [D(15 #$$5 (5:ax) atg (t a)ll" (t (f;a),, and h(t",) Egzzes H( (f;s). In addition, ((a) is a (N 1)-by-Co matrix as
(o) (X())aa) x(t):a) x()))
where each column x(t)aa) = (
Given the assumption that the proba- bilities of being source node P(e) > for v E R where R denotes the set of nodes from which i is reachable via directed path,
follow the same argument in Lemma 10 of [8] to show @. Then (2))) also strictly concave. The strong concavity of Eg E v [2(+))] in terms of can similarly proved.
Q is strictly concave in terms of
Information on computer resources
Numerical experiments in this paper carried out in Google Colab a single Nvidia A100 GPU with 80GB of memory available. Regardless of failed experiments, the entire numerical experiment section took approximately 50 hours of computation time.
In addition, to evaluate the benchmark algorithm ConNle, we used SNOPT (version 7.7) as the underlying constrained optimization solver [15 [6].
22

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims made in the abstract and introduction accurately reflect the paper's contributions and scope.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The paper discusses the limitations of this work in the Conclusion section.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: All assumptions are stated and all proofs are provided in either Methodology section or Appendix.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
The paper fully discloses all the information needed to reproduce the main experimental results. Code for experiments provided supplementary material.
Justification:
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
The paper provides open access to the data and code with sufficient instructions to faithfully reproduce the main experimental results.
Justification:
6. Experimental setting/details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: The paper specifies all the training and test details necessary to understand the results.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
The paper reports standard deviations where appropriate with suitable explana- tions.
Justification:
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed reproduce the experiments?
Question:
23

---

Answer: [Yes]
Justification: The paper provides information on computer resources in the Appendix.
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of
Answer: [Yes]
Justification: The research conducted in the paper conforms, in every respect, with the NeurIPS Code of Ethics
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [Yes]
Justification: In Real cascading data analysis section, we apply method analyze the research topic cascades in social science across U.S. universities and uncover the latent research topic diffusion networks among the top U.S. social science programs.
11. Safeguards
Question: Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: This paper uses open-domain data and code, properly crediting the license and creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: New code accompanying this paper is well documented and submitted as supplementary material.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots. if applicable, as well as details about compensation (if any)?
Answer: [NA]
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Declaration LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
24


Reviews:
[{"review_id": "VpRAkpMBDQ", "reviewer": "Reviewer_3jm1", "summary": "The paper proposes a method for inferring heterogeneous diffusion structure from network cascades. In particular, cascades can follow different diffusion pathways in a multi-layer network topology, where each transmission channel has a different transmission rate. Motivated by this, they propose to infer transmission rates and topologies for both the observed network (for instance a social network) and the latent diffusion network (based on latent social factors), and individual node probabilities that cascades follow these channels. The parameters of the resulting double network mixture model are constrained using sparsity and low-rank constraints, and optimized using the EM algorithm. Guarantees for convexity of the parameter sets are also obtained. Results are reported on three synthetic datasets where the proposed method outperforms existing work, as well as on a dataset of US social science programs where the method recovers meaningful topic diffusion networks.", "questions": "1. Please compare to references [A,B]?\n2. Discuss the applicability of the method to the case where more than two network pathways are involved?\n3. Provide timing analysis to address concerns about scalability?\n4. These are more for clarification- in Proposition 2, line 165, page 5, how does a larger value of the matrix 2-norm indicate lower rank? \nLikewise, in line 166, how does a larger value of the maxnorm indicate a lower sparsity level? Definition of accuracy in line 210 on page 7, how is higher the better?", "limitations": "Scalability of the method is a concern.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The paper is well written and easy to follow along.\n2. The problem of inferring heterogeneous diffusion network structure from cascade data is interesting.\n3. The proposed method's usefulness is demonstrated empirically.\n\nWeaknesses:\n1. Missing citations to more recent work on inferring network structure from cascades, [A, B].\n2. The case where more than two transmission channels may be involved should be discussed as well for a comprehensive treatment.\n3. As admitted by the authors, the method is not scalable beyond networks of a few hundred nodes. This limits its usefulness in a real epidemiological scenario where interventions are required.\n\n[A] Rong, Yu, Qiankun Zhu, and Hong Cheng. \"A model-free approach to infer the diffusion network from event cascade.\" Proceedings of the 25th ACM international on conference on information and knowledge management. 2016.\n[B] Gray, Caitlin, Lewis Mitchell, and Matthew Roughan. \"Bayesian inference of network structure from information cascades.\" IEEE Transactions on Signal and Information Processing over Networks (2020).", "was_revised": true, "final_justification": "The concerns I raised in the review were mostly addressed by the authors' rebuttal. The scalability issue remains, as also pointed out by other reviewers, at the same time, the low rank constraint which causes this scalability problem is arguably necessary for model identification. Given the methodological contributions and rigor of presentation outweigh the scalability downside by some margin, I would recommend a borderline accept rating.", "rebuttal": "We really appreciate the reviewer’s thoughtful feedback, which has significantly improved our paper! In the following response, we discuss and justify the implications of our assumptions, and demonstrate our method's generalizability with extensive experiments.\n\n**Q1**\n\nWe have thoroughly examined both the Bayesian‐MCMC framework [B] and the NPDC algorithm [A]. Although the Bayesian MCMC framework [B] can jointly infer both network topology and transmission‐rate posteriors under general diffusion models including Exp/Pow/Ray models, its computational cost is high. NPDC is a completely model‐free method that only infers network topology from a clustering perspective. Both methods thus support very general diffusion dynamics, yet neither was designed to handle the specific scenario we target—a mixture of propagation pathways spanning two overlapping networks. Due to their similarity in this way and time constraint, we limit our empirical evaluation to an illustrative comparison against NPDC [A].\n\nIn this experiment, we consider network size $N = 100$ and sample $M = 2000$. We try to estimate a true latent network $\\bf{\\Psi}$ whose first 75 singular values decay at a factor of $e^{-1}$ and the others are close to 0. We report the result in the following table.\n\n| Algorithm | $MAE_{\\\\boldsymbol{\\\\Theta}}$ | $MAE_{\\\\boldsymbol{\\\\Psi}}$ | $Acc_{\\\\boldsymbol{\\\\Psi}}$ | $Pre_{\\\\boldsymbol{\\\\Psi}}$ | $Recall_{\\\\boldsymbol{\\\\Psi}}$ |\n|:---------:|:---------------------------:|:---------------------------:|:---------------------------:|:---------------------------:|:-----------------------------:|\n| Our Alg   | 0.155(0.020)               | 0.283(0.041)               | 0.862(0.042)               | 0.761(0.063)               | 1.000(0.001)                 |\n| NPDC      | –                          | –                          | 0.076(0.018)               | 0.040(0.009)               | 0.715(0.248)                 |\n\n**Table 1** Performance of Our Method and NPDC (mean ± std)\n\nResults in the table show that our method outperforms NPDC under this setting. Note that NPDC returns a score for each pair of nodes in the network that indicates how likely there exists an edge between this pair of nodes. To compare the results of NPDC under the specified metrics, these scores are converted to binary values.\n\n**Q2**\n\nIn this paper, we impose sparse and low-rank structure to differentiate diffusion networks. In general, the proposed method can be extended to multiple diffusion networks as long as we have auxiliary information to avoid ambiguity from column permutation across different  networks. For example, we can typically access $p$-dimensional features from each unit as $\\{\\boldsymbol{X} _i = (X _{i1},\\cdots, X _{ip} )\\} _{i=1}^N$. Then we can model $p$ different diffusion networks  $\\boldsymbol{A}^{(k)} \\in \\mathbb{R} _+^{N\\times N},\\; k = 1,\\cdots,p$ as \n\n$\n (\\boldsymbol{A}^{(k)})_{ij} =  f \\big( \\text{sim}( X _{ik}, X _{jk}), \\alpha_k \\big),  \\; 1\\leq i \\neq j \\leq N,\\; k = 1,\\cdots, p,\n$\n\nwhere $\\text{sim}(\\cdot)$ can be feature similarity measurement between unit $i$ and $j$, and $f(\\cdot ; \\alpha_k)$ can be any function modeling the pairwise transmission rates in the $k$-th diffusion network $\\boldsymbol{A}^{(k)}$. Intuitively, $\\boldsymbol{A}^{(k)}$ can be interpreted as pathways via the proximity over the $k$-th features of units. Then the multiple diffusion networks $\\boldsymbol{A}^{(k)},k =1,\\cdots,p$ are identifiable. In addition, when we select model $f(\\cdot ; \\alpha_k)$ is convex in terms of parameters $\\alpha_k$, then the objective function is still convex, and therefore enjoy both the statistical and computational advantages of the original model.\n\n**Q3**\n\nWe agree that EM algorithm is not very scalable in general. However, although our optimization is formulated as EM-type update, we argue that the implicit convex nature of our objective function allows the optimization to be as efficient as gradient descent, which is fundamentally different to classic EM, and therefore scalable to large network setting. \n\nThe first reason for the poor scalability of classic EM algorithm is the posterior inference in the E-step, which may suffer from low computational efficiency due to sampling and numerical integration. However, as we pointed out in section 3.4, the E-step in our algorithm has an explicit form in equation (10), and therefore can be efficiently calculated without using parallel computing. \n\nThe second scalability issue of classic EM is that the optimization in M-step may be expansive. However, we have proved in Theorem 3.1 that the objective function $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ in the M-step is strictly concave, which makes the optimization as efficient as gradient descent. In addition, in each EM-iteration we only perform gradient ascent one time instead of maximizing $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ until convergence, where $\\boldsymbol{\\Omega} = \\\\{\\boldsymbol{\\Theta},\\boldsymbol{\\Psi},\\boldsymbol{\\pi}\\\\}$ are our model parameters. \nIt is been proved in [1] that the above EM-update is equivalent to gradient decent given strictly concave and well-behaved $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$, and therefore it can be proved that our optimization enjoys the geometric convergence rate as \n\n$$\n\\\\|  \\boldsymbol{\\Omega}^{(s)} - \\boldsymbol{\\Omega}^{\\star} \\\\| \\leq \\rho^s \\\\|  \\boldsymbol{\\Omega}^{(0)} - \\boldsymbol{\\Omega}^{\\star}\\\\|,\n$$\n\nfor some $\\rho<1$, where $\\boldsymbol{\\Omega}^{(s)} $ are the estimation from the $s$-th update, and $\\boldsymbol{\\Omega}^{\\star}$ are the true parameter values. The rate equivalently fast to that of gradient descent. \n\nLast but not least, when handling large networks, we can replace the standard SVD operation in the $M$-step by randomized SVD or Lanczos algorithm. Then, we can further reduce the computational cost from $\\mathcal{O}(N^3)$ to $\\mathcal{O}(rN^2)$ where $r<<N$ is the rank of network $\\boldsymbol{\\Psi}$ and $N$ is the size of network (or in practice, $r$ + $k$ where $k$ is usually 10 to 20 for more stable convergence). This is at the same order of matrix multiplication between two low-rank matrix. \\textbf{The overall computational complexity of our algorithm is: $T \\cdot O(N^2 \\cdot M) + T \\cdot O(N^2 \\cdot r)$}, where $T$ is the number of iterations executed and $M$ is sample size.\n\nBesides the above justification of our method's scalability from a theoretic perspective, we also perform numerical experiments to proof it. To handle possible memory constraints when network size is large, we use a mini-batch version of our algorithm such that each update of both E-step and M-step are conducted using $B$ samples where $B$ is batch size. In on complete iteration, we make $\\lceil M/B \\rceil$ updates so that all $M$ samples are fully used. \n\nIn our first experiment, we fix total number of iterations $T = 100$, sample size $M = 5000$, batch size $B = 40$, dimension $r + k = 20$ for randomized SVD, and only change $N = 500, 1000, 2000, 4000$. We report the mean and the standard deviation of execution time in second per iteration for all four different network sizes in the following table.\n\n| Size          | 500         | 1000        | 2000        | 4000        |\n|:-------------:|:-----------:|:-----------:|:-----------:|:-----------:|\n| Time (s/iter) | 0.032(0.002)| 0.952(0.004)| 3.685(0.001)| 14.574(0.006)|\n\n**Table 2** Execution Time Per Iteration under Different Node Sizes (mean ± std).\n\nWe see that when network size $N$ is large enough for GPU kernels to start fully synchronize, the execution time starts to agree with the theoretical computation complexity of our algorithm: As $N$ doubles, execution time per iteration increases by a factor of four. Note that if more memory is available, we can increase the batch size $B$ accordingly for better parallelization on GPU and thus shorter execution time per iteration. Please see response to reviewer 3 for more numerical experiments related to the scalability of our method.\n\n**Q4**\n\nIn Proposition 2 where a lower value for \n\n$\\max _{N \\in \\boldsymbol{\\Lambda _1}(\\boldsymbol{\\Theta}),\\|N\\| _{\\infty} \\leq 1}\\|N\\| _2$ \n\nindicates a higher rank of $N$. Intuitively, for $N \\in \\{ \\boldsymbol{\\Lambda_1}(\\boldsymbol{\\Theta}),\\|N\\|_{\\infty} \\leq 1\\}$, the Frobenius norm $\\| N \\|_F$ is up bounded due to support of $N$ is fixed and magnitude of each element in $N$ is also up bounded. Given the relation that \n\n$\\| N \\|^2_F = \\sum_{i=1}^N \\sigma^2_i$,\n\nif largest $\\sigma_i$, $\\| N \\|_2$ becomes larger, then the number of zero singular values can be larger, i.e., lower rank. \nIf $\\| N \\|_2$ become smaller, then the number of zero singular values increases, i.e., higher rank. Similarly, for term \n\n$\\max_{N \\in \\boldsymbol{\\Lambda _2}(\\boldsymbol{\\Psi}),\\|N\\| _2 \\leq 1}\\|N\\| _{\\infty}$,  consider $N \\in \\{\\boldsymbol{\\Lambda_2}(\\boldsymbol{\\Psi}),\\|N\\|_2 \\leq 1\\}$  basically up bounds the $\\| N \\|^2_F$ since \n\n$\\| N \\| _F \\leq \\sqrt{\\text{rank}}\\| N \\| _2$.  \n\nThen if $\\|N\\|_{\\infty}$ is larger then $N$ has to be sparser.\n\nReferences:\n\n[1] Statistical guarantees for the EM algorithm: From population to sample-based analysis"}, {"review_id": "5c4RIiZF9W", "reviewer": "Reviewer_fDKS", "summary": "This paper proposes a novel double network mixture model for inferring latent diffusion networks in heterogeneous cascade processes. The model represents cascade pathways as a distributional mixture of two networks, capturing complementary diffusion patterns without modeling complex individual states. The authors develop a data-driven optimization approach using EM algorithm with convex relaxations, supported by theoretical guarantees on identifiability and convergence. Experiments on synthetic and real-world data (U.S. university research topic cascades) demonstrate the model’s superiority over baselines in recovering diffusion networks and handling heterogeneity.", "questions": "1. For **Weakness 1**, how might the model be extended to learn or validate a **flexible, nonparametric baseline** hazard $h(\\cdot)$ rather than assuming one of the fixed Exp/Pow/Ray forms?\n\n2.  For **Weakness 2**, can the authors provide detailed scaling benchmarks for networks with **$N>200$, and discuss strategies** (e.g. distributed SVD, stochastic updates) to mitigate the polynomial—or effectively exponential—growth in the EM M-step?\n\n3.  For **Weakness 3**, is it possible to generalize the **binary indicator $Z^c_i\\in\\{0,1\\}$ to a time-varying or multinomial mixing variable $Z^c_i(t)$** so that nodes can switch networks mid-cascade or leverage more than two diffusion channels?\n\n4. For **Weakness 4**, have the authors empirically verified that\n\n$$\n\\frac{\\partial^2}{\\partial\\lambda^2}H\\bigl(t\\mid t',\\lambda\\bigr)\\;\\ge\\;0\n$$\n\nfor **Pow and Ray hazards across the $\\lambda$ values learned in their experiments**, to ensure the M-step convexity conditions truly hold in practice?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\nThis paper demonstrates **theoretical rigor** by ensuring identifiability through structural constraints (sparsity and low-rank assumptions) and offering statistical and computational guarantees via convex optimization. The column-wise network mixture effectively captures **diffusion patterns at an exponential scale**, outperforming single-network baselines in heterogeneous environments. **Extensive experiments** on synthetic datasets with varying topologies and transmission models, as well as on real research-topic cascades, convincingly validate the method’s effectiveness.\n\nWeaknesses:\n1. While the paper illustrates exponential (Exp), power‐law (Pow) and Rayleigh (Ray) models, it really only requires an additive‐risk form\n $$\n       S(t\\mid t_j,\\lambda_{ji})= \\exp\\Bigl(-\\int_{t_j}^{t}\\lambda_{ji}h(s - t_j)\\mathrm{d}s\\Bigr).\n     $$\n\n\nfor some **nonnegative baseline $h(\\cdot)$**. There is no discussion of **choosing or fitting a more flexible, nonparametric $h(\\cdot)$ to capture**, for example, time‐varying or abrupt changes in hazard.\n\n2. The EM algorithm’s computational complexity grows exponentially with network size N; for N > 200, the number of mixture patterns (2^N) may **preclude convergence**, as demonstrated by experiments limited to N = 200. Moreover, enforcing a low-rank constraint on Ψ via nuclear-norm relaxation is suboptimal when the true latent network is high-rank, potentially introducing bias in dense scenarios.\n\n3. The model fixes each node $i$ to one network per cascade ($Z^c_i\\in\\{0,1\\}$), so it **cannot capture mid‐cascade switches or interactions among three or more diffusion channels**.\n\n4. In the paper, Theorem 3.1 requires the hazard function to satisfy the convexity condition(line 192)\n\n$$\n\\frac{\\partial^2}{\\partial\\lambda^2}H\\bigl(t\\mid t',\\lambda\\bigr)\\\\ge0.\n$$\n\nHowever, while the Exp, Pow, and Ray models typically meet this requirement, the authors do not verify it over the actual $\\lambda$ values learned in their experiments. If this convexity assumption is violated, the M-step convergence guarantees may no longer hold. **I recommend numerically evaluating**\n\n$$\n\\frac{\\partial^2 H(t\\mid t',\\lambda)}{\\partial\\lambda^2}\n$$\n\nfor the Pow and Ray hazard functions across the trained $\\lambda$ range to confirm that the condition indeed holds in practice.", "was_revised": true, "final_justification": "The authors provided reasonable responses to the questions I raised and supplemented them with corresponding experiments. In their response, they addressed and clarified **Theorem 3.1’s assumption of the convexity of the hazard function** and the **computational complexity of the EM algorithm**, and further noted that in the picked experiments in Section 4.1, they numerically verified that\n\n$$\n\\frac{\\partial^2 H(t \\mid t^{\\prime}, \\lambda)}{\\partial \\lambda^2}\n$$\n\nfor the learned $\\lambda$ across the Pow/Ray models in these experiments is 0. This resolves the issues I mentioned in my review, and therefore, I am willing to raise my score to 4.", "rebuttal": "We really appreciate the reviewer’s thoughtful feedback, which has significantly improved our paper! In the following response, we discuss and justify the implications of our assumptions. We also conduct extensive experiments to demonstrate the generalizability of our method.\n\n**Q1**\n\nThe transmission model in this paper can be readily generalized to nonparametric model such as kernel function following [3]. Then we will add additional layer of optimization in the M-step for updating the nonparametric model. However, it is found that Exp/Pow/Ray forms can efficiently capture the activation time lag in many real cascade [2]. In addition, the focus of our paper is to model the complex cascade mixing pattern and recover the underlying diffusion networks, instead of fitting the observed activation time lags. Based on the above reasons, we use Exp/Pow/Ray model in this paper, and leave the extension to the nonparametric transmission model as a future study.  \n\n**Q2**\n\nWe agree that EM algorithm is not very scalable in general. However, although our optimization is formulated as EM-type update, we argue that the implicit convex nature of our objective function allows the optimization to be as efficient as gradient descent, which is fundamentally different to classic EM, and therefore scalable to large network setting. \n\nThe first reason for the poor scalability of classic EM algorithm is the posterior inference in the E-step, which may suffer from low computational efficiency due to sampling and numerical integration. However, as we pointed out in section 3.4, the E-step in our algorithm has an explicit form in equation (10), and therefore can be efficiently calculated without using parallel computing. \n\nThe second scalability issue of classic EM is that the optimization in M-step may be expansive. However, we have proved in Theorem 3.1 that the objective function $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ in the M-step is strictly concave, which makes the optimization as efficient as gradient descent. In addition, in each EM-iteration we only perform gradient ascent one time instead of maximizing $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ until convergence, where $\\boldsymbol{\\Omega} = \\\\{\\boldsymbol{\\Theta},\\boldsymbol{\\Psi},\\boldsymbol{\\pi}\\\\}$ are our model parameters. \nIt is been proved in [1] that the above EM-update is equivalent to gradient decent given strictly concave and well-behaved $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$, and therefore it can be proved that our optimization enjoys the geometric convergence rate as \n\n$$\n\\\\|  \\boldsymbol{\\Omega}^{(s)} - \\boldsymbol{\\Omega}^{\\star} \\\\| \\leq \\rho^s \\\\|  \\boldsymbol{\\Omega}^{(0)} - \\boldsymbol{\\Omega}^{\\star}\\\\|,\n$$\n\nfor some $\\rho<1$, where $\\boldsymbol{\\Omega}^{(s)} $ are the estimation from the $s$-th update, and $\\boldsymbol{\\Omega}^{\\star}$ are the true parameter values. The rate equivalently fast to that of gradient descent. \n\nLast but not least, when handling large networks, we can replace the standard SVD operation in the $M$-step by randomized SVD or Lanczos algorithm. Then, we can further reduce the computational cost from $\\mathcal{O}(N^3)$ to $\\mathcal{O}(rN^2)$ where $r<<N$ is the rank of network $\\boldsymbol{\\Psi}$ and $N$ is the size of network (or in practice, $r$ + $k$ where $k$ is usually 10 to 20 for more stable convergence). This is at the same order of matrix multiplication between two low-rank matrix. \\textbf{The overall computational complexity of our algorithm is: $T \\cdot O(N^2 \\cdot M) + T \\cdot O(N^2 \\cdot r)$}, where $T$ is the number of iterations executed and $M$ is sample size.\n\nBesides the above justification of our method's scalability from a theoretic perspective, we also perform numerical experiments to proof it. To handle possible memory constraints when network size is large, we use a mini-batch version of our algorithm such that each update of both E-step and M-step are conducted using $B$ samples where $B$ is batch size. In on complete iteration, we make $\\lceil M/B \\rceil$ updates so that all $M$ samples are fully used. \n\nIn our first experiment, we fix total number of iterations $T = 100$, sample size $M = 5000$, batch size $B = 40$, dimension $r + k = 20$ for randomized SVD, and only change $N = 500, 1000, 2000, 4000$. We report the mean and the standard deviation of execution time in second per iteration for all four different network sizes in the following table.\n\n| Size          | 500         | 1000        | 2000        | 4000        |\n|:-------------:|:-----------:|:-----------:|:-----------:|:-----------:|\n| Time (s/iter) | 0.032(0.002)| 0.952(0.004)| 3.685(0.001)| 14.574(0.006)|\n\n**Table 1** Execution Time Per Iteration under Different Node Sizes (mean ± std).\n\nWe see that when network size $N$ is large enough for GPU kernels to start fully synchronize, the execution time starts to agree with the theoretical computation complexity of our algorithm: As $N$ doubles, execution time per iteration increases by a factor of four. Note that if more memory is available, we can increase the batch size $B$ accordingly for better parallelization on GPU and thus shorter execution time per iteration. Please see response to reviewer 3 for more numerical experiments related to the scalability of our method.\n\n**Q3**\n\nIn this paper, we impose sparse and low-rank structure to differentiate diffusion networks. In general, the proposed method can be extended to multiple diffusion networks as long as we have auxiliary information to avoid ambiguity from column permutation across different  networks. For example, we can typically access $p$-dimensional features from each unit as $\\\\{\\boldsymbol{X}_ i = (X_ {i1},\\cdots, X_ {ip} )\\\\}_ {i=1}^N$. Then we can model $p$ different diffusion networks $\\boldsymbol{A}^{(k)} \\in \\mathbb{R}_+^{N\\times N}\\, k = 1,\\cdots,p$ as\n\n$$ (\\boldsymbol{A}^{(k)})_ {ij} = f\\big(\\text{sim}( X_ {ik}, X_ {jk}), \\alpha_k \\big),\\; 1\\leq i \\neq j \\leq N,\\; k = 1,\\cdots, p, $$\n\nwhere sim$(\\cdot)$ can be feature similarity measurement between unit $i$ and $j$, and $f(\\cdot ; \\alpha_k)$ can be any function modeling the pairwise transmission rates in the $k$-th diffusion network $\\boldsymbol{A}^{(k)}$. Intuitively, $\\boldsymbol{A}^{(k)}$ can be interpreted as pathways via the proximity over the $k$-th features of units. Then the multiple diffusion networks $\\boldsymbol{A}^{(k)},k =1,\\cdots,p$ are identifiable. In addition, when we select model $f(\\cdot ; \\alpha_k)$ is convex in terms of parameters $\\alpha_k$, then the objective function is still convex, and therefore enjoy both the statistical and computational advantages of the original model. \n\nThe proposed method is also straightforward to extend to time-varying inference by introducing time-varying weights $w_c(t)$ that penalize old cascades, i.e.,the older a cascade $c$, the smaller its $w_c(t)$. Then the corresponding time-varying diffusion networks and mixing parameters become $\\boldsymbol{\\Theta}(t),\\boldsymbol{\\Psi}(t), \\boldsymbol{\\pi}(t)$, both of which smoothly changes over time and that\nrecent cascades have higher importance in determining current network structure than old cascades.\n\n\n**Q4**\n\nSince we only use Exp/Pow/Ray model as explained in our response to the first question, $\\partial^2H(t |t^{\\prime}, \\lambda)/\\partial\\lambda^2$ is 0 and does not depend on $t$ or $\\lambda$ for all three models. We have also picked experiments in Section 4.1 and numerically verified that $\\partial^2H(t |t^{\\prime}, \\lambda)/\\partial\\lambda^2$ for learned $\\lambda$ across Pow/Ray models in these experiments are 0.\n\nReferences:\n\n[1] Statistical guarantees for the EM algorithm: From population to sample-based analysis\n\n[2] Uncovering the structure and temporal dynamics of information propagation\n\n[3] Learning networks of heterogeneous influence"}, {"review_id": "2bnoQoTP8L", "reviewer": "Reviewer_q7jT", "summary": "This paper proposes a novel double network mixture model to infer hidden diffusion networks when information spread in terms of cascading behavior is highly heterogenous. Authors provided theoretical basis and conducted several experiments to support their idea.", "questions": "Questions:\n1. As mentioned before, I am quite concerned about the scalability of the proposed method. Would it be easy to extend the method so that it could be applied to networks with at least thousands of nodes if not millions. \n2. Is it possible to extend the method so that it instead learns structures rather than making hard assumptions such as low-rank latent networks?\n3. Transmission rates are assumed over static networks. Can they be generalized to dynamic setting as seen in real world where nodes’ behavior may change over time?", "limitations": "yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. The core idea of the paper, i.e., double network mixture model and distributional mixture representation of cascades is novel and it allows the model to capture complex diffusion patterns between observed and latent channels. The paper reads well and the topic is interesting.\n2. Authors have provided rigorous theoretical basis to back the proposed method that guarantees the algorithm will eventually converge.\n\nWeaknesses:\n1. As also acknowledged by the authors, due to the presence of EM step, the proposed method is computationally expensive and this limits scalability of their algorithm to larger networks (N=200 is significantly smaller than most of the synthetic and real-world networks.\n2. Several strong assumptions make the paper less reliable IMO. Examples include but not limited to imposing low-rank assumption on the latent network.\n3. Accurate infection times are usually hard to obtain in real-world scenarios and the method seems to require precise infection times.", "was_revised": true, "final_justification": "I have read the rebuttal and seen their new experiments that have convinced me their work is ready for publication.", "rebuttal": "We really appreciate the reviewer’s thoughtful feedback, which has significantly improved our paper! In the following response, we discuss and justify the implications of our assumptions, and demonstrate our method's generalizability with extensive experiments.\n\n**Q1**\n\nWe agree with the reviewer that EM algorithm is not very scalable in general. However, although our optimization is formulated as EM-type update, we argue that the implicit convex nature of our objective function allows optimization to be as efficient as gradient descent, which is fundamentally different to classic EM, and therefore scalable to large network setting. \n\nThe first reason for the poor scalability of classic EM algorithm is the posterior inference in the E-step, which may suffer from low computational efficiency due to sampling and numerical integration. However, as we illustrate in section 3.4, E-step in our algorithm has an explicit form in equation (10), and therefore can be efficiently calculated.\n\nThe second scalability issue of classic EM is that the optimization in M-step may be expansive. However, we have proved in Theorem 3.1 that the objective function $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ in the M-step is strictly concave, which makes the optimization as efficient as gradient descent. In addition, in each EM-iteration we only perform gradient ascent once instead of maximizing $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ until convergence, where $\\boldsymbol{\\Omega} = \\\\{\\boldsymbol{\\Theta},\\boldsymbol{\\Psi},\\boldsymbol{\\pi}\\\\}$ are our model parameters. \nIt is been proved in [1] that the above EM-update is equivalent to gradient decent given strictly concave and well-behaved $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$, and therefore it can be proved that our optimization enjoys the geometric convergence rate as \n\\begin{align*}\n    \\\\|  \\boldsymbol{\\Omega}^{(s)} - \\boldsymbol{\\Omega}^{\\star} \\\\| \\leq \\rho^s \\\\|  \\boldsymbol{\\Omega}^{(0)} - \\boldsymbol{\\Omega}^{\\star}\\\\|,\n\\end{align*}\nfor some $\\rho<1$, where $\\boldsymbol{\\Omega}^{(s)} $ are the estimation from the $s$-th update, and $\\boldsymbol{\\Omega}^{\\star}$ are the true parameter values. The rate equivalently fast to that of gradient descent. \n\nLast but not least, when handling large networks, we can replace the standard SVD operation in the $M$-step by randomized SVD or Lanczos algorithm. Then, we can further reduce the computational cost from $\\mathcal{O}(N^3)$ to $\\mathcal{O}(rN^2)$ where $r<<N$ is the rank of network $\\boldsymbol{\\Psi}$ and $N$ is the size of network (or in practice, $r$ + $k$ where $k$ is usually 10 to 20 for more stable convergence). This is at the same order of matrix multiplication between two low-rank matrix. The overall computational complexity of our algorithm is: $T \\cdot O(N^2 \\cdot M) + T \\cdot O(N^2 \\cdot r)$, where $T$ is the number of iterations executed and $M$ is sample size.\n\nWe also perform numerical experiments to justify the above theoretical perspective. To handle possible memory constraints when network size is large, we use a mini-batch version of our algorithm such that each update of both E-step and M-step are conducted using $B$ samples where $B$ is batch size. In on complete iteration, we make $\\lceil M/B \\rceil$ updates so that all $M$ samples are fully used. \n\nIn our first experiment, we fix total number of iterations $T = 100$, sample size $M = 5000$, batch size $B = 40$, dimension $r + k = 20$ for randomized SVD, and only change $N = 500, 1000, 2000, 4000$. We report the mean and the standard deviation of execution time in second per iteration for all four different network sizes in the following table.\n\n| Size          | 500         | 1000        | 2000        | 4000        |\n|:-------------:|:-----------:|:-----------:|:-----------:|:-----------:|\n| Time (s/iter) | 0.032(0.002)| 0.952(0.004)| 3.685(0.001)| 14.574(0.006)|\n\n**Table 1** Execution Time Per Iteration under Different Node Sizes (mean ± std).\n\nIn our second experiment, we fix network size $N = 1000$, sample $M = 20000$, and batch size $B = 300$. We report the execution time and performance of our algorithm against NetRate[2] when $T = 50, 75, 100$ to demonstrate the performance and scalability of our algorithm, and to give a rough estimate of how much iterations are required for ideal performance.\n\n| Algorithm | Metric                             | $T = 50$       | $T = 75$       | $T = 100$      |\n|:---------:|:----------------------------------:|:--------------:|:--------------:|:--------------:|\n| Our Alg   | Time (s)                           | 181.389(0.563) | 270.689(0.491) | 373.010(0.718) |\n|           | $MAE_{\\\\boldsymbol{\\\\Theta}}$       | 0.222(0.001)   | 0.212(0.001)   | 0.207(0.001)   |\n|           | $MAE_{\\\\boldsymbol{\\\\Psi}}$         | 0.690(0.002)   | 0.702(0.001)   | 0.693(0.001)   |\n|           | $MAE_{\\\\pi}$                       | 0.061(0.001)   | 0.062(0.001)   | 0.063(0.002)   |\n|           | $Acc_{\\\\boldsymbol{\\\\Psi}}$         | 0.669(0.005)   | 0.730(0.005)   | 0.773(0.003)   |\n| NetRate   | Time (s)                           | 104.849(0.010) | 156.660(0.062) | 208.546(0.109) |\n|           | $MAE_{\\\\boldsymbol{\\\\Theta}}$       | 0.870(0.002)   | 0.876(0.001)   | 0.880(0.002)   |\n|           | $MAE_{\\\\boldsymbol{\\\\Psi}}$         | 0.841(0.001)   | 0.842(0.001)   | 0.843(0.001)   |\n|           | $Acc_{\\\\boldsymbol{\\\\Psi}}$         | 0.501(0.013)   | 0.502(0.015)   | 0.493(0.018)   |\n\n**Table 2** Performance of Our Method and NetRate under Different Total Iterations (mean ± std).\n\nBy the table, the ratio between our algorithm's execution time and that of NetRate [2] is a constant independent of $N$. This demonstrates the scalability of our algorithm since NetRate [2] is can be scaled efficiently to large networks. Also, our algorithm achieves better performance of that of NetRate [2].\n\n**Q2**\n\nFirst, we argue that low-rank assumption is not hard assumption since it is an intrinsic properties of many real-world networks. one can interpret $\\boldsymbol{\\Psi}$ as the diffusion over latent relations among the units that are not captured by the structural network $\\boldsymbol{\\Theta}$. The magnitude of $\\boldsymbol{\\Psi}_{ij}$ reflects the distance between unit $i$ and $j$ in terms of their latent factors. It has been found that large-scale networks typically has or can be approximated by low-rank structure [4][5][6] due to that high-dimensional latent factors can be always embeded into a low dimensional latent space that preserve original distances [4]. \n\nSecond, our method does not require network $\\boldsymbol{\\Psi}$ to be exact low-rank. For example, our method works even when all singular values are non-zero as long as the singular values decays quickly. In many real-world networks, the singular values indeed decays at exponential rate [3]. To demonstrate this point, we perform numerical experiments to show our method still recovers diffusion network violating low-rank assumption. \n\nIn this experiment, we consider network size $N = 1000$, sample $M = 40000$, and batch size $B = 300$. We try to estimate a true latent network $\\bf{\\Psi}$ whose first 500 singular values decay at a factor of $e^{-1}$ and the others are close to 0. We report the result in the following table.\n\n| Algorithm | $MAE_{\\\\boldsymbol{\\\\Theta}}$ | $MAE_{\\\\boldsymbol{\\\\Psi}}$ | $Acc_{\\\\boldsymbol{\\\\Psi}}$ |\n|:---------:|:-----------------------------:|:---------------------------:|:---------------------------:|\n| Our Alg   | 0.259(0.001)                  | 0.446(0.002)                | 0.903(0.010)               |\n| NetRate   | 0.768(0.007)                  | 0.769(0.002)                | 0.639(0.003)               |\n\n**Table 3** Performance of Our Method and NetRate when Low Rank Assumption is Violated (mean ± std).\n\nFinally, we impose sparse structure mainly for identifying different diffusion networks. The assumption can be droped as long as we have auxiliary information to avoid ambiguity from column permutation across different  networks. For example, we can typically access $p$-dimensional features from each unit as $\\\\{\\boldsymbol{X}_ i = (X_ {i1},\\cdots, X_ {ip} )\\\\}_ {i=1}^N$\n. Then we can model $p$ different diffusion networks $\\boldsymbol{A}^{(k)} \\in \\mathbb{R}_+^{N\\times N}\\, k = 1,\\cdots,p$ as\n\n$$ (\\boldsymbol{A}^{(k)})_ {ij} = f\\big(\\text{sim}( X_ {ik}, X_ {jk}), \\alpha_k \\big),\\; 1\\leq i \\neq j \\leq N,\\; k = 1,\\cdots, p, $$\n\nwhere sim$(\\cdot)$ can be feature similarity measurement between unit $i$ and $j$, and $f(\\cdot ; \\alpha_k)$ can be any function modeling the pairwise transmission rates in the $k$-th diffusion network $\\boldsymbol{A}^{(k)}$. Intuitively, $\\boldsymbol{A}^{(k)}$ can be interpreted as pathways via the proximity over the $k$-th features of units. Then the multiple diffusion networks $\\boldsymbol{A}^{(k)},k =1,\\cdots,p$ are identifiable. In addition, when we select model $f(\\cdot ; \\alpha_k)$ is convex in terms of parameters $\\alpha_k$, then the objective function is still convex, and therefore enjoy both the statistical and computational advantages of the original model.   \n\n**Q3**\n\nThe proposed method is also straightforward to extend to time-varying inference by introducing time-varying weights $w_c(t)$ that penalize old cascades, i.e.,the older a cascade $c$, the smaller its $w_c(t)$. Then the corresponding time-varying diffusion networks and mixing parameters become $\\boldsymbol{\\Theta}(t),\\boldsymbol{\\Psi}(t), \\boldsymbol{\\pi}(t)$, both of which smoothly changes over time and that recent cascades have higher importance in determining current network structure than old cascades.\n\nReferences:\n \n[1] Statistical guarantees for the EM algorithm: From population to sample-based analysis\n\n[2] Uncovering the Temporal Dynamics of Diffusion Networks\n\n[3] The low-rank hypothesis of complex systems\n\n[4] Why are big data matrices approximately low rank?\n\n[5] The link-prediction problem for social networks.\n\n[6] Link prediction via matrix factorization"}, {"review_id": "okXzLxLn4g", "reviewer": "Reviewer_LD4Y", "summary": "This paper presents a novel double network mixture model to infer latent diffusion structures from observed cascade data, in settings where cascading behavior is governed by heterogeneous and possibly unobserved pathways. The key idea is to model the diffusion process as a mixture over two networks—one possibly observed and one latent—by introducing node-wise, cascade-specific Bernoulli indicators. The authors propose an EM-based estimation framework with convex relaxation techniques (nuclear norm), derive identifiability conditions, and provide both statistical guarantees and empirical validation through simulations and real-world sociology publication data.", "questions": "1. **Posterior variance of $\\hat{\\pi}_i^c$:** The estimation of $\\hat{\\pi}_i^c$ is treated as a point estimate (Eq. 10). Can the authors elaborate on the uncertainty of these estimates, or propose a way to obtain variance/confidence intervals?\n\n2. **Scalability:** The authors acknowledge EM’s scalability limitations in the discussion. It would help if they can give a rough estimate of runtime or memory complexity as N grows. Also: is it possible to parallelize E-step across cascades?\n\n3. **Interpretation of π vector:** In real data, what does a high π\\_i mean in practice? Would be nice to provide a short discussion of how π correlates with observable university features (e.g., ranking, location).\n\n4. **ρ selection:** How is ρ (the nuclear norm constraint) chosen in practice? Is there a validation process, or is it fixed?\n\n5. **Real-world validation:** In the sociology case study, results are interesting but more qualitative validation (e.g., overlap with known collaborations or conference co-membership) would make the latent network interpretation more convincing.", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths: The paper tackles a relevant and underexplored problem: inferring heterogeneous diffusion structures when cascades spread through multiple unobserved mechanisms. The methodology is sound, well-formulated mathematically, and supported by solid theoretical guarantees (identifiability, convexity of optimization). The empirical evaluation is quite extensive, including different network topologies and transmission models. The clarity is generally clear, although not always very fluid to read. A few symbols (e.g., $\\Lambda$) are overloaded, and figure captions could be more descriptive. The contribution is significant to the field of network inference and cascade modeling.  Application to sociological topic diffusion also illustrates a valuable real-world use case.\n\nNo major weaknesses to me at this point. See questions below.", "was_revised": true, "final_justification": "I have read the rebuttal and seen their new experiments that have convinced me their work is ready for publication.", "rebuttal": "We really appreciate the reviewer’s thoughtful feedback, which has significantly improved our paper! In the following response, we discuss and justify the implications of our assumptions, and demonstrate our method's generalizability with extensive experiments.\n\n**Q1**\n\nUnlike model parameters $\\bf{\\Theta}$ and $\\bf{\\Psi}$ and $\\boldsymbol{\\pi}$, $\\hat{\\pi}_i^c$ indicates the propability of unit $i$ choosing diffusion netork $\\boldsymbol{\\Theta}$ at the $c$-th cascade, which is treated as part of the missing data in EM algorithm rather than parameters we truly wish to estimate. In other words, we do not have replicate sample to estimate the variance $\\hat{\\pi}_i^c$ since it changes over both unit sample $i$ and cascade sample $c$. On the other hand, we can estimate the uncertainty of $\\hat{\\pi}_i$, which is provided in the tables as standard deviations. \n\n**Q2**\n\nWe agree that EM algorithm is not very scalable in general. However, although our optimization is formulated as EM-type update, we argue that the implicit convex nature of our objective function allows the optimization to be as efficient as gradient descent, which is fundamentally different to classic EM, and therefore scalable to large network setting. \n\nThe first reason for the poor scalability of classic EM algorithm is the posterior inference in the E-step, which may suffer from low computational efficiency due to sampling and numerical integration. However, as we pointed out in section 3.4, the E-step in our algorithm has an explicit form in equation (10), and therefore can be efficiently calculated without using parallel computing. \n\nThe second scalability issue of classic EM is that the optimization in M-step may be expansive. However, we have proved in Theorem 3.1 that the objective function $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ in the M-step is strictly concave, which makes the optimization as efficient as gradient descent. In addition, in each EM-iteration we only perform gradient ascent one time instead of maximizing $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$ until convergence, where $\\boldsymbol{\\Omega} = \\\\{\\boldsymbol{\\Theta},\\boldsymbol{\\Psi},\\boldsymbol{\\pi}\\\\}$ are our model parameters. \nIt is been proved in [1] that the above EM-update is equivalent to gradient decent given strictly concave and well-behaved $\\boldsymbol{Q}(\\boldsymbol{\\Omega})$, and therefore it can be proved that our optimization enjoys the geometric convergence rate as \n\n$$\n\\\\|  \\boldsymbol{\\Omega}^{(s)} - \\boldsymbol{\\Omega}^{\\star} \\\\| \\leq \\rho^s \\\\|  \\boldsymbol{\\Omega}^{(0)} - \\boldsymbol{\\Omega}^{\\star}\\\\|,\n$$\n\nfor some $\\rho<1$, where $\\boldsymbol{\\Omega}^{(s)} $ are the estimation from the $s$-th update, and $\\boldsymbol{\\Omega}^{\\star}$ are the true parameter values. The rate equivalently fast to that of gradient descent. \n\nLast but not least, when handling large networks, we can replace the standard SVD operation in the $M$-step by randomized SVD or Lanczos algorithm. Then, we can further reduce the computational cost from $\\mathcal{O}(N^3)$ to $\\mathcal{O}(rN^2)$ where $r<<N$ is the rank of network $\\boldsymbol{\\Psi}$ and $N$ is the size of network (or in practice, $r$ + $k$ where $k$ is usually 10 to 20 for more stable convergence). This is at the same order of matrix multiplication between two low-rank matrix. \\textbf{The overall computational complexity of our algorithm is: $T \\cdot O(N^2 \\cdot M) + T \\cdot O(N^2 \\cdot r)$}, where $T$ is the number of iterations executed and $M$ is sample size.\n\nBesides the above justification of our method's scalability from a theoretic perspective, we also perform numerical experiments to proof it. To handle possible memory constraints when network size is large, we use a mini-batch version of our algorithm such that each update of both E-step and M-step are conducted using $B$ samples where $B$ is batch size. In on complete iteration, we make $\\lceil M/B \\rceil$ updates so that all $M$ samples are fully used. \n\nIn our first experiment, we fix total number of iterations $T = 100$, sample size $M = 5000$, batch size $B = 40$, dimension $r + k = 20$ for randomized SVD, and only change $N = 500, 1000, 2000, 4000$. We report the mean and the standard deviation of execution time in second per iteration for all four different network sizes in the following table.\n\n| Size          | 500         | 1000        | 2000        | 4000        |\n|:-------------:|:-----------:|:-----------:|:-----------:|:-----------:|\n| Time (s/iter) | 0.032(0.002)| 0.952(0.004)| 3.685(0.001)| 14.574(0.006)|\n\n**Table 1** Execution Time Per Iteration under Different Node Sizes (mean ± std).\n\nWe see that when network size $N$ is large enough for GPU kernels to start fully synchronize, the execution time starts to agree with the theoretical computation complexity of our algorithm: As $N$ doubles, execution time per iteration increases by a factor of four. Note that if more memory is available, we can increase the batch size $B$ accordingly for better parallelization on GPU and thus shorter execution time per iteration. Please see response to reviewer 3 for more numerical experiments related to the scalability of our method. \n\n**Q3**\n\nIf we code $Z=1$ to mean diffusion via network $\\boldsymbol{\\Theta}$, and $Z=0$ for network $\\boldsymbol{\\Psi}$, then high $\\pi_i$ indicates the $i$-th unit more prefers to interact with cascade via network $\\boldsymbol{\\Theta}$ and vice versa. In our application, a high $\\pi_i$ indicates that the institution $i$ have higher probability to communicate research topic via latent academic network instead of geographic network. We have studied the relation between $\\pi_i$ with university rankings and found that high-ranked universities have higher tendency of exchanging idea via latent network. The average $\\pi_i$ among the 104 universities is $0.601$, indicating that research topics diffuse via latent network more frequently than via geographic network. Due to the space limit, we will add the additional analysis in Supplementary. \n\n**Q4**\n\nWe choose the nuclear norm constraint based on cross-validation in practice. Here is also a short discussion on the general strategy to choose the best hyperparameters \n\n$[\\lambda_{\\boldsymbol{\\Theta}}, \\lambda_{\\boldsymbol{\\Psi}}, \\mu, \\rho_{\\boldsymbol{\\Theta}}, \\rho_{\\boldsymbol{\\Psi}}]$, \n\nwhere $\\lambda$ denotes learning rate, $\\mu$ denotes low-rank penalty, $\\rho$ denotes $l_1$ penalty, and the subscript denotes the network associated with. In general, as network size, sample size, or edge density increases, smaller \n\n$\\lambda_{\\boldsymbol{\\Theta}}, \\lambda_{\\boldsymbol{\\Psi}}$ \n\nare required for stable and good performance of the algorithm. \n\nFor fixed $\\{\\boldsymbol{\\Theta}, \\boldsymbol{\\Psi}, \\boldsymbol{\\pi}\\}$,  there exist best threshold $\\lambda_{\\boldsymbol{\\Psi}} \\cdot\\mu$ for singular value soft-thresholding operation and best thresholds \n\n$\\lambda_{\\boldsymbol{\\Theta}} \\cdot \\rho_{\\boldsymbol{\\Theta}}$, $\\lambda_{\\boldsymbol{\\Psi}} \\cdot \\rho_{\\boldsymbol{\\Psi}}$ \n\nfor $l_1$ soft-thresholding operation. Thus, for fixed $\\{\\boldsymbol{\\Theta}, \\boldsymbol{\\Psi}, \\boldsymbol{\\pi}\\}$, we need to adjust penalties inversely proportional to learning rates for the best performance of the proposed algorithm. $l_1$ penalties should be increased as transmission rates on $\\boldsymbol{\\Theta}, \\boldsymbol{\\Psi}$ increase, while $\\mu$ should be increased as rank of $\\boldsymbol{\\Psi}$ increases.\n\n**Q5**\n\nWe agree with the reviewer on this point in general. Notice that the dataset involves information and authors from $29,725$ unique articles from 23 top generalist sociology journal. Extracting all the authors and preprocessing their collaborations relations and conference co-membership can be intensively time-consuming. Therefore, we focus more on introducing and illustrating a new cascade analysis method in this paper, and will add more qualitative validation as future works.\n\nReferences:\n \n[1] Statistical guarantees for the EM algorithm: From population to sample-based analysis"}]

Author rebuttal:

Meta-review: The paper proposes a double mixture graphical model for inferring latent diffusion networks from heterogeneous cascade data, with theoretical guarantees and empirical validation on synthetic and real-world datasets.

The work makes a contribution by introducing a novel modeling framework that captures heterogeneous diffusion patterns through overlapping network structures. The approach is well-motivated, mathematically rigorous, and supported by clear identifiability and optimization guarantees. The empirical studies, including both synthetic cascades and sociological topic diffusion data, demonstrate ability to recover meaningful latent structures and outperform baselines. 

For the camera-ready version, the authors should expand the scalability analysis with clearer discussion of limitations and potential improvements, including practical constraints when networks grow larger. They should emphasize robustness to violations of the low-rank assumption and discuss how imprecise infection times may affect applicability. Adding more qualitative validation in the real-world case study and improving clarity of notation and figure captions would further strengthen the final version.
Decision: Accept (poster)