# Distinguishing simple and complex contagion processes on networks

Giulia Cencetti  Fondazione Bruno Kessler, Trento, Italy  Diego Andrés
Contreras  Aix-Marseille Univ, Université de Toulon, CNRS, Centre de Physique
Théorique, Turing Center for Living Systems, Marseille, France  Marco
Mancastroppa  Aix-Marseille Univ, Université de Toulon, CNRS, Centre de
Physique Théorique, Turing Center for Living Systems, Marseille, France  Alain
Barrat  Aix-Marseille Univ, Université de Toulon, CNRS, Centre de Physique
Théorique, Turing Center for Living Systems, Marseille, France

###### Abstract

Contagion processes, including disease spreading, information diffusion, or
social behaviors propagation, are often schematized as processes evolving on
networks of interactions, either as simple contagion, i.e. involving one
connection at a time, or as complex contagion, in which multiple interactions
are needed for a contagion event. Empirical data on spreading processes
however, even when available, do not easily allow to uncover which of these
underlying contagion mechanisms is at work. We propose here a strategy to
discriminate between several contagion mechanisms upon the observation of a
single instance of a spreading process. The strategy is based on the
observation of the order in which network nodes are infected, and on the
measure of its correlation with their local topology. The relation between the
order of infection and the topology depends indeed on the nodes susceptibility
to infection and hence on the contagion mechanism. We obtain a classifier
that, for each single realization of a spreading process, predicts with high
accuracy the category of the underlying contagion mechanism. In particular,
the classifier discriminates between processes of simple contagion, processes
involving threshold mechanisms and processes driven by group interactions
(i.e., by “higher-order” mechanisms). Its performance remains robust against
partial observation of the process and for processes unfolding on previously
unknown networks. Our results represent an important advance in the
understanding of contagion processes by providing an efficient method using
only limited information to distinguish between several possible contagion
mechanisms.

Phenomena that can be described as contagions, such as disease spreading,
information diffusion, or propagation of social behaviors, are widely present
in nature and in human societies [1, 2, 3, 4, 5, 6]. Modelling contagion
processes in a population typically includes two main steps. First, one
describes how the state of the hosts (individuals who receive and propagate
the disease/information/behavior) can evolve. For instance, one often assumes
that they can only be in one of few possible states, such as susceptible
(healthy), infectious (having the disease/information and able to transmit
it), or recovered (cured and immunized). Second, the propagation is described
along the structure of interactions between hosts, often encoded through a
network in which nodes represent hosts and links represent their interactions.
The resulting framework of network epidemiology has been applied to the spread
of human and computer viruses [7, 4, 6], of rumors [8, 9], of innovations [10,
11, 12, 13, 14] or of behavior [15].

Depending on the phenomenon to model however, the fundamental propagation
mechanisms are different. To describe the spread of infectious diseases,
models of simple contagion, in which a single interaction between a
susceptible and an infectious can lead to a transmission event, are adequate
[1, 6]. In other cases, such as social contagion of behaviors, peer influence
and reinforcement mechanisms can play an important role, and growing empirical
evidence indicate that single interactions are not sufficient to cause
transmission [15, 16, 17, 18, 19]. These cases are hence better described by
so-called complex contagion models, in which each transmission event requires
interactions with multiple infectious hosts. For instance, threshold models
assume that a susceptible node can be infected only if a given fraction of its
neighbours are infectious, in order to mimic the fact that an individual may
adopt an innovation only if enough friends are already adopters [12, 18].
Another type of models takes into account the fact that interactions can occur
not only in binary fashion (as implicit when using networks as a
representation) but also involve group (“higher-order”) interactions [20]:
this can lead to reinforcement mechanisms, in which the tendency of a host to
become infected depends non-linearly on the number of infectious hosts in the
group [21, 22, 23].

For both simple and complex contagions, most studies start from a propagation
mechanism and design models to represent it and study how the structure of the
interaction network impacts the spread [6]. Here, we consider the reverse
question: given an observed single instance of a contagion process on a
network, can we determine whether it is governed by a simple or complex
contagion, whether transmission events are influenced by the state of other
nodes, and whether threshold or group effects are at work? Such issues are
particularly challenging, as highlighted by the fact that patterns and
phenomenology of complex contagion processes at the population level can be
indistinguishable from those of two simple contagion processes in interaction
[24, 25]. Previous works in this direction have shown the limitations of
simple contagion processes by investigating for instance how the probability
of infection of a host depends on the number and repartition of infectious
hosts in their neighborhood [26, 17, 19, 18], or by showing that the temporal
evolution of the number of infectious can in specific cases be better
represented by a complex contagion model [18, 27, 28]. Deep learning
approaches have also been used [29], being however rather data demanding.
Other works have highlighted how different network structures favor different
types of contagion: complex contagion is typically favored by local
cohesiveness, while simple contagion is favored by small network diameters
[15]. This approach however requires to observe the process under
investigation on different networks with different characteristics [30]. In
general thus, we still lack a clear identification of the main features
distinguishing simple and complex contagion processes.

Figure 1:  The first row reports a toy network with 4 infected nodes (black)
and how the different models of propagation would imply contagion of further
nodes (colored). The second row gives scatterplots of $k_{i}$ vs. $o_{i}$ and
$k_{i}^{\Delta}/k_{i}$ vs. $o_{i}$ for single numerical realizations of each
model on the workplace data set, where $o_{i}$ is the order in which the node
$i$ with degrees $k_{i}$ and $k_{i}^{\Delta}$ has been reached by the
propagation in that realization. The values of the corresponding correlation
coefficients are given in the plots. Parameters: $\beta=0.005$ in panels (e)
and (f), $\theta=0.007$ in (g) and (h), $\beta=0.005$, $\beta^{\Delta}=0.8$ in
(i) and (j), $\nu=4$, $\lambda=0.001$ in (k) and (l). The colored curves
report the mean of all $k_{i}$ or $k_{i}^{\Delta}/k_{i}$ for each $o_{i}$.

Here, we put forward a new method based on the correlations between the order
in which successive nodes in the network are reached by the spread and their
basic local properties. We show that different paradigmatic models of simple
and complex contagion lead to different correlation patterns, and how to
exploit them to build a classification tool able to determine whether a given
instance of a spread is due to a simple contagion, a threshold model or a
higher-order contagion model. We emphasize that our method applies to single
instances of a spread. Moreover, it does not require any information on the
state of a host’s neighbours when it becomes infectious. We investigate the
robustness of our classifier with respect to incomplete data, as well as the
possibility to apply it to a process taking place on a different network than
the one on which the classifier has been calibrated.

#### Models of contagion.

We consider four models of contagion process on networks. For simplicity, we
use SI models, i.e., each node can only be in two states, susceptible (S) or
infected (I), and infected nodes do not recover. We consider processes in
discrete time, and the models differ in the mechanism determining how a node
can switch from the S to the I state.

We first consider a simple contagion process (Fig. 1(a)): every susceptible
node can be infected independently by each of its infected neighbors with a
probability per unit time $\beta$. The disease spreads thus following the
network structure of pairwise links among nodes.

The second process (Fig. 1(b)) is the threshold model [12]. In this
deterministic model, a susceptible node becomes infected when the fraction of
its neighbors that are infected reaches a threshold $\theta$. The example
sketched in Fig. 1(b) shows the case where $\theta=0.5$, so only the nodes
which have at least half of their neighbors infected (black nodes) will be
infected. This can represent for instance a situation where being infected
means having a certain opinion, and a node can be convinced by its peers,
eventually conforming itself to its neighborhood.

We also consider two models of complex contagion that involve higher-order
contagion mechanisms, i.e., interactions among groups of nodes [20]. These
models are defined on hypergraphs or simplicial complexes, in which a
hyperedge of size $m$ represents an interaction among a group of $m$ nodes
(simplicial complexes are hypergraphs $\mathcal{H}$ such that, for each
hyperedge –simplex– $e=\\{i_{1},...,i_{m}\\}$, all subsets of $e$ are also
hyperedges of $\mathcal{H}$). In particular, we use the simplicial contagion
[21] model, the simplest generalization of simple contagion to group
interactions. We only consider simplices up to the second order, i.e. three
nodes interactions, neglecting structures of higher orders (which thus appear
only as decomposed into second order simplices). Each susceptible node can be
infected by an infectious neighbour with probability $\beta$ per time step (as
in the simple contagion), but also if it belongs to a simplex of $3$ nodes in
which the two others are infectious. This second process happens with
probability $\beta^{\Delta}$ per time step (see Fig. 1(c)). This model is thus
a mix of simple and complex contagion processes. In addition, we consider the
non-linear hypergraph (NL-hyper) model [23], which includes contagion
processes in interactions of arbitrary sizes: if in a hyperedge of size $m$
there are $n$ infected individuals, each of the remaining $m-n$ susceptible
nodes is independently infected with probability $\lambda n^{\nu}$ at each
time step (Fig.1(d)).

Given an observed single realization of one of these models on a network, our
goal is to devise a method to determine which model it corresponds to. To this
aim, we consider several empirical data sets as the substrates on which the
processes unfold. We use data representing physical or online interactions
between individuals in several contexts: a workplace [31], a primary school in
France [32], a high school in France [33] a scientific conference [31], a
hospital [34], a University campus [35], an elementary and a middle school in
Utah [36], and an email dataset [37]. We refer to the Supplementary Material
(SM) for details on these data. In each case, we consider the aggregated
network or hypergraph: nodes represent individuals and a link between two
nodes represents the fact that they have been in contact at least once during
the data collection (similarly, a hyperedge of size $m$ corresponds to the
fact that a simultaneous clique of size $m$ has been observed at least once,
see SM and [21]). Links and hyperedges can also be weighted by their numbers
of appearance in the data. We consider here for simplicity unweighted networks
and hypergraphs and discuss the case of weighted structures in the SM. In the
main text, we moreover give mostly results obtained with the workplace data
set and refer to the SM for the other data sets.

#### Results.

Simple contagion processes on networks are characterized by a hierarchical
dynamics: the large degree nodes are reached early, and a progressive cascade
follows towards small degree nodes [38]. The order in which nodes are infected
is thus influenced by their degrees. This is illustrated in Fig. 1(e) for an
instance of a simulation of a simple contagion on the workplace data set
network. Figure 1(i) and (k) show that similar cascades are observed for the
simplicial and NL-hyper contagions models. The situation is different for the
threshold model. As the infection of a node does not depend on its number of
infected neighbors but on the corresponding fraction, nodes with large and
small degree can have the same infection risk: no correlation is observed
between the order of infection and the degree (Fig. 1(g)).

Figure 2: Results for the workplace data set. (a): Distributions of the
correlation $C_{1}$ between the order of infection and the degree of nodes,
for the four contagion models. (b): ROC curve when using $C_{1}$ to classify
threshold model processes against the other three. (c): Distributions of the
correlation $C_{2}$ between the order of infection and the ratio
$k^{\Delta}/k$ of infected nodes, for the four contagion models. (d): ROC
curve when using $C_{2}$ to classify the simplicial model against simple and
threshold. (e): ROC curve when using $C_{2}$ to classify the NL-hyper model
against simple and threshold. For the stochastic models (simple, simplicial
and NL-hyper) 1000 realizations are implemented for each parameter setting,
always starting with one random infected node. For the threshold model, due to
its deterministic nature, we use only one realization for each different
initial condition, i.e. one for each network node, for each parameter value.
Parameter values: $\beta\in\\{0.005,0.008,0.014,0.023,0.039\\}$ for both
simple and simplicial models, $\beta^{\Delta}=0.8$,
$\lambda\in\\{0.0001,0.001,0.006,0.011,0.015\\}$, $\nu=4$,
$\theta\in\\{0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10\\}$.

We explore this point further in Fig. 2(a), where we show the distribution of
the Spearman correlation $C_{1}=\text{corr}(o,k)$ between the order in which
nodes are infected and their degree $k$ [39], computed for each numerical
realization of each model and for a wide range of parameter values. The
distributions are similar for the simple and higher-order contagion processes,
with similar ranges of (only negative) correlation values, while the
distribution of $C_{1}$ obtained for the threshold process has a broader
support including also positive values. We thus consider the possibility to
use the value of $C_{1}$ to identify whether a given realization results from
the threshold model or from another model. We assess the quality of this
classification through the parametric Receiver Operating Characteristic (ROC)
curve, built as follows: given $c\in[-1,1]$, we classify each instance having
$C_{1}\geq c$ as resulting from a threshold model. If the instance was really
produced by the threshold model, it is a true positive (TP), and else a false
positive (FP). If the instance instead has a correlation $C_{1}<c$, it is
classified as resulting from one of the simple, simplicial or NL-hyper
processes: if this is correct, it is a true negative (TN), while if it was a
threshold process it is a false negative (FN). The ROC curve presents, as $c$
varies, the true positive ratio $TPR=TP/(TP+FN)$, i.e. the fraction of
instances of threshold model that are correctly classified, versus the false
positive ratio $FPR=FP/(FP+TN)$, i.e. the fraction of instances of the other
models that are wrongly classified. The quality of the classifier is
summarized by the area under the ROC curve (AUC), which is $0.5$ for a random
classification and $1$ for a perfect one. Figure 2(b) shows the ROC curve for
the simulations on the workplace data set, with a very high AUC of $0.95$ for
the classification of processes using $C_{1}$.

To identify processes involving higher-order mechanisms, we need to take into
account not only the degree of nodes but their participation to higher-order
structures, and, specifically, whether they are mostly involved in structures
conveying simple contagion (links) or favoring complex higher-order contagion
processes (hyperedges of larger sizes). We thus consider the correlation
between the order in which nodes are infected and the ratio $k^{\Delta}/k$
between the number of hyperedges of size $3$ to which they belong and their
degree: $C_{2}=\text{corr}(o,k^{\Delta}/k)$ 111It is in principle possible to
consider the number of hyperedges of larger sizes to which each node
participates; for the sake of simplicity, we limit here to size $3$, and, if a
node belongs to a hyperedge of size $m>3$, we decompose the hyperedge into
triangles to compute $k^{\Delta}$.. Figure 1 (f), (h), (j) and (l) illustrate
the correlation $C_{2}$ on specific instances of each process, and Fig. 2(c)
shows its distribution over instances of each model. The distributions are
similar for the simplicial and NL-hyper contagion cases on the one hand, and
for the simple and threshold models on the other hand. Figure 2 shows the ROC
curves obtained when using $C_{2}$ to classify instances of the simplicial
model (panel (d)) or of the NL-hyper model (panel (e)) against instances of
simple and threshold processes. A very good performance is obtained in both
cases.

The results of such classification tasks may however depend on the model
parameters. Indeed, while the distributions of correlation values are rather
stable when changing the parameter of the threshold and simple model (see SM
section S6), this is not the case for the higher-order contagion models, which
mix simple and complex processes. In the simplicial model, in particular, for
small $\beta^{\Delta}$ and large $\beta$ transmission events along links
dominate and higher-order effects become less relevant. The quality of the
classifications based on $C_{2}$ decreases then (see SM, Figs. S12 and S15).
We note however that this is in fact not a weakness of the classification,
since in this limit the mechanisms at work in simplicial and simple model do
become more similar to each other (e.g., in the limit $\beta^{\Delta}\to 0$,
the simplicial model reduces to the simple contagion). A similar (but smaller)
effect is observed in NL-hyper model for large values of $\lambda$ and small
values of $\nu$ (SM, Figs. S12 and S16)).

We now combine the previous results using $C_{1}$ and $C_{2}$ to build a
global classifier. We consider in addition the correlations between the order
$o$ of infection of nodes with $k^{\Delta}$ and with $k^{|}$, their number of
purely pairwise links (excluding connections part of higher-order
interactions): respectively $C_{3}=\text{corr}(o,k^{\Delta})$ and
$C_{4}=\text{corr}(o,k^{|})$. We use a Random Forest (RF) classifier [41], a
standard Machine Learning method, to perform the overall classification of
instances of the four models. The performance of this classification task can
be visually assessed by the confusion matrix depicted in Fig. 3(a): it gives
in row $x$ and column $y$ the number of instances of a model $x$ that are
classified as resulting from model $y$. The matrix shows that simple and
threshold model instances are identified almost perfectly, while simplicial
and NL-hyper model instances are confused more easily. When merging the
higher-order model realizations as one unique class of Complex Higer-Order
(C-HO) processes, the resulting confusion matrix (Fig. 3(b)) shows a very high
performance. The accuracy of the classification task with three classes
(simple, threshold and complex higher-order) is indeed very high, as shown in
Fig. 3(c) at varying $\beta$.

We show in the SM that the classifier yields similar results for instances of
the various models simulated on each of the nine interaction networks we
consider. Moreover, while we have up to now considered only SI processes on
unweighted networks, we also show in the SM that the classifier can be
extended to the case of weighted networks by simply replacing the degrees of
nodes with their strengths (see SM). In addition, the SI model does not take
into account that nodes do not remain infectious at all times, and we thus
perform in the SM a preliminary study using simulations of the SIR framework,
in which hosts recover after having been infected, and cannot any more
participate to transmission events. We obtain similar results (see SM),
leaving further investigations of the role of the recovery parameters and of
more complex models to future works.

Figure 3: (a): Confusion matrix with 4 classes: S (simple), C-S (complex
simplicial), C-NL (complex NL-hyper), and C-T (complex threshold). Training is
implemented using 3288 instances and testing using 1096 instances (with
approximately the same number of instances for each model). (b): Confusion
matrix merging C-S and C-NL into C-HO (complex higher-order). Each element of
the confusion matrices reports the number of instances of each specific model
(Actual classes) that have been predicted as belonging to a specific class
(Predicted classes). Thus, the diagonal elements count the number of well
predicted cases and the off-diagonal elements the cases that are not well
predicted. 2466 instances were used for training and 822 for testing. (c):
Accuracy of Random Forest classification with classes S, C-T, C-HO at varying
$\beta$ in simple and simplicial. 2247 instances used for training and 750 for
testing for each value of $\beta$. (d): Confusion matrix obtained by combining
results on 9 different networks. The relative accuracy assembling in classes
S, C-T, C-HO is 0.84. 51609 instances used for training and 17203 for testing.
(e): Confusion matrix obtained by training the classifier on the high school
network and testing it on the conference network. The accuracy when
considering three classes S, C-T, C-HO is 0.79. 6116 instances used for
training and 12928 for testing. The parameters are set as in Fig. 2.

We have thus built a classification method that works for a given substrate:
the classifier is trained using data coming from one network, and applied on
processes run on the same network. Moreover, we have assumed that the order of
contagion and the local properties of all nodes are known. We now examine the
robustness of the method when moving away from such idealized conditions.
First, we consider the more general problem in which instances of the
contagion processes occurring on different data sets are mixed. The resulting
confusion matrix is shown in Fig. 3(d): the classifier is still able to
distinguish between simple, threshold-based and higher-order based contagion
processes even if the processes come from different substrates (see SM, Fig.
S13(a) for accuracy values).

Moreover, we consider the case of a process observed on a “new” network, while
the classifier has been trained on a different data set. We show in the SM
(Fig. S9(a)) that the resulting accuracy depends on the data sets chosen for
training and testing, and remains high in many cases. The best case (accuracy
of $0.84$), is obtained with a classifier trained with instances run on the
high school data set, and tested on instances run on the conference one: the
resulting confusion matrix is displayed in Fig 3(e) (the accuracy at varying
parameters is reported in SM, Fig. S13(b)). Similarly, the classifier can be
trained making use of instances run on multiple data sets and then tested on
intances from a further data set. The best case corresponds to using the
workplace data set for testing while all the others have been used for
training (accuracy of $0.81$, see Fig. S9(b) of the SM).

We finally report in the SM results obtained when relaxing the assumption of a
complete observation of the spreading instance. First, we assume that only a
fraction of the nodes (chosen arbitrary at random) can be observed. The values
of the correlations are thus computed using only the order of infection and
degrees of the observed nodes (both for training and testing). In this case,
the accuracy of the classifier remains high, with values above $0.7$ even when
only $20\%$ of the nodes are observed. This can be understood by the
robustness of the computation of the correlations when randomly removing a
fraction of the data points. Second, we consider a situation in which only the
initial part of the spread can be observed: only the order of infection and
the degrees of the first $h$ infected are used to compute the correlations.
The classifier performances are more impacted by this restriction, with low
accuracy for small values of $h$. In order to observe e.g. the occurrence of a
cascade from large to small degrees, one needs indeed to observe more than the
first data points, and having information beyond the initial phase of the
spread brings more accurate results.

#### Discussion.

We have here developed a framework for classifying contagion processing on
networks through the observation of correlations between the order in which
nodes are reached and their local structural properties. We emphasize that the
classification task (i) uses only local information, without the need to
access the whole network structure, (ii) does not use any information on the
infection status of the nodes’ neighbours, (iii) is applied on single
instances of a process, and (iv) can distinguish between a simple contagion
process, a process driven by a threshold mechanism and a process involving
contagion on higher-order interactions. The proposed classifier yields a very
good accuracy on several real-world networks, remaining robust against partial
observation of the process, and can be trained on known networks and applied
on a process occurring on a different network.

Our work has several limitations worth discussing. First, we have assumed to
have access to the precise values of the degrees for all observed nodes, as
well as the precise ordering of infection. To be more realistic, one could
simulate the use of estimated values by inserting noise in the degree values.
As the classifier is based on the measure of correlations, we expect that its
accuracy should not be strongly impacted. It could however be affected by e.g.
extreme errors such as hubs classified mistakenly as low degree nodes or vice-
versa. Second, we have mostly considered SI models, where nodes do not
recover, and all nodes are finally reached if the network is connected. The
preliminary results shown in the SM for the more realistic SIR model are
promising but deserve further investigations of the impact of the recovery
rate values. The case of more complex models including a latency period
(between the contagion event and the transition to the infectious state)
and/or asymptomatic infectious individuals, are also worth considering.

Finally, we have considered a limited series of data sets as substrate for the
contagion processes, and the classifier performance depends on the network
characteristics. For instance, the classification accuracy is lower for the
schools data sets. This could be due to the strong community structure of
these networks, as the nodes belonging to the community of the seed are then
the first to be infected, regardless of the process. Future works will involve
investigating how different network structures affect the classifier’s
performance, and how using a well chosen set of networks with different
structures could enhance its overall performance when applied to any data.

## References

  * Anderson and May [1992] R. M. Anderson and R. M. May,  _Infectious Diseases of Humans: Dynamics and Control_ (Oxford University Press, Oxford, 1992).
  * Keeling and Rohani [2011] M. J. Keeling and P. Rohani,  _Modeling infectious diseases in humans and animals_ (Princeton University Press, 2011).
  * Centola and Macy [2007] D. Centola and M. W. Macy, American Journal of Sociology 113, 702 (2007).
  * Barrat  _et al._ [2008] A. Barrat, M. Barthelemy, and A. Vespignani,  _Dynamical Processes on Complex Networks_ (Cambridge University Press, Leiden, 2008).
  * Castellano  _et al._ [2009] C. Castellano, S. Fortunato, and V. Loreto, Rev. Mod. Phys. 81, 591 (2009).
  * Pastor-Satorras  _et al._ [2015] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Reviews of modern physics 87, 925 (2015).
  * Pastor-Satorras and Vespignani [2001] R. Pastor-Satorras and A. Vespignani, Physical review letters 86, 3200 (2001).
  * Daley and Kendall [1964] D. J. Daley and D. G. Kendall, Nature 204, 1118 (1964).
  * Moreno  _et al._ [2004] Y. Moreno, M. Nekovee, and A. F. Pacheco, Phys. Rev. E 69, 066130 (2004).
  * Valente [1995] T. W. Valente,  _Network models of the diffusion of innovations_ (Hampton Press, 1995).
  * Rogers [2003] E. Rogers,  _Diffusion of Innovations, 5th Edition_ (Free Press, 2003).
  * Watts [2002] D. J. Watts, Proceedings of the National Academy of Sciences 99, 5766 (2002).
  * Watts and Dodds [2007] D. J. Watts and P. S. Dodds, Journal of consumer research 34, 441 (2007).
  * Iacopini  _et al._ [2018] I. Iacopini, S. Milojević, and V. Latora, Physical review letters 120, 048301 (2018).
  * Centola [2010] D. Centola, Science 329, 1194 (2010).
  * Weng  _et al._ [2012] L. Weng, A. Flammini, A. Vespignani, and F. Menczer, Sci. Rep. 2, 335 (2012).
  * Ugander  _et al._ [2012] J. Ugander, L. Backstrom, C. Marlow, and J. Kleinberg, Proceedings of the national academy of sciences 109, 5962 (2012).
  * Karsai  _et al._ [2014] M. Karsai, G. Iniguez, K. Kaski, and J. Kertész, Journal of The Royal Society Interface 11, 20140694 (2014).
  * Mønsted  _et al._ [2017] B. Mønsted, P. Sapieżyński, E. Ferrara, and S. Lehmann, PloS one 12, e0184148 (2017).
  * Battiston  _et al._ [2020] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Iacopini  _et al._ [2019] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nature communications 10, 1 (2019).
  * Ferraz de Arruda  _et al._ [2021] G. Ferraz de Arruda, M. Tizzani, and Y. Moreno, Commun. Phys. 4, 24 (2021).
  * St-Onge  _et al._ [2022] G. St-Onge, I. Iacopini, V. Latora, A. Barrat, G. Petri, A. Allard, and L. Hébert-Dufresne, Commun. Phys. 5, 1 (2022).
  * Hébert-Dufresne  _et al._ [2020] L. Hébert-Dufresne, S. V. Scarpino, and J.-G. Young, Nature Physics 16, 426 (2020).
  * Lucas  _et al._ [2022] M. Lucas, I. Iacopini, T. Robiglio, A. Barrat, and G. Petri, arXiv , arxiv:2206.07645 (2022).
  * Aiello  _et al._ [2012] L. M. Aiello, A. Barrat, C. Cattuto, R. Schifanella, and G. Ruffo, EPJ Data Science 1, 12 (2012).
  * Fink  _et al._ [2016] C. Fink, A. C. Schmidt, V. Barash, J. Kelly, C. Cameron, and M. Macy, in  _Tenth International AAAI Conference on Web and Social Media_ (2016).
  * Sprague and House [2017] D. A. Sprague and T. House, PloS one 12, e0180802 (2017).
  * Murphy  _et al._ [2021] C. Murphy, E. Laurence, and A. Allard, Nature Communications 12, 1 (2021).
  * Horsevad  _et al._ [2022] N. Horsevad, D. Mateo, R. E. Kooij, A. Barrat, and R. Bouffanais, Nature Communications 13, 1 (2022).
  * G’enois and Barrat [2018] M. G’enois and A. Barrat, EPJ Data Science 7, 11 (2018).
  * Stehlé  _et al._ [2011] J. Stehlé, N. Voirin, A. Barrat, C. Cattuto, L. Isella, J. Pinton, M. Quaggiotto, W. Van den Broeck, C. Régis, B. Lina, and P. Vanhems, PLOS ONE 6, e23176 (2011).
  * Mastrandrea  _et al._ [2015] R. Mastrandrea, J. Fournet, and A. Barrat, PloS one 10, e0136497 (2015).
  * Vanhems  _et al._ [2013] P. Vanhems, A. Barrat, C. Cattuto, J.-F. Pinton, N. Khanafer, C. Régis, B.-a. Kim, B. Comte, and N. Voirin, PLoS ONE 8, e73970 (2013).
  * Sapiezynski  _et al._ [2019] P. Sapiezynski, A. Stopczynski, D. D. Lassen, and S. Lehmann, Scientific Data 6, 1 (2019).
  * Toth  _et al._ [2015] D. J. Toth, M. Leecaster, W. B. Pettey, A. V. Gundlapalli, H. Gao, J. J. Rainey, A. Uzicanin, and M. H. Samore, Journal of The Royal Society Interface 12, 20150279 (2015).
  * Klimt and Yang [2004] B. Klimt and Y. Yang, in  _European conference on machine learning_ (Springer, 2004) pp. 217–226.
  * Barthélemy  _et al._ [2004] M. Barthélemy, A. Barrat, R. Pastor-Satorras, and A. Vespignani, Phys. Rev. Lett. 92, 178701 (2004).
  * Spearman [1987] C. Spearman, The American journal of psychology 100, 441 (1987).
  * Note [1] It is in principle possible to consider the number of hyperedges of larger sizes to which each node participates; for the sake of simplicity, we limit here to size $3$, and, if a node belongs to a hyperedge of size $m>3$, we decompose the hyperedge into triangles to compute $k^{\Delta}$.
  * Ho [1995] T. K. Ho, in  _Proceedings of 3rd international conference on document analysis and recognition_ , Vol. 1 (IEEE, 1995) pp. 278–282.
  * SP [2008] “Sociopatterns collaboration,” http://www.sociopatterns.org/ (2008), accessed: 2022-11-25.
  * Gemmetto  _et al._ [2014] V. Gemmetto, A. Barrat, and C. Cattuto, BMC infectious diseases 14, 695 (2014).

## Distinguishing simple and complex contagion processes on networks:
Supplementary Material

##  Appendix S1 Data description and generation of interaction networks and
hypergraphs

The data sets we consider describe interactions between individuals, either in
terms of physical proximity or online, and are temporally resolved. We thus
preprocess these data in order to build the three types of structures needed
for the numerical simulations: the temporally aggregated pairwise network, the
aggregated hypergraph with all orders of interactions, and the aggregated
triangle network.

First we build the aggregated network: we consider each interaction between
two individuals as a pairwise link of a static network with weight given by
the number of times (consecutive or not) that the connection has appeared in
the temporal data. The aggregated network is used to simulate the simple and
threshold contagion models, and the contagion events of the simplicial
contagion model involving pairwise links.

We then generate the hypergraph by considering as hyperedges all the groups of
simultaneous all-to-all interactions in the temporal data. Each hyperedge is
stored with its weight given by the number of temporal repetitions. This is
used to simulate the NL-hyper contagion process, where each hyperedge is
characterized by its own infection probability.

Finally, we consider each hyperedge of order higher than 2 (involving more
than 3 nodes) and we decompose it into multiple groups of three nodes (in
other words we project these interactions on hyperedges of size 3). The set of
the obtained three-nodes interactions, together with the hyperedges of order
2, forms the aggregated triangle network. This is used to simulate the
infections taking place on triangles in the simplicial contagion model (those
occurring with probability $\beta^{\Delta}$). The union of the aggregated
triangle network and the pairwise network forms indeed a simplicial complex.

The unweighted networks are obtained by setting all the weights equal to 1.

data set | Nb. of nodes | Nb. of pairwise | Nb. of 2nd order | Max. hyperedge
---|---|---|---|---
|  | connections | connections (projected) | size
Workplace | 217 | 4274 | 772 | 4
Hospital | 75 | 1139 | 703 | 5
Conference | 403 | 9565 | 2634 | 9
Primary school France | 242 | 8317 | 5138 | 5
High school France | 327 | 5818 | 2370 | 5
University campus (48h) | 495 | 3486 | 2981 | 6
Elementary school Utah (1d) | 339 | 8763 | 6868 | 8
Middle school Utah (1d) | 591 | 13550 | 4453 | 7
Email | 143 | 1800 | 6578 | 18
Table S1: Size details about the used data sets. The third column reports the
number of pairwise connections in the aggregated network, i.e. all the
connections of the networks reduced to pairwise interactions. The fourth
column reports the number of 2nd order connections in the triangle network,
i.e. the simplicial complex obtained by projecting higher orders on simplices
of 3 nodes. The fifth column reports the number of nodes in the largest
hyperedge.

We use the following publicly available data sets of social interactions (see
also Table S1):

  * •

Sociopatterns: data sets describing pairwise face-to-face interactions among
individuals, collected using RFID devices [42]. These data have a temporal
resolution of 20 seconds and describe contacts between individuals collected
during several successive days. We consider five of these data sets: a
workplace [31], a conference [31], a hospital [34], a high school [33], and a
primary school [43, 32].

  * •

University campus [35]: data set describing pairwise co-location interactions
among students at the campus of the Technical University of Denmark for one
month, recorded by the exchange of Bluetooth radio packets between smartphones
every 5 minutes for 1 month. We consider a reduced version of the original
data set, involving only the first 48 hours, in order to obtain an aggregated
network with a moderate density.

  * •

Contacts among Utah’s School-age Population [36]: data set describing pairwise
interactions among students in an elementary and a middle school in Utah. The
interactions were collected with a temporal resolution of 20 seconds for 2
days via Wireless Ranging Enabled Nodes. We consider a limited version of the
original data set, involving only the first day.

  * •

Email Enron [37]: data set describing group interactions among a company
employees via email messages: each node corresponds to an email address and a
hyperedge includes the sender and all receivers of an email.

##  Appendix S2 Models of contagion

We describe in details the four different contagion models used in the main
text. For the sake of generality we define the models on weighted networks,
making use of the weighted adjacency matrix $W$, which reduces to the
adjacency matrix for non-weighted networks. In order to make clear the
distinction between models we give for each of them the node infection
probability.

In simple contagion processes, infections take place independently on each
link connecting infected to susceptible nodes at rate $\beta$ multiplied by
the weight of the link (normalized such that $\max_{(i,j)}W_{ij}=1$). The
probability that a susceptible node $i$ is infected at a specific time $t$ is
hence given by:

$p_{i}(t)=1-\prod_{j\in I(t)}(1-\beta W_{ij}),$ (S1)

where $I(t)$ is the set of infected nodes at time $t$, so the second term is
the product of the independent probabilities that node $i$ is not infected by
each of its infected neighbors.

The threshold model is a deterministic model where a susceptible node is
infected when the fraction of its neighbors that are infected exceeds a
threshold $\theta$. More generally, if the network is weighted, a susceptible
node is infected when the weight of its connections with infected nodes
divided by the total weight of its connections exceeds the threshold. The
general probability can be described by a Heaviside function:

$p(t)=H\biggl{(}\frac{1}{s_{i}}\sum_{j\in
I(t)}W_{ij}-\theta\biggr{)}=\begin{cases}1\text{ if }\frac{1}{s_{i}}\sum_{j\in
I(t)}W_{ij}\geq\theta\\\ 0\text{ if }\frac{1}{s_{i}}\sum_{j\in
I(t)}W_{ij}<\theta\end{cases}\,.$ (S2)

where $s_{i}=\sum_{j}W_{ij}$ is the total strength of node $i$.

The simplicial contagion model can be considered as a generalization of the
simple model on simplicial complexes, where nodes are connected via simplices
instead of only pairwise links. We only consider simplices up to the second
order, i.e. three nodes interactions, neglecting structures of higher orders
(which will only appear decomposed into edges of first and second order). Each
susceptible node can be infected via its first-order connections, at rate
$\beta$ (like in simple contagion), or via its second-order connections, if
present, at rate $\beta^{\Delta}$. Second order contagions can only happen if
the susceptible node is part of one (or more) closed triangles where the other
two nodes are both infected. Notice that in this case, as the model is defined
on a simplicial complex, the two nodes additionally contribute independently
to the first-order contagion. All infection probabilities are multiplied by
the weight of the involved simplex, via the matrix $W$ at first order and
tensor $W^{\Delta}$ at second order, where the element $W^{\Delta}_{ikl}$ is
the weight of the closed triangle $(i,k,l)$, normalized by the same factor as
the link weights. The probability per unit time that a susceptible node is
infected is hence given by the combination of its first and second order
interactions with infected nodes:

$p_{i}(t)=1-\prod_{j\in I(t)}(1-\beta W_{ij})\prod_{k,l\in
I(t)}(1-\beta^{\Delta}W^{\Delta}_{ikl})\ .$ (S3)

The non-linear hypergraph contagion model is defined on hypergraphs rather
than on simplicial complexes, and multi-body interactions are modelled through
a non-linear infection term. In this case, we consider the SI epidemic model
on a hypergraph with hyperedges of arbitrary size (i.e. links, triangles and
also higher-order structures up to the maximum hyperedge size $M$): if in a
hyperedge $h=(i_{1},i_{2},...i_{m})$ of size $m$ there are $n_{h}$ infected
individuals, each of the remaining $(m-n_{h})$ susceptible nodes is
independently infected with probability $\lambda n_{h}^{\nu}W^{m}_{h}$. Note
that $\nu$ sets the non-linearity, $\lambda\in[0,1]$ and $m\in[2,M]$. All the
infection terms are multiplied by the weight of the involved hyperedge, via
the tensor $W^{m}$ at the $(m-1)$th order, whose element
$W^{m}_{h}=W^{m}_{i_{1},i_{2},...,i_{m}}$ is the weight of the hyperedge
$h=(i_{1},i_{2},...i_{m})$. By assuming that each hyperedge produces an
independent contagion process, the probability per unit time that a
susceptible node $i$ is infected is hence given by:

$p_{i}(t)=1-\prod_{m=2}^{M}\prod_{h\in\mathcal{E}^{m}_{i}}(1-\lambda
n_{h}^{\nu}W^{m}_{h})$ (S4)

where $\mathcal{E}^{m}_{i}$ is the set of all hyperedges of size $m$
containing $i$ and $n_{h}$ is the number of infected individuals in $h$.

##  Appendix S3 Classification tasks for processes on all the data sets
considered

###  S3.1 SI models on unweighted networks

The analyses reported in the main text using the workplace network can be
repeated with the networks described in section S1. The results for SI model
on unweighted networks are shown in Figs. S1 and S2, which report the
distributions of the correlations $C_{1}$ and $C_{2}$, the resulting ROC
curves and the confusion matrices, similarly to Fig. 2 and Fig. 3(a,b) in the
main text. It is evident that the results obtained in the main text can be
generalized to other networks, showing the versatility of the proposed
classifier (with a classification accuracy always above 0.7). We however
notice that the classification performances are worse for the networks of the
two schools in France and the elementary school in Utah, with in particular a
lower ability to distinguish between simple and simplicial model when using
$C_{2}$ (even if the AUC is still above 0.7). This could be due to the
existence of a strong community structure in these networks, but future work
is needed to investigate how different network topologies affect the
classification based on the measure of correlations.

Figure S1: Unweighted SI model. Analogous of Fig. 2 and Fig. 3(a) of the main
text for four other data sets: conference, hospital, high school, and primary
school. For each data set, the top panels report the distributions of $C_{1}$
and $C_{2}$. The three panels below give: in orange, the ROC curve when using
$C_{1}$ to classify threshold model processes against the other three; in
light green, the ROC curve when using $C_{2}$ to classify the simplicial model
against simple and threshold; in green, the ROC curve when using $C_{2}$ to
classify the NL-hyper model against simple and threshold. The bottom panels
report the confusion matrices with 4 classes (simple, complex threshold,
complex simplicial and complex NL-hyper) and with 3 classes (simple, complex
threshold and complex higher-order). For the stochastic models (simple,
simplicial and NL-hyper) 1000 realizations are implemented for each parameter
setting, always starting with one random infected node. For the threshold
model, due to its deterministic nature, we use only one realization for each
different initial condition, i.e. one for each network node, for each
parameter value. Confusion matrices are obtained using the same number of
instances for each model. Parameter values:
$\beta\in\\{0.005,0.008,0.014,0.023,0.039\\}$ for both simple and simplicial
models, $\beta^{\Delta}=0.8$,
$\lambda\in\\{0.0001,0.001,0.006,0.011,0.015\\}$, $\nu=4$,
$\theta\in\\{0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10\\}$. Additional
data sets are shown in Fig. S2.

Figure S2: Unweighted SI model. Same as Fig. S1 for the data sets: University
campus, elementary school in Utah, middle school in Utah, and email.

###  S3.2 SI models on weighted networks

The classifier can be applied to spreading processes occurring on weighted
networks: the processes spread according to the links’ and hyperedges’
weights, as described in section S2. The classifier takes as input the
correlations between the order of infection of the nodes and their strengths
instead of their degrees.

The results are reported in Figs. S3, S4 and S5. Overall, the classifier
performs well in these cases too (accuracy always above 0.7). For some data
sets, we observe however a reduced performance with respect to the unweighted
case, especially for the ability to distinguish between higher-order complex
contagion and simple contagion. We have checked that the weights of hyperedges
of size larger than 2 are typically much smaller than the weights of the
links, so that the spreading processes along links dominate even for higher
order processes: as a result, higher-order complex and simple contagion become
more similar processes.

Figure S3: Weighted SI model. Analogous to Figs. 2 and 3 for the workplace
data set, using the weighted aggregated networks for the spreading processes.
The parameters are chosen as in Fig. S1, except for $\theta$, which becomes a
threshold for nodes strength instead of nodes degree and is consequently
adapted, taking the values
$\\{0.04,0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.22\\}$.

Figure S4: Weighted SI model. Analogous of Fig. S3 for other networks:
conference, hospital, high school in France, and primary school in France.
Additional data sets are shown in Fig. S5.

Figure S5: Weighted SI model. Same as Fig. S4 for the data sets: University
campus, elementary school in Utah, middle school in Utah, and email.

###  S3.3 SIR models on unweighted networks

The above analyses can be repeated using an SIR instead of an SI model, i.e.
including the possibility for nodes to recover (with a fixed probability per
unit time $\mu$) after they have been infected. A recovered node cannot infect
other nodes or be infected again. The process ends when none of the nodes is
infected (they are all susceptible or recovered).

The results of the classifier in this case are reported for unweighted
networks in Figs. S6, S7 and S8. The classifier performances are slightly
reduced with respect to the SI case. This can be understood by the fact that
the nodes remain infected only a limited amount of time (independent of their
topological properties and of their neighbors’ properties). When a node
recovers, it does not anymore contribute to the risk of infection of its
neighbours. When reinforcement effects are needed (complex contagion), this
can thus alter substantially the relationship between the structure of a
node’s neighborhood and its risk of becoming infected, especially if the
recovery is fast.

Here, as a preliminary analysis, we have considered simulations with
$\mu=0.1$, and the classifier performances remain good, with accuracy values
always above 0.7. A future deeper analysis will investigate in more details
the impact of the recovery parameter on the classifier performance.

Figure S6: Unweighted SIR model. Analogous of Figs. 2 and 3 for the SIR model
with $\mu=0.1$ with the workplace network. Additional data sets are shown in
Figs. S7 and S8.

Figure S7: Unweighted SIR model. Analogous of Fig. S1 for the SIR model with
$\mu=0.1$. Results for additional data sets are shown in Fig. S8.

Figure S8: Unweighted SIR model. Same as Fig. S7 for additional data sets:
University campus, elementary school in Utah, middle school in Utah, and
email.

##  Appendix S4 Classification across data sets

We report here the results of classification experiments in which we train and
test the classifier on instances of the spreading processes taking place on
distinct networks.

In a first experiment, we train the classifier using instances of spreading
processes run on one data set, and we test on instances simulated on a second
data set. The classification accuracy depends on the networks chosen for
training and testing, as reported in the table in Fig. S9(a), where the
results above 0.7 are highlighted (diagonal elements represent the cases where
training and testing are implemented using the same data set). The same
results are also depicted on the right as a network of data sets: arrows start
from the data set used for training and end at the data set used for testing.
Only arrows associated to an accuracy above 0.6 are drawn, in purple for
accuracy between 0.6 and 0.7, in turquoise above 0.7.

The second experiment instead involves training on eight different data sets
and testing on the ninth one. The table in Fig. S9(b) reports the accuracy
obtained for testing on each data set. The highest accuracy is obtained when
the testing data set is the workplace (we report the corresponding confusion
matrix on the right of the figure).

In all these experiments we consider the same model parameter values as in
Fig. 2:

$\beta\in\\{0.005,0.008,0.014,0.023,0.039\\}$ for both simple and simplicial
models,

$\beta^{\Delta}=0.8$,

$\lambda\in\\{0.0001,0.001,0.006,0.011,0.015\\}$,

$\nu=4$,

$\theta\in\\{0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10\\}$.

Figure S9: (a): Train on one data set, test on another one. The table reports
the accuracy for each pair of data sets, while the network only reports the
pairs with accuracy above 0.6, where arrows go from the training to the
testing data set. (b): Train on eight data sets, test on the ninth. The table
reports the accuracies and the confusion matrix is shown for the best case:
testing on workplace with training implemented on the other data sets.

##  Appendix S5 Robustness with respect to incomplete observations

The results presented in the main text have been obtained by observing the
entire infection chains, i.e. all the infected nodes with their infection
order. Here we explore a more realistic case where the access to infection
data is limited and only a partial observation of infection chains is
possible. In the two following sections we show how the proposed
classification changes with respect to a reduction in the number of observed
infections: first observing just a random percentage of nodes, secondly
observing just the first infected.

###  S5.1 Observing only a fraction of the nodes

Let us assume that only a random subset of nodes can be observed. We rely on
the same numerical simulations that are presented in the main text, with the
difference that the correlation between infection order and topological
features is computed only for an a priori chosen subset of nodes, ignoring all
the others. We repeated the experiment for 10 realizations of the random
subsets, and for each of the following subset sizes: $20\%$, $40\%$, $60\%$,
$80\%$ and $100\%$ of the network nodes. For each experiment we computed the
Random Forest accuracy of the classification in three classes
$\\{$Simple/Threshold/Higher-order$\\}$, the area under the ROC curve when
using $C_{1}$ to single out the threshold model vs. the others and the AUC
when using $C_{2}$ to single out the higher-order models vs. the simple and
threshold ones. The results are reported in Fig. S10 for the nine data sets
considered: the classification remains very efficient in all cases.

Figure S10: Random Forest accuracy and AUC when using $C_{1}$ and $C_{2}$, as
a function of the number of observed nodes in the different data sets. The
$x$-axes report the number of observed nodes, corresponding to $20\%$, $40\%$,
$60\%$, $80\%$ and $100\%$ of the total. The dots and bars show average and
standard deviation for results collected on 10 realizations of the random
subsets for each size.

###  S5.2 Observing only the first infected nodes

As a second experiment, we try to classify contagion models by observing only
the first infected nodes. The number of nodes to observe is chosen as in the
previous experiment, corresponding to $20\%$, $40\%$, $60\%$, $80\%$ and
$100\%$ of the total size. In this case, Figure S11 shows that the performance
depends more strongly on the fraction of observed nodes, and is smaller with
respect to the case where the same number of nodes have been chosen at random.
Combined with the results of Fig. S10, this suggests that a reduction of the
number of observed nodes does not have a drastic effect if they are
distributed along the entire infection chain, but the method is less robust if
only the first phase of the infection pattern is observed.

Figure S11: Random Forest accuracy and AUC when using $C_{1}$ and $C_{2}$,
obtained by observing only the first infected nodes in each infection chain.
The number of observed nodes is reported on the $x$-axes, corresponding to
$20\%$, $40\%$, $60\%$, $80\%$ and $100\%$ of the total number of nodes.

##  Appendix S6 Influence of model parameters

The classification experiments shown in the main text take into account the
four different models of contagion with several choices of each model’s
parameters. In some cases the classifier performances depend on these
parameters, while in others the parameters choice is less relevant, as
discussed in the following.

###  S6.1 Parameters of the higher-order models

The parameters’ choice is particularly relevant when the variation of these
parameters can make two models’ contagion mechanisms more similar. For
instance the simplicial model is an extension of the simple one and reduces to
a simple propagation for $\beta^{\Delta}\to 0$. The dynamics – and hence the
order at which nodes are infected – in the simplicial model appears thus more
similar to that of simple model for smaller values of $\beta^{\Delta}$ and
larger values of $\beta$: in such cases, contagion events take place more via
pairwise links than via triangles, weakening the importance of the higher-
order mechanism. This naturally affects the classification performances, as
shown in Fig. S12(a), which shows the area under the ROC curve when using
$C_{2}$ to distinguish instances of the simplicial model from instances of the
simple propagation, for varying values of $\beta$ and $\beta^{\Delta}$ (the
instances of the simple model correspond to values of $\beta$ ranging between
0.005 and 0.1). The AUC decreases when $\beta$ increases and when
$\beta^{\Delta}$ decreases.

Analogously, in Fig. S12(b) we show the AUC obtained when using $C_{2}$ to
classify instances of the NL-hyper model and instances of the simple model, at
varying $\nu$ and $\lambda$. In this case, the AUC decreases slightly when
$\lambda$ increases and $\nu$ decreases, but remains very high for all
explored parameter values.

Figure S12: AUC of the ROC curve when using $C_{2}$ to distinguish simplicial
vs. simple contagion models (a) and NL-hyper vs. simple contagion (b), for
varying parameters. Each colored pixel is the result of 1000 experiments. The
spreading processes are simulated on the workplace data set.

In Fig. 3(d) of the main text we showed the confusion matrix obtained when
classifying the four contagion models using all the nine considered networks.
We depict here the confusion matrix obtained when considering the two higher-
order models as a single class, and we show the related accuracy obtained at
varying the parameter $\beta$ in simple and simplicial contagion (Fig.
S13(a)). Analogously we show the confusion matrix and the accuracy obtained at
varying $\beta$ relative to the experiment reported in Fig. 3(e), where the
classifier is trained on the High School France data set and tested on the
conference data set (Fig. S13(b)).

Figure S13: Confusion matrix and classification accuracy vs. $\beta$ in simple
and simplicial model in the case of classification with all data sets (a) and
in the case of the classification on the conference data set after training on
the High School France network (b). For (a) five different splits of the
correlation measures in training and testing sets have been considered, the
plot reports mean and standard deviation.

In Figs. S14, S15, and S16 the distributions of $C_{1}$ and $C_{2}$ are shown
at varying parameters $\beta$ and $\lambda$. While the correlation
distributions obtained with the simple contagion model are not particularly
sensitive to the value of its unique parameter $\beta$, the effect of $\beta$
is instead important in the simplicial model, where it plays a role in the
competition between contagion via links and via triangles. Smaller values of
$\beta$ (i.e. more contagions on triangles) lead to smaller values of the
correlation $C_{2}$ (larger in absolute value), and thus to larger differences
with respect to the simple contagion model. A similar effect is observed in
the NL-hyper model at decreasing $\lambda$.

Figure S14: Simple contagion model. Distributions of correlations $C_{1}$ and
$C_{2}$ at varying parameter $\beta$ in unweighted SI simulations on the nine
considered networks. Figure S15: Simplicial contagion model. Distributions of
correlations $C_{1}$ and $C_{2}$ at varying parameter $\beta$ (while
$\beta^{\Delta}$ is fixed at 0.8) in unweighted SI simulations on the nine
considered networks. Figure S16: NL-hyper model. Distributions of correlations
$C_{1}$ and $C_{2}$ at varying parameter $\lambda$ (while $\nu$ is fixed at 4)
in unweighted SI simulations on the nine considered networks.

###  S6.2 Threshold parameter in threshold model

The correlations measured in the threshold model are almost always around
zero, since the nodes susceptibility usually does not depend on their local
neighborhood. This is hence almost independent on the threshold parameter
$\theta$, as shown in Fig. S17. We however notice that for certain networks
smaller values of $\theta$ correspond to smaller values of correlations
$C_{1}$ and $C_{2}$. This may be due to the fact that the smaller is $\theta$
the less important is the threshold effect, and the more similar the threshold
model becomes to the simple contagion one. The effect is however small so
overall $\theta$ does not affect significantly the classifier performance.

Figure S17: Threshold model. Distributions of correlations $C_{1}$ and $C_{2}$
at varying parameter $\theta$ in unweighted SI simulations on the nine
considered networks.


11institutetext: Moscow Institute of Physics and Technology, Russia
22institutetext: Weierstrass Institute for Applied Analysis and Stochastics,
Germany

# Stochastic optimization for dynamic pricing††thanks: The research is
supported by the Ministry of Science and Higher Education of the Russian
Federation (Goszadaniye) №075-00337-20-03, project No. 0714-2020-0005.

Dmitry Pasechnyuk 11 0000-0002-1208-1659 Pavel Dvurechensky 22
0000-0003-1201-2343 Sergey Omelchenko 11 0000-0002-7266-3917 Alexander
Gasnikov 11 0000-0002-7386-039X

###### Abstract

We consider the problem of supply and demand balancing that is stated as a
minimization problem for the total expected revenue function describing the
behavior of both consumers and suppliers. In the considered market model we
assume that consumers follow the discrete choice demand model, while suppliers
are equipped with some quantity adjustment costs. The resulting optimization
problem is smooth and convex making it amenable for application of efficient
optimization algorithms with the aim of automatically setting prices for
online marketplaces. We propose to use stochastic gradient methods to solve
the above problem. We interpret the stochastic oracle as a response to the
behavior of a random market participant, consumer or supplier. This allows us
to interpret the considered algorithms and describe a suitable behavior of
consumers and suppliers that leads to fast convergence to the equilibrium in a
close to the real marketplace environment.

###### Keywords:

Automatic pricing Expected revenue function Stochastic convex optimization
Supply and demand balancing.

##  1 Introduction

With the development of platforms for online trading and services provision,
the problem of dynamic pricing becomes more and more urgent. The environment
of online marketplaces and financial services raises the question of
establishing the most relevant prices for the items presented, and in such
cases, relevance is understood as the possibility of market clearing, i.e.
balancing supply and demand.

There are various approaches to finding equilibrium prices. In this paper, we
consider an approach based on the characterization of the market state by a
potential function of total expected revenue, similar to the function of total
excessive revenue from [13]. In [9], the authors, using this approach, analyze
the case in which consumers follow a discrete choice model with imperfect
behavior introduced by random noise in the assessment of utility, and
suppliers seek to maximize their profit, taking into account the quantity
adjustment costs. This paper substantiates the convexity of the used potential
as a function of prices. Hence, it follows that there are prices that minimize
the potential, and it turns out that this minimum satisfies the condition of
market clearing.

Thus, if there is an intermediary responsible for the formation of prices, it
only needs to build a sequence of price values leading to a minimum. At the
same time, it seeks to use the smallest possible amount of observations of
market participants. Note that the introduced formulation is very convenient
for applying the results of optimization theory to find the optimal prices.
Indeed, the function under consideration is to be minimized and it
simultaneously includes the characteristics of both consumers and suppliers,
although in the classical game-theoretic approach the opposition of the
interests of these parties leads through the Nash equilibrium scheme to
minimax problems. Nevertheless, it turns out to be possible to propose exactly
the potential function that describes the whole system and reduces its dynamic
to only tendency to an extreme point, similar to physical systems. On the
other hand, the considered function has the form of a sum including together
the terms that characterize the behavior of both parties and each agent. This
means that we can consider both consumers and suppliers as the same market
participants and uniformly take into account their interests, without
distinguishing between their types. Together, the described advantages lead to
the ability to use the developments of convex optimization to find methods for
constructing a sequence of prices, and analyze them in terms of the oracle
complexity, which in this setting means the number of observations of the
agents’ behavior.

In this paper, continuing the ideas of [9], the problem of finding equilibrium
prices is posed in a stochastic setting, and the number of observations of the
intermediary for single market participants is determined as a measure of the
effectiveness of the method for solving this problem. Section 2 describes the
problem statement and introduces the stochastic oracle for considered function
to optimize. Section 3 discusses various algorithms for stochastic convex
optimization, proposes theoretical estimates of the efficiency of the
algorithms, and describes the practical advantages of each of them. Section 4
generalizes the stochastic setting for the case of an infinite number of
consumers, which makes it possible to take into account the previously
unobserved consumers, and describes the algorithms for this setting. Section 5
considers a special case of the problem with zero quantity of the adjustment
costs, in which the potential loses its smoothness property, and proposes
several approaches to get around this difficulty. Section 6 shows the results
of modeling of the behavior of the proposed algorithms for a synthetic problem
and demonstrates the improvement in the efficiency of the methods, which can
be achieved by considering the stochastic formulation of the problem, in
comparison with the dynamics from [9].

##  2 Problem statement

Let us imagine the marketplace environment. There are a number of suppliers
(for example, shops) and a number of consumers (buyers). And every supplier
offers some products (goods) in the assortment. Products are grouped by their
type, every group contains some alternatives. In turn, every consumer can
choose to buy one of the alternatives, guided by its subjective utility and
price (with some element of randomness). For reasons of increasing profit,
suppliers can change the assortment, taking into account the costs of these
changes themselves. We use the corresponding mathematical models for the
consumer and the supplier, and characterize them by the values of expected
surplus and maximal revenue, respectively. Summing these terms up we obtain
the function of total expected revenue. It describes the current imbalance in
the market system. Its minimum corresponds to the prices at which the market
is cleared. This means that we can formulate an optimization problem for
finding equilibrium prices.

Let us consider such an optimization problem (in a simplified form) for the
case of $n$ product alternatives divided into $m$ disjoint groups
$G_{i}\subset\\{1,...,n\\}$, $S$ suppliers with convex costs functions
$c_{s}:\mathbb{R}^{n}_{+}\rightarrow\mathbb{R}_{+}$, closed and convex sets of
capacity constraints $Y_{s}\subset\mathbb{R}^{n}_{+}$, typical supplies
$\hat{y}_{s}\in Y_{s}$ and quantity adjustment costs equal to
$-\Gamma_{s}\cdot\|y-\hat{y}_{s}\|_{2}^{2}$ for some parameter $\Gamma_{s}>0$,
$D$ consumers with matrix $A=\\{a_{id}>0\\}_{1,1}^{n,D}$ of alternatives
subjective utility for each of them, and prices vector $p$:

$\min_{p\in\mathbb{R}^{n}_{+}}\left\\{f(p):=\sum_{s=1}^{S}\pi_{s}(p)+\sum_{d=1}^{D}E_{d}(p)\right\\},$
(1)

where

$\pi_{s}(p)=\max_{y\in Y_{s}}\left\\{\langle y,p\rangle-
c_{s}(y)-\Gamma_{s}\cdot\|y-\hat{y}_{s}\|_{2}^{2}\right\\}$ (2)

is the maximal revenue of supplier with taking into account costs
$c_{s}(\cdot)$ and quantity adjustment costs parametrized by $\Gamma_{s}$ for
given prices, and

$E_{d}(p)=\ln\left(\sum_{j=1}^{m}\left(\sum_{i\in
G_{j}}e^{(a_{id}-p_{i})/\mu_{j}}\right)^{\mu_{j}}\right)$ (3)

is the expected surplus $\displaystyle
E_{d}(p)=\mathbb{E}_{\epsilon}\left[\max_{i}\left\\{a_{id}-p_{i}+\epsilon_{i}\right\\}\right]$
for the discrete choice demand model with noise and for corresponding nested
logit distribution [7], with some correlation parameters $0<\mu_{j}\leq 1$. We
can also provide an explicit expression for the gradient of the introduced
objective function:

$\nabla f(p)=\sum_{s=1}^{S}y_{s}(p)-\sum_{d=1}^{D}x_{d}(p),$ (4)

where $y_{s}(p)$ is the optimal solution of the optimization problem (2) (by
Demyanov–Danskin theorem [2]), and

$[x_{d}(p)]_{i}=\frac{\displaystyle e^{(a_{id}-p_{i})/\mu_{j}}\left(\sum_{k\in
G_{j}}e^{(a_{kd}-p_{k})/\mu_{j}}\right)^{\mu_{j}-1}}{\displaystyle\sum_{h=1}^{m}\left(\sum_{k\in
G_{h}}e^{(a_{kd}-p_{k})/\mu_{h}}\right)^{\mu_{h}}},\quad i\in\\{1,...,n\\}$

are the probabilities to choose certain alternative by the consumer111Since
$\displaystyle
E_{d}(p)=\mathbb{E}_{\epsilon}\left[\max_{i}\left\\{a_{id}-p_{i}+\epsilon_{i}\right\\}\right]=\sum_{i}\mathbb{P}\left[a_{id}-p_{i}+\epsilon_{i}=\max_{i}\left\\{a_{id}-p_{i}+\epsilon_{i}\right\\}\right]\cdot\mathbb{E}_{\epsilon}\left[a_{id}-p_{i}+\epsilon_{i}\right]$,
and therefore $\displaystyle\frac{\partial E_{d}(p)}{\partial
p_{i}}=-\mathbb{P}\left[a_{id}-p_{i}+\epsilon_{i}=\max_{i}\left\\{a_{id}-p_{i}+\epsilon_{i}\right\\}\right]$.

To use optimization methods for dynamic pricing we can utilize that $f$ is of
the form of sum. If the number of suppliers $S$ and consumers $D$ is big,
evaluation of all the term in both sums is too expensive in term of working
time and algorithmic complexity. But if we randomly pick up only one of
$(S+D)$ terms at iteration, computing resources are used much more sparingly,
which leads to faster operation of the algorithm. So, considering only one
term of sum from (4), we introduce the stochastic gradient oracle:

$\nabla f_{i}(p)=\begin{cases}y_{s}(p)&i=s\leq S\\\
-x_{d}(p)&i=S+d\end{cases},\quad i\in\\{1,...,S+D\\}$ (5)

Note, that the use of stochastic oracle is very natural. Indeed, using the
gradient dynamics of prices from [9], we evaluate all of the terms in (5) at
iteration, i.e. it is necessary to consider the behaviour of all the consumers
and suppliers to make one step of dynamic. But in practice we cannot guarantee
that we collect all this information in short time (consumers may be
impermanent, so the waiting time may be arbitrarily long). At the same time,
the decisions of consumers and suppliers are not rigidly connected: we can
observe some number of consumer’s sales daily (represented by $x_{d}(\cdot)$),
and much less often and independently the periodic store assortment changes
(represented by $y_{d}(\cdot)$). Therefore, using the dynamics with stochastic
oracle in the form of (5) we can immediately take into account newly observed
behaviour and make iteration, that is now cheap both in computation and in the
required downtime. However, in real-life environment we alsp cannot estimate
the probability of customer’s choice, represented by the vector $x_{d}(p)$
included in second case of (5). We can obtain only some single sales, those
are in fact the random samples of the form $\mathcal{X}_{d}(p)=(0\;\dotsi
0\;1\;0\;\dotsi 0)$, where $[\mathcal{X}_{d}(p)]_{i}=1$ w.p. $[x_{d}(p)]_{i}$,
and $[\mathcal{X}_{d}(p)]_{j}=0$ for $i\neq j$. Therefore,
$\mathbb{E}[\mathcal{X}_{d}(p)]=x_{d}(p)$, and we can introduce another, more
practical, stochastic gradient oracle for constructing our dynamics:

$\widetilde{\nabla}f_{i}(p)=\begin{cases}y_{s}(p)&i=s\leq S\\\
-\mathcal{X}_{d}(p)&i=S+d\end{cases},\quad i\in\\{1,...,S+D\\}.$ (6)

Now, let us clarify properties of the objective function and introduced
oracles. In our simplified setting the following result holds:

###### Lemma 1

(Theorem 3.7 [9]) $f$ has $L$-Lipschitz continuous gradient w.r.t.
$\|\cdot\|_{2}$ with

$L=\sum_{s=1}^{S}\frac{1}{\Gamma_{s}}+\sum_{d=1}^{D}\frac{1}{\min_{j}\mu_{j}},$

and each $\nabla f_{i}(p)$ for all $i\in\\{1,...,S+D\\}$ is $L_{i}$-Lipschitz
continuous w.r.t. $\|\cdot\|_{2}$ with

$L_{i}=\begin{cases}\displaystyle\frac{1}{\Gamma_{i}}&i\leq S\\\
\displaystyle\frac{1}{\min_{j}\mu_{j}}&i>S\end{cases},\quad
i\in\\{1,...,S+D\\}$

##  3 Algorithms and theoretical guarantees

###  3.1 Stochastic gradient descent

Let us consider the simplest dynamic of prices, based on the classical
stochastic gradient descent method. A very natural interpretation for this
dynamics is that in every iteration we can observe only one of the market
participants, supplier or consumer. As soon as participant makes an economical
decision (consumer chooses the product or supplier modifies the supply plan),
we evaluate the stochastic oracle $\widetilde{\nabla}f_{i}(p_{t})$, and make a
step to the equilibrium prices. This dynamic is listed as Algorithm 1. We use
the $[\cdot]_{+}$ notation for the positive part function, i.e. for
$a=[b]_{+}$ we have $a_{i}=\max\\{0,b_{i}\\}$, and denote by
$i\sim\mathcal{U}\\{1,...,S+D\\}$ the i.i.d. random variables from discrete
uniform distribution.

We analyse this dynamic as the projected stochastic gradient method with
Polyak–Ruppert averaging with tunable parameter $C>0$ to control the step size
of the method. One practical advantage of this method is the robustness to the
choice of $C$: it may be chosen regardless of theoretical value of $L$.
However, the analysis additionally requires the condition of stochastic
gradient’s boundedness, but due to the smoothness of $f$ we can bound the norm
of (5), while the additional randomization in the second case of (6) acts on
the standard simplex and therefore is also bounded. Unfortunately, I didn’t
understand the previous sentence.

Algorithm 1 SGD dynamic

1:$p_{0}$ — starting prices values, $N$ — number of iterations, $C$ —
parameter to control the step size

2:for $t=0,\,1,\,\ldots,\,N-1$ do

3: $i\sim\mathcal{U}\\{1,...,S+D\\}$

4: $\displaystyle
p_{t+1}=\left[p_{t}-\frac{C}{\sqrt{t+1}}\widetilde{\nabla}f_{i}(p_{t})\right]_{+}$

5:end for

6:$\displaystyle\widetilde{p}_{N}=\frac{1}{N}\sum_{t=1}^{N}p_{t}$

7:return $\widetilde{p}_{N}$

###### Theorem 3.1

(Theorem 7 [8]) Let us assume that used stochastic oracle is uniformly
bounded: $\|\widetilde{\nabla}f_{i}(p_{t})\|_{2}\leq B$ for all
$i\in\\{1,...,S+D\\}$ and $t\in\\{1,...,N\\}$. The suboptimality of prices
$\widetilde{p}_{N}$ given by SGD dynamic (Algorithm 1) is decreasing as
follows

$\mathbb{E}[f(\widetilde{p}_{N})]-f(p_{*})\leq\frac{\|p_{0}-p_{*}\|_{2}^{2}+CB^{2}(1+C\ln{N})}{2C\sqrt{N}}.$

Moreover, to obtain the prices satisfying the suboptimality bound

$\mathbb{E}[f(\widetilde{p}_{N})]-f(p_{*})\leq\varepsilon,$

it is sufficient to call $\widetilde{\nabla}f_{i}$ oracle
$\displaystyle\mathcal{O}\left(\frac{1}{\varepsilon^{2}}\right)$ times.

Therefore, considered SGD dynamic with Polyak–Ruppert averaging obtains
convergence rate of $\displaystyle\mathcal{O}\left(N^{-1/2}\right)$, up to a
logarithmic term. This also matches the result from [13].

###  3.2 Adaptive stochastic gradient method

In this section we describe a slightly different AdaGrad [3] dynamics that has
a different step size policy. More precisely, the stepsize is chosen based on
the stochastic subgradients on the trajectory of the method. This allow the
algorithm to adapt to the local information and possibly make longer steps.

Algorithm 2 AdaGrad dynamic

1:$p_{0}$ — starting prices values, $N$ — number of iterations, $\eta$ — step
size parameter, $\epsilon$ — small term to prevent zero division

2:$H_{0}=0$

3:for $t=0,\,1,\,\ldots,\,N-1$ do

4: $i\sim\mathcal{U}\\{1,...,S+D\\}$

5: $g_{t}=\widetilde{\nabla}f_{i}(p_{t})$

6: $H_{t+1}=H_{t}+\langle g_{t},g_{t}\rangle$

7: $\displaystyle
p_{t+1}=\left[p_{t}-\frac{\eta}{\sqrt{H_{t+1}+\epsilon}}g_{t}\right]_{+}$

8:end for

9:$\displaystyle\widetilde{p}_{N}=\frac{1}{N}\sum_{t=1}^{N}p_{t}$

10:return $\widetilde{p}_{N}$

###### Theorem 3.2

(Corollary 4.3.8 [4]) Let $\|p_{t}-p_{0}\|_{\infty}\leq R$ for all
$t\in\\{1,...,N\\}$, $\eta$ is proportional to $R$. The suboptimality of the
prices $\widetilde{p}_{N}$ given by AdaGrad dynamic (Algorithm 2) is
decreasing as follows

$\mathbb{E}[f(\widetilde{p}_{N})]-f(p_{*})\leq\frac{3R}{2N}\cdot\sum_{i=1}^{n}\mathbb{E}\left[\left(\sum_{t=1}^{N}[g_{t}]_{i}^{2}\right)^{1/2}\right].$

To obtain the prices satisfying the suboptimality bound

$\mathbb{E}[f(\widetilde{p}_{N})]-f(p_{*})\leq\varepsilon,$

it is sufficient to call $\widetilde{\nabla}f_{i}$ oracle
$\displaystyle\mathcal{O}\left(\frac{1}{\varepsilon^{\color[rgb]{0,0,0}2}}\right)$
times.

Hence the proposed AdaGrad dynamic demonstrates convergence rate of the order
$\displaystyle\mathcal{O}\left(N^{-1/2}\right)$. This asymptotic is similar to
that for SGD dynamic, but in practice such a simple modification of the step
size allows to discernibly improve the convergence rate. At the same time, the
step size hyperparameter $\eta$ is still free and it allows to manually tune
the algorithm for the best practical efficiency.

##  4 The case of infinite number of consumers

In general, the total expected revenue framework considered in [9] allows one
to describe not only the individual consumers, but also the groups of
consumers with close behavior. At the same time, their model covers only the
setting of finite number of agents, which may be not completely practical.
Indeed, if the number of agents is huge or new agents can enter the
marketplace as time goes, it makes sense to consider the limit when the number
of agents tends to infinity. This leads to a general, non-finite-sum objective
given as an expectation.

To be more specific, we consider that consumers are represented by a random
variable with some unknown distribution, i.e. $d\sim\mathcal{D}$. Then the
characteristic vector $(a_{id})_{i=1}^{n}$ is also random. To maintain the
property of market clearing at optimal prices we also assume that the number
of suppliers is infinite and that they are represented by another random
variable, i.e. $s\sim\mathcal{S}$. The next step is to take the limit in (1)
and make transition to expectation w.r.t. distributions
$\mathcal{D},\mathcal{S}$ instead of sums. For the sake of normalization we
introduce the parameter $0<\beta<1$ that is equal to the fraction of suppliers
among all market participants. Informally,
$\beta=\lim_{(S+D)\rightarrow\infty}S/(S+D)$.

In this way, considering function $\frac{1}{S+D}f(p)$ instead of $f(p)$ and
taking the limit as $(S+D)\rightarrow\infty$, we have the new optimization
problem in the form of expectation:

$\min_{p\in\mathbb{R}^{n}_{+}}\left\\{\widetilde{f}(p):=\beta\cdot\mathbb{E}_{s\sim\mathcal{S}}[\pi_{s}(p)]+(1-\beta)\cdot\mathbb{E}_{d\sim\mathcal{D}}[E_{d}(p)]\right\\}.$
(7)

Note, that to preserve the convergence properties of the SGD dynamic
considered in Section 3.1 it is sufficient just to generalize the used
stochastic oracle (6) to the proposed setting by defining

$\widetilde{\nabla}\widetilde{f}(p):=\begin{cases}y_{s}(p)&\text{for
}s\sim\mathcal{S}\text{ w.p. }\beta\\\ -\mathcal{X}_{d}(p)&\text{for
}d\sim\mathcal{D}\text{ w.p. }1-\beta\end{cases}.$ (8)

Thus, to generate the stochastic gradient, we first with probability $\beta$
choose to choose among suppliers or with probability $1-\beta$ we choose to
choose among consumers. Then, in the former case we sample supplier from the
distribution $\mathcal{S}$ and in the latter case we sample consumer from the
distribution $\mathcal{D}$. Finally, for the chosen agent the stochastic
gradient is defined as in (6).

Since problem (7) is a general stochastic optimization problem, we can apply
Algorithm 1 and obtain the dynamic listed as Algorithm 3. Its interpretation
is quite similar to that of Algorithm 1: at the every iteration we observe the
behaviour of one market participant, consumer or supplier, and change the
prices in a proper way. As in the SGD dynamic, samplings $d\sim\mathcal{D}$,
$s\sim\mathcal{S}$ and switching w.p. $\beta$ are provided by natural flow of
participants, we assume that information about participants decisions arrives
uniformly. Since Lemma 1 and the conditions of Theorem 3.1 still hold, the
convergence rate of dynamic below is similar to the one given in Theorem 3.1.

Algorithm 3 SGD dynamic (online setting)

1:$\beta$ — fraction of suppliers, $p_{0}$ — starting prices values, $N$ —
number of iterations, $C$ — parameter to control learning rate

2:for $t=0,\,1,\,\ldots,\,N-1$ do

3: $\displaystyle
p_{t+1}=\left[p_{t}-\frac{C}{\sqrt{t+1}}\widetilde{\nabla}\widetilde{f}(p_{t})\right]_{+}$

4:end for

5:$\displaystyle\widetilde{p}_{N}=\frac{1}{N}\sum_{t=1}^{N}p_{t}$

6:return $\widetilde{p}_{N}$

##  5 The case of zero quantity adjustment costs

In this section we return to the setting of Section 2 and consider the
limiting case when in (2) $\Gamma_{s}=0$. In this case there is no typical
supply $\hat{y}_{s}$ and all the necessary information about the costs
incurred by the supplier is set by the function $c_{s}(\cdot)$. In this case
there is no guarantee that the function $f$ in (1) is smooth (cf. Lemma 1) and
the algorithms have to be properly modified to guarantee the convergence. This
section is devoted to such modifications.

###  5.1 Mirror descent dynamic

Stochastic mirror descent (SMD) [5, 10, 11] is widely used and theoretically
optimal algorithm for stochastic convex non-smooth optimization problems. We
consider here a particular case of the SMD dynamic with Euclidean proximal
setup, in which it looks quite similar to SGD dynamic, but uses a different
step size policy. SMD dynamic is listed below as Algorithm 4.

Algorithm 4 SMD dynamic

1:$p_{0}$ — starting prices values, $N$ — number of iterations, $C$ —
parameter to control step size

2:for $t=0,\,1,\,\ldots,\,N-1$ do

3: $i\sim\mathcal{U}\\{1,...,S+D\\}$

4: $\displaystyle
p_{t+1}=\left[p_{t}-\frac{CR}{M\sqrt{t+1}}\widetilde{\nabla}f_{i}(p_{t})\right]_{+}$

5:end for

6:$\displaystyle\widetilde{p}_{N}=\frac{1}{N}\sum_{t=1}^{N}p_{t}$

7:return $\widetilde{p}_{N}$

Theoretical guarantee below utilizes the Lipschitz continuity of function $f$,
more precisely the boundedness of the stochastic gradient norm. Considering
the expression (4), we can bound the first term by applying some economic
reasoning, namely that the supply is always limited or scarcity principle [1].
The second term is also bounded since all $x_{d}$ are probability vectors and
belong to the standard simplex. Stochastic gradient given by (6) is bounded by
the same reason.

###### Theorem 5.1

(Proposition 1 [6]) Let $\|p_{t}-p_{0}\|_{2}\leq R$ and
$\mathbb{E}[\|\widetilde{\nabla}f(p_{t})\|^{2}]\leq M^{2}$ for all
$t\in\\{1,...,N\\}$. The suboptimality of prices $\widetilde{p}_{N}$ given by
SMD dynamic (Algorithm 4) is decreasing as follows

$\mathbb{E}[f(\widetilde{p}_{N})]-f(p_{*})\leq\frac{\max\\{C,C^{-1}\\}RM}{\sqrt{N}}.$

Moreover, to obtain the prices satisfying the suboptimality bound

$\mathbb{E}[f(\widetilde{p}_{N})]-f(p_{*})\leq\varepsilon,$

it is sufficient to call $\widetilde{\nabla}f_{i}$ oracle
$\displaystyle\mathcal{O}\left(\frac{1}{\varepsilon^{2}}\right)$ times.

So, the described SMD dynamic has convergence rate of
$\displaystyle\mathcal{O}\left(N^{-1/2}\right)$ in non-smooth case that takes
place if $\Gamma_{s}=0$ for some $s$. This bound matches that of the SGD
dynamic from Theorem 3.1.

###  5.2 Dual smoothing

The special structure of the function $\pi_{s}(p)$ in (2) allows us to use the
Nesterov’s smoothing technique [12]. The idea is, in the case when
$\Gamma_{s}=0$, to replace the $-\Gamma_{s}\cdot\|y-\hat{y}_{s}\|_{2}^{2}$
term in (2) with the synthetic penalty $-\eta\cdot\|y-y_{0}\|_{2}^{2}$, where
$\eta=\varepsilon/(2R^{2})$ for some fixed target suboptimality $\varepsilon$
and $R$ such that $\|y_{*}-y_{0}\|_{2}\leq R$. With this substitution for all
$s$, following the argumentation of Lemma 1, we have that the modified
function $f_{\eta}$ has Lipschitz-continuous gradient with constant

$L=\frac{2SR^{2}}{\varepsilon}+\sum_{d=1}^{D}\frac{1}{\min_{j}\mu_{j}}.$

At the same time, if we minimize $f_{\eta}$ up to accuracy
$\frac{\varepsilon}{2}$, i.e. we have for some $\widetilde{p}$ that

$f_{\eta}(\widetilde{p})-\min_{p\in\mathbb{R}^{n}_{+}}f_{\eta}(p)\leq\frac{\varepsilon}{2},$

then it holds that

$f(\widetilde{p})-f(p_{*})\leq\varepsilon.$

It means that we can obtain the solution of the problem (1) satisfying the
target suboptimality bound by optimizing the modified function $f_{\eta}$ that
is smooth and therefore allows us to apply some of the methods described in
previous sections.

Note that the transition from non-smooth setting to the smooth one with the
described approach is not free in terms of convergence rate due to the
dependence of the constant $L$ on the target suboptimality $\varepsilon$.

##  6 Numerical experiments

In this section, we focus on the problem (1), which is motivated by important
applications to offering smart pricing options by online marketplaces and
management of demand and supply by such financial intermediaries like brokers
[9].

We provide the results of our numerical experiments, which are performed on a
PC with processor Intel Core i7-8650U 1.9 GHz using pure Python 3.7.3 (without
C code) under managing OS Windows 10 (64-bits). Numpy.float128 data type with
precision $1e-18$ and with max element $\approx 1.19e+4932$ is used. Random
seed is set to 17.

We compare Algorithm 1 and Algorithm 2 with pricing dynamics 4.2 and 4.4 from
[9] on the problem (1) with the following settings: number of suppliers $S=5$,
number of consumers $D=10$, number of products $n=20$, number of groups $m=5$.
$\Gamma_{s}$ is chosen like $10^{-4}$ (some small value that affects to
convergence of pricing schemes from [9]). We generated also $S$ vectors
$\hat{y}_{s}$ from uniform distribution $\mathcal{U}[0.01,2]$ of size $n$.
Initialization point $p$ is chosen from uniform distribution
$\mathcal{U}[0.01,5]$ of size $n$, but it is scaled by the maximum element.
Matrix $A$ consists of columns with each one is drawn from uniform
distribution And $\mu$ values are drawn i.i.d. from uniform distribution
$\mathcal{U}[0.1,1]$ of size $m$.

We consider $c_{s}(y)$ equal to $||y||_{2}^{2}$. Hence, the closed form
solution of (2) is

$y_{s}(p)=\frac{p+2\Gamma_{s}\hat{y}_{s}}{2(1+\Gamma_{s})}.$

To estimate the suboptimality we perform pricing dynamic 4.2 [9] (Gradient
descent method) with stopping criterion $||p_{t+1}-p_{t}||_{2}\leq 10^{-10}$
and obtain lower bound for $f_{*}=f(p_{*})$ from Theorem 4.3 [9].

During the experiments we store objective suboptimality and number of oracle
calls for each algorithm at every iteration. The results are presented on the
Figure 1.

Figure 1: Dependence of $f-f_{*}$ from the number of oracle calls.

As we can see from the Figure 1, there is almost no difference between
accelerated and non-accelerated pricing schemes from [9] and that the
stochastic algorithms converge to the optimal value much faster than
accelerated gradient method with respect to chosen parameters.

## Conclusion

We propose a stochastic version of the formulation of the problem of finding
equilibrium prices by minimizing the potential function of total expected
revenue, proposed in a deterministic form in [9]. Thanks to the introduction
of the stochastic oracle, we analyze and interpret in terms of observing the
real-life marketplace environment several dynamics based on the efficient
stochastic gradient optimization methods. We also analyze the case of a non-
smooth potential function in the case of zero quantity adjustment costs. In
addition, we propose a generalized setting allowing an infinite number of
market participants and equally suitable application of optimization methods.
Numerical experiments show that stochastic methods turn out to be discernibly
more efficient in comparison with full-gradient dynamics.

## References

  * [1] Burke, E.: Thoughts and details on scarcity: originally presented to the Right Hon. William Pitt, in the month of November, 1795. No. 17937, T. Gillet (1800)
  * [2] Danskin, J.M.: The theory of max-min and its application to weapons allocation problems, vol. 5. Springer Science & Business Media (2012)
  * [3] Duchi, J., Hazan, E., Singer, Y.: Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning research 12(7) (2011)
  * [4] Duchi, J.C.: Introductory lectures on stochastic optimization. The mathematics of data 25, 99 (2018)
  * [5] Duchi, J.C., Agarwal, A., Johansson, M., Jordan, M.I.: Ergodic mirror descent. SIAM Journal on Optimization 22(4), 1549–1578 (2012)
  * [6] Lan, G., Nemirovski, A., Shapiro, A.: Validation analysis of mirror descent stochastic approximation method. Mathematical programming 134(2), 425–458 (2012)
  * [7] McFadden, D.: Modeling the choice of residential location. spatial interaction theory and planning models, a. karlquist, l. lundquist, f. snickbars, and jw weibull, eds. North-Holland, Amsterdam 1(9), 7 (1978)
  * [8] Moulines, E., Bach, F.: Non-asymptotic analysis of stochastic approximation algorithms for machine learning. Advances in neural information processing systems 24, 451–459 (2011)
  * [9] Müller, D., Nesterov, Y., Shikhman, V.: Dynamic pricing under nested logit demand. arXiv preprint arXiv:2101.04486 (2021)
  * [10] Nemirovski, A., Juditsky, A., Lan, G., Shapiro, A.: Robust stochastic approximation approach to stochastic programming. SIAM Journal on Optimization 19(4), 1574–1609 (2009). https://doi.org/10.1137/070704277
  * [11] Nemirovskij, A.S., Yudin, D.B.: Problem complexity and method efficiency in optimization. Wiley-Interscience (1983)
  * [12] Nesterov, Y.: Smooth minimization of non-smooth functions. Mathematical programming 103(1), 127–152 (2005)
  * [13] Nesterov, Y., Shikhman, V.: Distributed price adjustment based on convex analysis. Journal of Optimization Theory and Applications 172(2), 594–622 (2017)