††thanks: These authors contributed equally to the work††thanks: These authors
contributed equally to the work††thanks: Corresponding author

# Tensorial and bipartite block models
for link prediction in layered networks and temporal networks

Marc Tarrés-Deulofeu<EMAIL_ADDRESS>Departament d’Enginyeria Química,
Universitat Rovira i Virgili, 43006 Tarragona, Catalonia  Antonia Godoy-Lorite
<EMAIL_ADDRESS>Department of Mathematics, Imperial College London,
London SW7 2AZ, United Kingdom  Roger Guimerà<EMAIL_ADDRESS>Departament d’Enginyeria Química, Universitat Rovira i Virgili, 43006
Tarragona, Catalonia  ICREA, 08010 Barcelona, Catalonia  Marta Sales-Pardo
<EMAIL_ADDRESS>Departament d’Enginyeria Química, Universitat Rovira i
Virgili, 43006 Tarragona, Catalonia

###### Abstract

Many real-world complex systems are well represented as multilayer networks;
predicting interactions in those systems is one of the most pressing problems
in predictive network science. To address this challenge, we introduce two
stochastic block models for multilayer and temporal networks; one of them uses
nodes as its fundamental unit, whereas the other focuses on links. We also
develop scalable algorithms for inferring the parameters of these models.
Because our models describe all layers simultaneously, our approach takes full
advantage of the information contained in the whole network when making
predictions about any particular layer. We illustrate the potential of our
approach by analyzing two empirical datasets—a temporal network of email
communications, and a network of drug interactions for treating different
cancer types. We find that modeling all layers simultaneously does result, in
general, in more accurate link prediction. However, the most predictive model
depends on the dataset under consideration; whereas the node-based model is
more appropriate for predicting drug interactions, the link-based model is
more appropriate for predicting email communication.

###### pacs:

##  I Introduction

Imagine a team of researchers looking for promising drug combinations to treat
a specific cancer type for which current treatments are ineffective. The team
has data on the effect of certain pairs of drugs on other cancer types, but
the data are very sparse—only a few drug pairs have been tested on each cancer
type, and each drug pair is tested in a few cancer types, at best, or has
never been tested at all. The challenge is to select the most promising drug
pairs for testing with the target cancer type, so as to minimize the cost
associated to unsuccessful tests.

We can formalize this challenge as the following inference problem: We have a
partial observation of the pairwise interactions between a set of nodes
(drugs) in different “network layers” (cancer types), and we need to infer
which are the unobserved interactions within each layer (drug interactions in
each cancer type). This challenge is relevant for the many systems that can be
represented as multilayer networks De Domenico  _et al._ (2013); Bassett and
Sporns (2017); Kéfi  _et al._ (2016); Pilosof  _et al._ (2017), and is also
formally analogous to the challenge of predicting the existence of
interactions between nodes in time-resolved networks Iribarren and Moro
(2009); Mucha  _et al._ (2010); Gauvin  _et al._ (2014); Peixoto (2015);
Ghasemian  _et al._ (2016); Sapienza  _et al._ (2017); Li  _et al._ (2017).
For instance, we would face the same situation if we had data about the daily
e-mail or phone communications between users, and wanted to infer the
existence of interactions between pairs of users on a certain unobserved day;
in this case each layer would be a different day.

Here, we introduce new generative models that are suitable to address the
challenge above. We model all layers concurrently, so that our approach takes
full advantage of the information contained in all layers to make predictions
for any one of them. Our approach relies on the fact that having information
on the interactions in different layers aids the inference process; in other
words, that the interactions in layers different from the one we are
interested in are informative about the interactions in the query layer. For
instance, biologically similar cancer types are likely to show similar
responses to the same drug pairs, and similar days of the week (for instance
weekdays versus weekends) are also likely to display similar communication
patterns for pairs of users.

Our approach is based on recent results on probabilistic inference on
stochastic block models, which has been successful at modeling the structure
of complex networks Guimerà and Sales-Pardo (2009); Peixoto (2014); Vallès-
Català  _et al._ and at predicting the behavior in biological Guimerà and
Sales-Pardo (2013) and social Rovira-Asenjo  _et al._ (2013); Godoy-Lorite
_et al._ (2016a) systems. In particular, we focus on mixed-membership
stochastic block models Airoldi  _et al._ (2008), in which nodes are allowed
to belong to multiple groups simultaneously. With these models it possible to
model large complex networks with millions of links and, because they are more
expressive than their fixed-membership counterparts, their predictive power is
often superior Godoy-Lorite  _et al._ (2016a). We propose two different mixed-
membership multi-layer network models—a tensorial model that takes nodes as
the basic unit to describe interactions in different layers, and a bipartite
model that takes links (or pairs of nodes) as the basic unit. In our models,
layers, as well as nodes or links, are grouped based on the similarities among
the interaction patterns observed in them. This is in contrast to existing
approaches, which do not take full advantage of the information that each
layer carries about the structure of some other layers.

We illustrate our models and inference approaches by analyzing two datasets—a
network of drug interactions in different cancer types, and a temporal network
of email communications Godoy-Lorite  _et al._ (2016b). We find that modeling
all layers simultaneously, and assuming that they can be grouped, results in
link predictions that are more accurate (in terms of standard metrics such as
precision and recall) than those of single-layer models and of simpler
multilayer models. However, the most predictive model (node-based or link-
based) depends on the dataset under consideration. Indeed, whereas for drug
interactions drug groups are very informative and, therefore, node-based
models are most predictive, temporal email networks are best described in
terms of links, that is, in terms of the relationships between pairs of
individuals rather than the individuals themselves.

##  II Tensorial and bipartite mixed-membership block models for layered
networks

We aim to model $N$ nodes interacting by pairs in $M$ different layers; these
layers correspond to the different contexts in which the nodes interact (for
example, different cancer types or time windows). We represent these
interactions as a layered graph $G$ whose links $(i,j,\ell)$ represent
interactions between nodes $i$ and $j$ in layer (or at time) $\ell$. Moreover,
we allow for multi-valued interactions so that $(i,j,\ell)$ can be of
different types $r_{ij\ell}\in R$, where $R$ is a finite set. Note that we can
use this formalism to model labels, attributes or ratings associated to the
interactions Guimerà and Sales-Pardo (2013); Godoy-Lorite  _et al._ (2016a);
graphs with binary interactions are therefore a particular case within this
general framework in which $r_{ij\ell}=1$ if the interaction occurs and
$r_{ij\ell}=0$ if it does not.

We consider two types of generative models—one that takes individual nodes as
its basic unit, and one that models links (or node pairs). The first
generative model, based on individual nodes, is as follows. There are $K$
groups of nodes and $L$ groups of layers. We assume that the probability that
a node in group $\alpha$ has an interaction of type $r$ with a node in group
$\beta$ in a layer in group $\gamma$ is $p_{\alpha\beta\gamma}(r)$.
Furthermore, we assume that both nodes and layers can belong to more than one
group. To model such mixed group memberships Airoldi  _et al._ (2008), to each
node $i$ we assign a vector $\theta_{i}\in\mathbb{R}^{K}$, where
$\theta_{i\alpha}\in[0,1]$ denotes the probability that node $i$ belongs to
group $\alpha$. Similarly, to each layer $\ell$ we assign a vector
$\eta_{\ell\gamma}\in\mathbb{R}^{L}$. These vectors are normalized so that
$\sum_{\alpha}\theta_{i\alpha}=\sum_{\gamma}\eta_{\ell\gamma}=1$. The
probability that link $(i,j,\ell)$ is of type $r$ is then

$\Pr[r_{ij\ell}=r]=\sum_{\alpha\beta\gamma}\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r)\,.$
(1)

Note that if link types are exclusive (i.e. each edge can be of only one
type), the probability tensor must satisfy the constraint $\sum_{r\in
R}p_{\alpha\beta\gamma}(r)=1$. Since this model is an extension of the mixed-
membership stochastic block model Airoldi  _et al._ (2008); Godoy-Lorite  _et
al._ (2016a) where the probability matrices become tensors because of the
multiple layers De Domenico  _et al._ (2013), we call it the tensorial mixed-
membership stochastic block model (T-MBM).

Our second generative model for layered networks is as follows. Instead of
assuming that nodes belong to groups, we assume that it is links (or pairs of
nodes, rather than individual nodes) that belong to groups Peixoto (2015). In
this model we have $J$ groups of links, and the probability that a link
$e_{ij}\equiv e$ in group $\alpha$ is of type $r$ in a layer $\ell$ in group
$\gamma$ is $p_{\alpha\gamma}(r)$. We also assume that links can belong to
more than one group so that $\zeta_{e\alpha}$ is the probability that link $e$
belongs to group $\alpha$ and $\sum_{\alpha}\zeta_{e\alpha}=1$. As before, to
each layer $\ell$ we also assign a vector $\eta_{\ell}\in\mathbb{R}^{L}$ of
group memberships. Then, the probability that a given link in a particular
layer is of type $r$ is

$\Pr[r_{ij\ell}=r]=\Pr[r_{e\ell}=r]=\sum_{\alpha\gamma}\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r)\,,$
(2)

where, as before, if link types are exclusive the probability matrices satisfy
the condition $\sum_{r\in R}p_{\alpha\gamma}(r)=1$. This model can be seen as
a bipartite model with two types of elements, links and layers. In this
representation, a link $e_{ij}$ has a connection of type $r$ to a layer $\ell$
if $r_{ij\ell}=r_{e\ell}=r$. Therefore, we call this model the bipartite
mixed-membership stochastic block model (B-MBM).

These models are novel in a number of ways. First, unlike other models of
multi-layer networks Peixoto (2015); Bacco  _et al._ (2017), they do not
assume any particular order in the layers, and therefore do not impose any
restrictions to how layers should be grouped. This is in contrast to
approaches for temporal networks that can only group layers corresponding to
consecutive times. While such restriction simplifies the task of grouping
layers, it also eliminates the possibility of identifying, for example,
periodicities in temporal networks. More importantly, this restriction
prevents models from being applicable to non-temporal multilayer networks. Our
models eliminate this restriction.

Second, unlike other models Stanley  _et al._ (2016); Bacco  _et al._ (2017),
ours assume that group memberships (of nodes or links) do not change from
layer to layer. Rather, we argue that in many relevant situations membership
is determined by intrinsic properties of nodes or links and so long as these
properties do not change, group membership should not change either. For
example, membership of individuals to groups in social networks is related to
demographic and socio-economic characteristics, which are unlikely to change
in periods of months or even a few years. Or in drug-interaction networks,
membership of drugs to groups is related to the mechanism of action and the
targets of the drug Guimerà and Sales-Pardo (2013), which do not change
regardless of the situation in which the drug is used.

Third, our models naturally deal with situations in which links have metadata,
that is, situations in which nodes are not only connected or disconnected, but
rather can be connected with links of different types.

Finally, unlike other models of multilayer and temporal networks Peixoto
(2015); Bacco  _et al._ (2017); Stanley  _et al._ (2016), in our models
nodes/links and layers do not belong to a single group, but rather to a
mixture of groups Airoldi  _et al._ (2008); Godoy-Lorite  _et al._ (2016a).
This allows us to develop efficient expectation-maximization algorithms that
can be massively parallelized Jefries (2017) and, at the same time, provide
better predictions than single-group models Godoy-Lorite  _et al._ (2016a).

##  III Inference equations and expectation-maximization algorithms

Given a set $G^{O}$ of observed links types, our goal is to predict the types
$r_{ij\ell}$ of links $(i,j,\ell)\not\in G^{O}$ whose type is unknown. Because
marginalizing over the parameters in our models (Eqs. (1) and (2)) is too
time-consuming, here we present a maximum likelihood approach (and the
corresponding expectation-maximization algorithms) for the two models above.

###  III.1 Tensorial model

Given the generative T-MBM model in Eq. (1), and abbreviating its parameters
as $\boldsymbol{\theta},\boldsymbol{\eta},\mathbf{p}$, the likelihood of the
model is

$P(G^{O}|\boldsymbol{\theta},\boldsymbol{\eta},\mathbf{p})=\prod_{(ij\ell)\in
G^{O}}\sum_{\alpha\beta\gamma}\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r_{ij\ell})\,.$
(3)

As we show below (Appendix A), the values of the parameters that maximize this
likelihood satisfy the following equations

$\displaystyle\theta_{i\alpha}$ $\displaystyle=$
$\displaystyle\frac{\sum_{(j\ell)\in\partial
i}\sum_{\beta\gamma}\omega_{ij\ell}(\alpha,\beta,\gamma)}{d_{i}}\,,$ (4)
$\displaystyle\eta_{\ell\gamma}$ $\displaystyle=$
$\displaystyle\frac{\sum_{(ij)\in\partial\ell}\sum_{\alpha\beta}\omega_{ij\ell}(\alpha,\beta,\gamma)}{d_{\ell}}\,,$
(5) $\displaystyle p_{\alpha\beta\gamma}(r)$ $\displaystyle=$
$\displaystyle\frac{\sum_{(i,j,\ell)\in
G^{O}|r_{ij\ell}=r}\omega_{ij\ell}(\alpha,\beta,\gamma)}{\sum_{(i,j,\ell)\in
G^{O}}\omega_{ij\ell}(\alpha,\beta,\gamma)}\,.$ (6)

Here, $\partial i=\\{(j,\ell)|(i,j,\ell)\in G^{O}\\}$ are the set of observed
layer-specific neighbors of node $i$ and $d_{i}=|\partial i|$ is the total
degree of the node in all the layers. Similarly,
$\partial\ell=\\{(i,j)|(i,j,\ell)\in G^{O}\\}$ is the set of observed links in
layer $\ell$ and $d_{\ell}=|\partial\ell|$. Finally,
$\omega_{ij\ell}(\alpha,\beta,\gamma)$ is the estimated probability that the
type of a given link $r_{ij\ell}$ is due to $i$, $j$ and $\ell$ belonging to
groups $\alpha$, $\beta$ and $\gamma$ respectively, and is given by

$\omega_{ij\ell}(\alpha,\beta,\gamma)=\frac{\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r_{ij\ell})}{\sum_{\alpha^{\prime}\beta^{\prime}\gamma^{\prime}}\theta_{i\alpha^{\prime}}\theta_{j\beta^{\prime}}\eta_{\ell\gamma^{\prime}}p_{\alpha^{\prime}\beta^{\prime}\gamma^{\prime}}(r_{ij\ell})}\,.$
(7)

These equations can be solved iteratively with an expectation-maximization
algorithm, starting with an initial estimate of $\boldsymbol{\theta}$,
$\boldsymbol{\eta}$, and $\mathbf{p}$ and, then, repeating the following
steps: (i) use Eq. (7) to compute $\omega_{ij\ell}(\alpha,\beta,\gamma)$ for
$(i,j,\ell)\in G^{O}$ (expectation step); (ii) use Eqs. (4)-(6) to compute
$\boldsymbol{\theta}$, $\boldsymbol{\eta}$, and $\mathbf{p}$ (maximization
step).

###  III.2 Bipartite model

Similarly, the likelihood of the B-MBM is

$P(G^{O}|\boldsymbol{\zeta},\boldsymbol{\eta},\mathbf{p})=\prod_{(e,\ell)\in
G^{O}}\sum_{\alpha\gamma}\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r_{e\ell})\,,$
(8)

and the maximum likelihood estimators of the parameters satisfy

$\displaystyle\zeta_{e\alpha}$ $\displaystyle=$
$\displaystyle\frac{\sum_{\ell\in\partial
e}\sum_{\beta\gamma}\phi_{e\ell}(\alpha,\gamma)}{d_{e}}\,,$ (9)
$\displaystyle\eta_{\ell\gamma}$ $\displaystyle=$
$\displaystyle\frac{\sum_{e\in\partial\ell}\sum_{\alpha}\phi_{e\ell}(\alpha,\gamma)}{d_{\ell}}\,,$
(10) $\displaystyle p_{\alpha\gamma}(r)$ $\displaystyle=$
$\displaystyle\frac{\sum_{(e\ell)\in
G^{O}|r_{e\ell}=r}\phi_{e\ell}(\alpha,\gamma)}{\sum_{(e\ell)\in
G^{O}}\phi_{e\ell}(\alpha,\gamma)}\,.$ (11)

Here $\partial e=\\{\ell|(e,\ell)\in G^{O}\\}$ are the observations of link
$e_{ij}$ in all layers and $d_{e}=|\partial e|$. As before,
$\partial\ell=\\{e|(e,\ell)\in G^{O}\\}$ are the observed links in layer
$\ell$ and $d_{\ell}=|\partial\ell|$. Finally, $\phi_{e\ell}(\alpha,\gamma)$
is the estimated probability that the type of a specific link $r_{e\ell}$ is
due to $e$ and $\ell$ belonging to groups $\alpha$ and $\gamma$ respectively;
we can compute $\phi_{e\ell}(\alpha,\gamma)$ as

$\phi_{e\ell}(\alpha,\gamma)=\frac{\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r_{e\ell})}{\sum_{\alpha^{\prime}\gamma^{\prime}}\zeta_{e\alpha^{\prime}}\eta_{\ell\gamma^{\prime}}p_{\alpha^{\prime}\gamma^{\prime}}(r_{e\ell})}\,.$
(12)

Like in the tensorial model, these equations can be solved iteratively using
an expectation-maximization algorithm.

##  IV Model comparison on real data

###  IV.1 Datasets

We perform experiments on two different datasets: the time-resolved email
network of an organization spanning one year Godoy-Lorite  _et al._ (2016b),
and a network of drug-drug interactions in different cancer cell lines Menden
_et al._ (2018). In the email dataset, we represent each day as a different
layer of the multi-layer network, and two users are considered to interact in
a given day if they send at least one email in either direction during that
day. We consider several e-mail networks that correspond to e-mail
communications within organizational units (see Table 1).

In the drug-drug interactions dataset, each layer corresponds to a different
cancer cell line and we have information on the effects of some drug pair
combinations on some cancer cell lines Menden  _et al._ (2018). In contrast to
the email dataset, in which all the interactions (or lack of interaction) are
observed, this dataset is sparsely observed—we have information about 1.5% of
the drug pairs. Specifically, the available experimental data is a real-valued
magnitude representing the combined efficiency of two drugs on a particular
cell line. These magnitudes range from large absolute values, in which case
the interaction is said to be synergistic (if it is positive) or antagonistic
(if it is negative), to small absolute values, in which case the interaction
it is said to be additive. In an additive interaction, the application of the
two drugs together has an efficiency equal or similar to the sum of the
efficiencies of each drug administered separately. By contrast, in a
synergistic (antagonistic) interaction the efficiency of the two drugs
administered together is significantly higher (lower) than the sum of the
efficiencies of each drug administered separately.

In Table 1 we show the characteristics of each dataset in terms of the types
of links $R$, the total number of nodes, the total number of layers, the total
number of possible links, and the number and fraction of actually observed
links. In all cases, we fitted and validated our models using a 5-fold cross-
validation scheme (Appendix C).

Table 1: Dataset characteristics. The email networks we consider are complete networks where no-links are treated as links of type $0$, so all potential links in the network are observed links. The drug-drug interaction network is a sparse dataset, where we only have information about $1.4\%$ of the links. Each observed drug-drug interaction can be of three types: antagonistic - ANT, additive - ADD, or synergistic - SYN. Dataset | Types of links $R$ | #Nodes | #Layers | #Observables | Fraction observed | #Observed
---|---|---|---|---|---|---
Email Unit 1 | $\\{0,1\\}$ | 104 | 365 | 1,954,940 | 100% |  $|G^{O}|_{1}$ = 20,807
Email Unit 2 | $\\{0,1\\}$ | 114 | 365 | 2,350,965 | 100% |  $|G^{O}|_{1}$ = 27,180
Email Unit 3 | $\\{0,1\\}$ | 116 | 365 | 2,434,550 | 100% |  $|G^{O}|_{1}$ = 23,979
Email Unit 4 | $\\{0,1\\}$ | 118 | 365 | 2,519,595 | 100% |  $|G^{O}|_{1}$ = 17,508
Email Unit 5 | $\\{0,1\\}$ | 141 | 365 | 3,602,550 | 100% |  $|G^{O}|_{1}$ = 23,923
Email Unit 6 | $\\{0,1\\}$ | 161 | 365 | 4,701,200 | 100% |  $|G^{O}|_{1}$ = 20,790
Email Unit 7 | $\\{0,1\\}$ | 225 | 365 | 9,198,000 | 100% |  $|G^{O}|_{1}$ = 60,238
Drug-drug interactions | {ANT, ADD, SYN} | 69 | 85 | 199,410 | 1.37% |  $|G^{O}|_{ANT}$ = 385
Drug-drug interactions | {ANT, non-ANT} | 69 | 85 | 199,410 | 1.37% |  $|G^{O}|_{ADD}$ = 1,543
Drug-drug interactions | {SYN, non-SYN} | 69 | 85 | 199,410 | 1.37% |  $|G^{O}|_{SYN}$ = 863

###  IV.2 Baseline models

We compare our models to three different baselines. The naive baseline takes
into account all the observations of a link $(i,j)$ in the training set. Then,
it makes predictions for the unobserved link types $r_{ijl}$ based on the
fraction of times link $(i,j)$ has been observed to be of type $r$ in the
training set

$P_{\rm naive}[r_{ijl}=r]=\frac{\sum_{s|(ijs)\in
G^{O}}\delta{r,r_{ijs}}}{N_{ij}}~{},$ (13)

where $N_{ij}$ is the number of times the $(i,j)$ is observed in the training
set.

The other two baselines help us to assess to what extent our models are able
to exploit correlations between different layers. First, the independent-layer
naive estimates the probability of link $(i,j,l)$ being of type $r$ as the
fraction of links of type $r$ observed in layer $l$,

$P_{\rm naive-IL}[r_{ijl}=r]=\frac{\sum_{(kn)|(knl)\in
G^{O}}\delta{r,r_{knl}}}{N_{l}}~{},$ (14)

where $N_{l}$ is the number of links of any type observed in layer $l$ of the
training set. Second, we consider an independent-layer mixed-membership
stochastic block model for each layer so that,

$P_{\rm
IL}[r_{ijl}=r]=\sum_{\alpha\beta}\theta_{i\alpha}^{l}\theta_{j\beta}^{l}p_{\alpha\beta}^{l}(r)~{},$
(15)

where the superindex $l$ denotes that each layer has its own set of
parameters. As in the tensorial and bipartite layered models, parameters are
subject to the constraints $\sum_{\alpha}\theta^{l}_{i\alpha}=1,\,\forall l$
and $\sum_{r}p^{l}_{\alpha\beta}(r)=1\,\forall l$. The parameters for this
model are obtained using the same method as in the tensorial and bipartite
mixed-membership models, but considering each layer separately (see also
Godoy-Lorite  _et al._ (2016a)).

###  IV.3 Email networks

We first consider the ability of each model and baseline to predict unobserved
links in the email networks listed in Table 1. To assess the performance of
each model for each network, we calculate the area under the ROC curve (AUC),
the precision, and the recall in 5-fold cross-validation experiments (see
Appendix C for details). The AUC measures how well a model separates active
links (type-1, for which there is communication between the individuals) from
inactive links (type-0, with no communication). In particular, it measures the
frequency with which an active unobserved link is assigned a higher
probability to be active than an inactive unobserved link. Precision accounts
for the fraction of links predicted to be active that are indeed active.
Recall gives the fraction of active links that are predicted to be active. To
calculate both precision and recall, we need to set a threshold $T$ that
allows to map probabilities $P[r_{ijl}=1]$ into a binary variable. We do it as
follows: if $P[r_{ijl}=1]\geq T$ then the model predicts that $r_{ijl}=1$,
otherwise it predicts $r_{ijl}=0$. In what follows, we choose $T$ as the
density of active links in the training set. The reasoning behind this
decision is that, because we are splitting data at random into a training and
a test set, the test set should have a fraction of active links close to that
of the training set. 111This particular choice of threshold leads, when models
are properly calibrated in a frequentist sense, to both precision and recall
having very similar values (see Supplementary Information).

In Fig. 1, we show that the bipartite link-based model outperforms the
tensorial and baseline models in all metrics (see Fig. S1 for all other email
units).

Figure 1: Predictive performance of the models for email networks. Top: AUC;
Middle: Precision; Bottom: Recall. Each bar represents the average of the
5-fold cross-validation for a given model (see Appendix C). The error bars
(shown as a vertical line, which is small and not visible in some cases)
represent the standard error of the mean.

In these email networks, the AUC is quite high even for the naive baseline
because most pairs of individuals never exchange an e-mail and therefore it is
easy to predict links for which $r_{ijl}=0$ in all observed layers. The
situation changes when we look at precision and recall, which clearly show
that the bipartite model is consistently and significantly superior at
predicting links that are active. Somewhat surprisingly, we also find that the
tensorial node-based model gives slightly lower values than the naive baseline
model. The explanation lies in the fact that, contrary to both the naive
baseline and the bipartite models, the tensorial model focuses primarily on
nodes rather than on links and is thus less likely to account for the fact
that many pairs of nodes in the network never communicate. More precisely, the
probabilities assigned by the tensorial model depend on the product of the
membership of the involved nodes, and these memberships are rarely equal to
zero. Hence, according to the tensorial model most links have a non-zero
probability of existing, including those that are inactive for all
observations in the training set.

To further investigate the workings of each approach, we analyze whether they
are properly calibrated in a frequentist sense, that is whether the fitted
models are able to reproduce statistical features of the training dataset
Gneiting  _et al._ (2007). In particular, we consider the marginal and
probabilistic calibration of all models. A model is probabilistically
calibrated if events to which the model assigns a probability $p$ are observed
with frequency $p$ Gneiting  _et al._ (2007). In our case, a model is
calibrated if a fraction $p$ of the links for which $P[r_{ijl}=1]=p$ actually
exist. A model is marginally calibrated if, on average, each type of event is
assigned a probability that is equal to the actual frequency of such events in
the training set. In our case, a model is calibrated if the mean
$P[r_{ijl}=1]$ assigned to links coincides with the density of the observed
network Gneiting  _et al._ (2007).

In Fig. 2, we show that all models are relatively well calibrated
probabilistically (higher probabilities correspond to higher frequencies),
although the calibration is noticeably worse for the network obtained for Unit
1 (see Fig. S2 for the remaining units). In general, the bipartite model is
better calibrated than the tensorial model, which is consistent with the
higher predictive accuracy of the bipartite, link-based model. Perhaps
surprisingly, the naive baseline model appears to have an even better
probabilistic calibration across all units. Figure 3 also shows that all
models are marginally calibrated.

In light of these observations, the difference in performance between
bipartite and naive models must come from the fact that the bipartite model is
able to detect temporal patterns that are relevant for the prediction of
active links. Indeed, we find that for all the email networks we consider,
temporal layers (days) are classified either as week days or as weekend days
(and holidays), so that it is more likely for any link to be active on a week
day. Interestingly, this is all the temporal information required to be able
to accurately predict whether a specific link is going to be active or not on
a certain day 222Note that our results do not depend on the number of latent
dimensions allowed for the temporal layers L, since for $L>2$ we also find
that temporal layers have $\eta_{\ell\gamma}\neq 0$ only for two latent groups
$\gamma$.

Figure 2: Probabilistic calibration of the models for email networks. Each
point in each line represents the average of the 5-fold cross-validation for a
given model, with error bars representing the standard error of the mean. The
Naive I.L. model is not included as it only assigned tiny probabilities that
resulted in a single data point near the origin.

Figure 3: Marginal calibration of the models for email networks. Each line
corresponds to a different model, and each point in a line corresponds to a
different email network (see Table 1). Each point represents the average over
the 5-fold cross-validation for a given model. Error bars are smaller than
symbols.

###  IV.4 Drug-drug interactions in cancer

Links in the drug-drug interaction network are of three different types:
synergistic, antagonistic, and additive; we trained the models considering the
three types of interactions. However, because the interesting question is
whether synergistic or antagonistic interactions can be predicted, we
evaluated the performace of each model for each task of these two tasks. For
instance, to evaluate the accuracy of a model at predicting synergistic
interactions, we binarized model predictions into synergistic and non-
synergistic. We then computed the metrics over this binary outcome as we did
for e-mail networks (Figs. 4, 5, and 6; Fig. S3 shows that all of the results
below are qualitatively similar when training our models on networks with only
two types of interactions: synergistic/non-synergistic or antagonistic/non-
antagonistic) 333Due to the sparsity of observations in these networks many
interactions were never observed in the training sets, and thus no group
memberships could be assigned to the links ($e_{ij}$) corresponding to those
interactions. We solved this cold start problem by, at each iteration,
assigning them the average membership of the observed interactions
$\zeta_{e}={\left\langle\zeta_{f}\right\rangle}_{f\in G^{O}}$. Analogously, if
a node $i$ had no observed interactions in the training set, at each iteration
we set its membership vector as the average of membership vectors for nodes
with observed interactions
$\theta_{i}={\left\langle\theta_{k}\right\rangle}_{k\in G^{O}}$..

Contrary to what we observed for the email networks, we find that the
tensorial model performs better than the bipartite model. Our results thus
suggest that for this dataset, grouping nodes (drugs) into groups summarizes
more parsimoniously the information relevant for prediction. This is
consistent with previous findings that show that mechanisms of action and
target pathways of drugs are related to the effect they display when combined
with other drugs, an information that is best captured by node memberships
than by link memberships Guimerà and Sales-Pardo (2013).

Interestingly, we observe differences in performance at detecting antagonistic
and synergistic interactions. For the synergistic interaction network, we find
that the tensorial model consistently outperforms the bipartite and baseline
models in all metrics (AUC, precision and recall), although its marginal
calibration is slightly worse that that of the other models. For antagonistic
interactions, the tensorial model also performs better than the bipartite and
baseline models in terms of AUC. However, the tensorial model has a precision
and recall that are similar to those of the independent-layers baseline model.
The generalized decrease in precision and recall with respect to synergistic
network does not come as a surprise since none of the models is perfectly
calibrated for probabilities lower than the density of the training set (Fig.
5). In fact, we observe that the fraction of antagonistic interactions for
which $P[r_{ijk}=1]<T$ is larger than desired. As a result, some antagonistic
interactions are counted as non-antagonistic interactions in terms of
precision and recall. This effect is exhacerbated by the fact that, due to the
sparsity of the network, a large fraction of interactions are assigned low
probability values by all of the models.

The fact that the independent-layer model has prediction and recall values
similar to those of the tensorial model can be explained by the fact that
antagonistic interactions are more localized to specific layers than
synergistic interactions are (see Fig. S6). This situation makes it easier for
the independent-layer baseline model to make more accurate predictions for
these layers. Note however, that if more information on antagonistic
interactions was available, the performance of the tensorial model would
likely be comparable to that of the synergistic case.

Figure 4: Predictive performance of the models for drug-drug interaction
networks. (a) AUC statistic. (b) Precision. (c) Recall. Each bar represents
the average of the 5-fold cross-validation for a given model, with error bars
(shown as a vertical line for clarity) representing the standard error of the
mean.

Figure 5: Probabilistic calibration of the models for drug-drug interaction
networks. Each point in each line represents the average of the 5-fold cross-
validation for a given model. Error bars represent the standard error of the
mean. The vertical dashed lines show the density of each training set.

Figure 6: Marginal calibration of the models for drug-drug interaction
networks. For each of the models we consider (see legend) we plot the average
probability for links being of a certain type (antagonistic or synergistic)
with respect to the density of links of that type in the training set. Each
point represents the average over the 5 training sets for a given model. Error
bars are smaller than symbols.

##  V Discussion

We have presented two mixed-membership multi-layer network models that can be
applied to any multi-layer networks, with layers representing temporal
snapshots of the interactions or different contexts for the interactions. By
extending the mixed-membership paradigm to the layers themselves, and by not
making any prior assumption about them, our models can detect and take
advantage of inter-layer correlations in the network of interactions to make
better predictions. As a result, both our multi-layer models outperform the
baseline models in almost all the studied cases, except for the cases in which
information is too sparse for the multi-layer model to recover unobserved
interactions with precision.

Importantly, none of the models we present—the tensorial node-based model or
the bipartite link-based model—is intrinsically better than the other; however
they can hold clues as to the mechanisms that are predictive of interaction
types. Our results precisely illustrate this fact. We find that the bipartite
model works better for email networks in which the communication between pairs
of users (links), rather than the users themselves, together with their
temporal evolutions are the relevant description unit for prediction. This
could be due to the fact that we are analyzing communication at a rather small
scale (people working within the same unit of an organization), and it is
possible that a node-based model could be better for communication between
users at a larger scale. Moreover, as the network grows the number of $\theta$
parameters for the tensorial model scales linearly with the number of nodes,
whereas the number of $\eta$ parameters for the bipartite model scales
quadratically. In really big networks it is then plausible that the tensorial
becomes more parsimonious.

Conversely, our results show that for the drug-drug interaction network the
relevant unit of description are drugs (nodes). This is consistent with the
fact that the mechanism of action/target that determines how a drug will
interact with another one; this information is encapsulated in the node (and
its observed interactions). The use of the interactions of nodes in different
cancer types (layers) boosts our ability to predict the type of type-dependent
interactions more precisely. In contrast, the description of these networks in
terms of drug-pair interactions completely misses the drug-specific
information that is relevant for prediction in this context.

Our results unambiguously show that using the information of the interactions
on other layers helps obtain better models. Remarkably, the flexibility of the
models we propose make this approach suitable to analyze multi-layer networks
in any context. A natural step to further improve the model and prediction
accuracy would be to include auxiliary data (i.e. metadata such as node or
link attributes) into the modeling process. This problem has just started
being explored in the literature Newman and Clauset (2016); Hric  _et al._
(2016); Sapienza  _et al._ (2017), so there is no general framework on how to
introduce auxiliary data into the inference process yet. Nonetheless, recent
results show that single-layer mixed-membership models are suitable models to
incorporate specific types of auxiliary data into the inference process
without adding methodological complexity Cobo  _et al._ (2017), thus opening
the window to developing general inference frameworks that consider different
types of metadata also in multilayer contexts.

###### Acknowledgements.

This work was supported by the Spanish Ministerio de Economia y
Comptetitividad (MINECO) Grant FIS2016-78904-C3-1-P, and European Union FET
Grant 317532 (MULTIPLEX). The authors acknowledge the computer resources at
MinoTauro and the technical support provided by the Barcelona Supercomputing
Center (FI-2015-3-0024, FI-2016-3-0026). The authors acknowledge AstraZeneca
UK Limited, Sanger and Sage Bionetworks-DREAM for providing the cancer drug
interaction data.

##  Appendix A Derivation of the expectation maximization equations for the
T-MBM

In the tensorial mixed-membership stochastic block model, we assign membership
vectors $\theta_{i\alpha}$, $\eta_{\ell\gamma}$ to each node $i$ and each
layer $\ell$, respectively. These membership vectors are properly normalized,
therefore represent the probability that each node/layer belongs to a specific
node/layer group:

$\forall
i:\sum_{\alpha=1}^{K}\theta_{i\alpha}=1\,,\quad\forall\ell:\sum_{\gamma=1}^{L}\eta_{\ell\gamma}=1\,.$
(A.1)

Because we consider that links cant take different values $r\in R$, to ensure
that each observed interaction has probability 1 of receiving any rating, we
normalize probability matrices $p_{\alpha\beta\gamma}(r)$

$\forall\alpha,\beta,\gamma:\sum_{r\in R}p_{\alpha\beta\gamma}(r)=1\,.$ (A.2)

Note that if $R=\\{0,1\\}$, then
$p_{\alpha\beta\gamma}(0)=1-p_{\alpha\beta\gamma}(1)$.

We maximize the likelihood (3) as a function of
$\boldsymbol{\theta},\boldsymbol{\eta},\mathbf{p}$ using an expectation
maximization (EM) algorithm. We start with a standard variational and use
Jensen’s inequality $\log\bar{x}\geq\overline{\log x}$ in order to transform
the logarithm of a sum into a sum of logarithms

$\displaystyle\log$ $\displaystyle
P(G^{O}|\boldsymbol{\theta},\boldsymbol{\eta},\mathbf{p})=\sum_{(ij\ell)\in
G^{O}}\log\sum_{\alpha\beta\gamma}\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r_{ij\ell})$
$\displaystyle=\sum_{(ij\ell)\in
G^{O}}\log\sum_{\alpha\beta\gamma}\omega_{ij\ell}(\alpha\beta\gamma)\,\frac{\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r_{ij\ell})}{\omega_{ij\ell}(\alpha\beta\gamma)}$
$\displaystyle\geq\sum_{(ij\ell)\in
G^{O}}\sum_{\alpha\beta\gamma}\omega_{ij\ell}(\alpha\beta\gamma)\log\frac{\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r_{ij\ell})}{\omega_{ij\ell}(\alpha\beta\gamma)}\,.$
(A.3)

Here we have introduced the auxiliary variable
$\omega_{ij\ell}(\alpha\beta\gamma)$, which is the estimated probability that
a given link’s type $r_{ij\ell}$ is due to $i$, $j$ and $\ell$ belonging to
groups $\alpha$, $\beta$ and $\gamma$ respectively. Note in the expression
above, equality holds when

$\omega_{ij\ell}(\alpha\beta\gamma)=\frac{\theta_{i\alpha}\theta_{j\beta}\eta_{\ell\gamma}p_{\alpha\beta\gamma}(r_{ij\ell})}{\sum_{\alpha^{\prime}\beta^{\prime}\gamma^{\prime}}\theta_{i\alpha^{\prime}}\theta_{j\beta^{\prime}}\eta_{\ell\gamma^{\prime}}p_{\alpha^{\prime}\beta^{\prime}\gamma^{\prime}}(r_{ij\ell})}\,.$
(A.4)

This is precisely the equation for the expectation step.

For the maximization step, we derive update equations for the parameters
$\boldsymbol{\theta},\boldsymbol{\eta},\mathbf{p}$ by taking derivatives of
the log-likelihood (A.3). Including Lagrange multipliers for the normalization
constraints (A.1), we obtain for $\theta_{i\alpha}$

$\theta_{i\alpha}=\frac{\sum_{j\ell\in\partial
i}\sum_{\beta\gamma}\omega_{ij\ell}(\alpha\beta\gamma)}{\sum_{j\ell\in\partial
i}\sum_{\alpha\beta\gamma}\omega_{ij\ell}(\alpha\beta\gamma)}=\frac{\sum_{j\ell\in\partial
i}\sum_{\beta\gamma}\omega_{ij\ell}(\alpha\beta\gamma)}{d_{i}}\,,$ (A.5)

where $\partial i=\\{j,\ell|(ij\ell)\in G^{O}\\}$ and $d_{i}=|\partial i|$ is
the degree of node $i$ in all the layers for any type of link. Similarly, for
$\eta_{\ell\gamma}$ we obtain

$\eta_{\ell\gamma}=\frac{\sum_{ij\in\partial\ell}\sum_{\alpha\beta}\omega_{ij\ell}(\alpha\beta\gamma)}{\sum_{ij\in\partial\ell}\sum_{\alpha\beta\gamma}\omega_{ij\ell}(\alpha\beta\gamma)}=\frac{\sum_{ij\in\partial\ell}\sum_{\alpha\beta}\omega_{ij\ell}(\alpha\beta\gamma)}{d_{\ell}}\,,$
(A.6)

where $\partial\ell=\\{i,j|(ij\ell)\in G^{O}\\}$ and $d_{\ell}=|\partial\ell|$
is the number of observed links of any type in layer $\ell$.

Finally, including a Lagrange multiplier for (A.2), we have for
$p_{\alpha\beta\gamma}(r)$

$p_{\alpha\beta\gamma}(r)=\frac{\sum_{(ij\ell)\in
G^{O}|t_{ij\ell}=t}\omega_{ij\ell}(\alpha\beta\gamma)}{\sum_{(ij\ell)\in
G^{O}}\omega_{ij\ell}(\alpha\beta\gamma)}\,.$ (A.7)

##  Appendix B Derivation of the expectation maximization equations for the
B-MBM

As in the tensorial model, we assign normalized membership vectors
$\zeta_{e\alpha}$, $\eta_{\ell\gamma}$ to links and layers, respectively. We
also consider probability matrices $p_{\alpha\gamma}(r)$ that are as well
normalized $\sum_{r\in R}p_{\alpha\gamma}(r)=1$).

In order to maximize the likelihood, we again use Jensen’s inequality to
transform the the logarithm of a sum into a sum of logarithms and introduce an
auxiliary variable $\phi_{e\ell}(\alpha,\gamma)$:

$\displaystyle\log$ $\displaystyle
P(G^{O}|\boldsymbol{\zeta},\boldsymbol{\eta},\mathbf{p})=\sum_{(ij\ell)\in
G^{O}}\log\sum_{\alpha\gamma}\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r_{e\ell})$
$\displaystyle=\sum_{(ij\ell)\in
G^{O}}\log\sum_{\alpha\gamma}\phi_{e\ell}(\alpha,\gamma)\,\frac{\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r_{e\ell})}{\phi_{e\ell}(\alpha,\gamma)}$
$\displaystyle\geq\sum_{(e\ell)\in
G^{O}}\sum_{\alpha\gamma}\phi_{e\ell}(\alpha,\gamma)\log\frac{\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r_{e\ell})}{\phi_{e\ell}(\alpha,\gamma)}\,.$
(A.8)

where again the equality holds when

$\phi_{e\ell}(\alpha,\gamma)=\frac{\zeta_{e\alpha}\eta_{\ell\gamma}p_{\alpha\gamma}(r_{e\ell})}{\sum_{\alpha^{\prime}\gamma^{\prime}}\zeta_{e\alpha^{\prime}}\eta_{\ell\gamma^{\prime}}p_{\alpha^{\prime}\gamma^{\prime}}(r_{e\ell})}\,,$
(A.9)

giving us the update equation (A.9) for the expectation step.

For the maximization step, we derive update equations for the parameters
$\boldsymbol{\zeta},\boldsymbol{\eta},\mathbf{p}$ by taken derivatives of the
log-likelihood (A.8). Including Lagrange multipliers for the normalization
constraints, we obtain

$\zeta_{e\alpha}=\frac{\sum_{\ell\in\partial
e}\sum_{\gamma}\phi_{e\ell}(\alpha,\gamma)}{\sum_{\ell\in\partial
e}\sum_{\alpha\gamma}\phi_{e\ell}(\alpha,\gamma)}=\frac{\sum_{\ell\in\partial
e}\sum_{\beta\gamma}\phi_{e\ell}(\alpha,\gamma)}{d_{e}}\,,$ (A.10)

where $\partial e=\\{\ell|(e\ell)\in G^{O}\\}$ are the set of layers in which
we observe link $e_{ij}$ and $d_{e}=|\partial e|$ is the total number of
layers in which we observe link $e_{ij}$. Similarly,

$\eta_{\ell\gamma}=\frac{\sum_{e\in\partial\ell}\sum_{\alpha}\phi_{e\ell}(\alpha,\gamma)}{\sum_{e\in\partial\ell}\sum_{\alpha\gamma}\phi_{e\ell}(\alpha,\gamma)}=\frac{\sum_{e\in\partial\ell}\sum_{\alpha}\phi_{e\ell}(\alpha,\gamma)}{d_{\ell}}\,,$
(A.11)

where $\partial\ell=\\{e|(e\ell)\in G^{O}\\}$ and $d_{\ell}=|\partial\ell|$.
Finally, including a Lagrange multiplier for the normalization of
$p_{\alpha\gamma}(r)$, we have

$p_{\alpha\gamma}(r)=\frac{\sum_{(e\ell)\in
G^{O}|r_{e\ell}=r}\phi_{e\ell}(\alpha,\gamma)}{\sum_{(e\ell)\in
G^{O}}\phi_{e\ell}(\alpha,\gamma)}\,.$ (A.12)

Equations (A.9)-(A.12) are solved iteratively with an EM algorithm following
the same procedure as in the tensorial model. The bipartite model also scales
linearly with the size of the dataset, but in this case the number of
parameters of the model is $IK+ML+|G^{O}|K\cdot L$, where the number of links
$I\leq N\cdot(N-1)/2$, thus, even though it increases the number of parameters
(number of nodes $N$ is typically smaller than number of links $I$), there is
one dimension less to run over all observed links in all layers $|G^{O}|$.

##  Appendix C Experimental details

With regards to the drug-drug interactions dataset, we divided the continuous
values of efficiency into three categories (synergistic, additive and
antagonistic) by setting two thresholds as suggested in the original
experimental data. These thresholds are -20.0 and 20.0, so that interactions
with an efficiency lower than -20.0 are classified as antagonistic, those with
an efficiency higher than 20.0 are classified as synergistic, and those in
between are considered additive Menden  _et al._ (2018).

For both datasets, we fitted and validated our models using a 5-fold cross-
validation scheme. We first divided the data into five equal splits. Then for
each fold we considered 4 splits as the training set to which we fitted the
model, and the remaining split was kept as the test set on which we made
predictions. In order to select the number of latent groups $K$, $J$, and $L$,
we used the smallest values for which the prediction accuracy had already
reached saturation values. These values were $K=5$, $L=5$ for the tensorial
model, and $J=2$, $L=2$ for the bipartite model.

For each fold, we repeated the fitting processes between 50 and 100 times with
different random initializations. The results we present correspond to the
average over the results for the five folds.

## References

  * De Domenico  _et al._ (2013) M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas, “Mathematical Formulation of Multilayer Networks,” Phys. Rev. X 3, 041022 (2013).
  * Bassett and Sporns (2017) Danielle S Bassett and Olaf Sporns, “Network neuroscience,” Nature Neuroscience 20, 353–364 (2017).
  * Kéfi  _et al._ (2016) Sonia Kéfi, Vincent Miele, Evie A. Wieters, Sergio A. Navarrete, and Eric L. Berlow, “How structured is the entangled bank? The surprisingly simple organization of multiplex ecological networks leads to increased persistence and resilience,” PLoS Biology 14, e1002527 (2016).
  * Pilosof  _et al._ (2017) Shai Pilosof, Mason A. Porter, Mercedes Pascual, and Sonia Kéfi, “The multilayer nature of ecological networks,” Nature Ecology & Evolution 1, 0101 (2017).
  * Iribarren and Moro (2009) José Luis Iribarren and Esteban Moro, “Impact of human activity patterns on the dynamics of information diffusion,” Phys. Rev. Lett. 103, 038702 (2009).
  * Mucha  _et al._ (2010) P. J. Mucha, T. Richardson, K. Macon, M. A. Porter, and J.-P. Onnela, “Community structure in time-dependent, multiscale, and multiplex networks.” Science 328, 876–878 (2010).
  * Gauvin  _et al._ (2014) Laetitia Gauvin, André Panisson, and Ciro Cattuto, “Detecting the Community Structure and Activity Patterns of Temporal Ne tworks: A Non-Negative Tensor Factorization Approach,” PLoS ONE 9, e86028 (2014).
  * Peixoto (2015) Tiago P. Peixoto, “Inferring the mesoscale structure of layered, edge-valued, and time-varying networks,” Phys. Rev. E 92, 042807 (2015).
  * Ghasemian  _et al._ (2016) Amir Ghasemian, Pan Zhang, Aaron Clauset, Cristophe r Moore, and Leto Peel, “Detectability Thresholds and Optimal Algorithms for Community Structur e in Dynamic Networks,” Phys. Rev. X 6, 031005 (2016).
  * Sapienza  _et al._ (2017) Anna Sapienza, Alain Barrat, Ciro Cattuto, and Laetitia Gauvin, “Estimating the outcome of spreading processes on networks with incomplete information: a mesoscale approach,” (2017), arXiv:1709.01806.
  * Li  _et al._ (2017) A. Li, S. P. Cornelius, Y. Y. Liu, L. Wang, and A. L. Barabási, “The fundamental advantages of temporal networks,” Science 358, 1042–1046 (2017).
  * Guimerà and Sales-Pardo (2009) R. Guimerà and M. Sales-Pardo, “Missing and spurious interactions and the reconstruction of complex networks.” Proc. Natl. Acad. Sci. U. S. A. 106, 22073–22078 (2009).
  * Peixoto (2014) T. P. Peixoto, “Hierarchical block structures and high-resolution model selection in large networks,” Phys. Rev. X 4, 011047 (2014).
  * (14) T Vallès-Català, T.P. Peixoto, R. Guimerà, and M. Sales-Pardo, “On the consistency between model selection and link prediction in networks,” ArXiv:1705.07967 [stat.ML].
  * Guimerà and Sales-Pardo (2013) R. Guimerà and M. Sales-Pardo, “A network inference method for large-scale unsupervised identification of novel drug-drug interactions,” PLoS Comput. Biol. 9, e1003374 (2013).
  * Rovira-Asenjo  _et al._ (2013) N. Rovira-Asenjo, T. Gumí, M. Sales-Pardo, and R. Guimerà, “Predicting future conflict between team-members with parameter-free models of social networks.” Sci. Rep. 3, 1999 (2013).
  * Godoy-Lorite  _et al._ (2016a) Antonia Godoy-Lorite, Roger Guimerà, Cristopher Moore, and Marta Sales-Pardo, “Accurate and scalable social recommendation using mixed-membership stochastic block models,” Proc. Natl. Acad. Sci. U.S.A. 113, 14207 –– 14212 (2016a).
  * Airoldi  _et al._ (2008) Edoardo M Airoldi, David M Blei, Stephen E Fienberg, and Eric P Xing, “Mixed membership stochastic blockmodels,” J. Mach. Learn. Res. 9, 1981–2014 (2008).
  * Godoy-Lorite  _et al._ (2016b) A Godoy-Lorite, R Guimerà, and M Sales-Pardo, “Long-term evolution of email networks: Statistical regularities, predictability and stability of social behaviors,” PLoS ONE 11, e0146113 (2016b).
  * Bacco  _et al._ (2017) Caterina De Bacco, Eleanor A. Power, Daniel B. Larremore, and Cristopher Moore, “Community detection, link prediction, and layer interdependence in multilayer networks,” CoRR abs/1701.01369 (2017).
  * Stanley  _et al._ (2016) Natalie Stanley, Saray Shai, Dane Taylor, and Peter J. Mucha, “Clustering Network Layers with the Strata Multilayer Stochastic Block Model,” IEEE Transactions on Network Science and Engineering 3, 95–105 (2016), 1507.01826 .
  * Jefries (2017) Bill Jefries, “MMSBM for sparks,” https://github.com/billjeffries/mixMemRec (2017).
  * Menden  _et al._ (2018) Michael Patrick Menden, Dennis Wang, Yuanfang Guan, Michael Mason, BenceSzalai, Krishna C Bulusu, Thomas Yu, Jaewoo Kang, Minji Jeon, Russ Wolfinger, Tin Nguyen, Mikhail Zaslavskiy, AstraZeneca-Sanger Drug Combination DREAM Consorti, In Sock Jang, Zara Ghazoui, Mehmet Eren Ahsen, Robert Vogel, EliasChaibub Neto, Thea Norman, Eric KY Tang, Mathew J Garnett, Giovanni Di Veroli, Stephen Fawell, Gustavo Stolovitzky, Justin Guinney, Jonathan R Dry, and Julio Saez-Rodriguez, “Community assessment of cancer drug combination screens identifies strategies for synergy prediction,” bioRxiv , 200451, https://doi.org/10.1101/200451 (2018).
  * Note (1) This particular choice of threshold leads, when models are properly calibrated in a frequentist sense, to both precision and recall having very similar values (see Supplementary Information).
  * Gneiting  _et al._ (2007) Tilmann Gneiting, Fadoua Balabdaoui, and Adrian E. Raftery, “Probabilistic forecasts, calibration and sharpness,” J. R. Statist. Soc. B 69, 243–268 (2007).
  * Note (2) Note that our results do not depend on the number of latent dimensions allowed for the temporal layers L, since for $L>2$ we also find that temporal layers have $\eta_{\ell\gamma}\not=0$ only for two latent groups $\gamma$.
  * Note (3) Due to the sparsity of observations in these networks many interactions were never observed in the training sets, and thus no group memberships could be assigned to the links ($e_{ij}$) corresponding to those interactions. We solved this cold start problem by, at each iteration, assigning them the average membership of the observed interactions $\zeta_{e}={\left\delimiter 0"426830A\zeta_{f}\right\delimiter 0"526930B}_{f\in G^{O}}$. Analogously, if a node $i$ had no observed interactions in the training set, at each iteration we set its membership vector as the average of membership vectors for nodes with observed interactions $\theta_{i}={\left\delimiter 0"426830A\theta_{k}\right\delimiter 0"526930B}_{k\in G^{O}}$.
  * Newman and Clauset (2016) M. E. J. Newman and Aaron Clauset, “Structure and inference in annotated networks,” Nat. Comm. 7, 11863 (2016).
  * Hric  _et al._ (2016) Darko Hric, Tiago P. Peixoto, and Santo Fortunato, “Network structure, metadata, and the prediction of missing nodes and annotations,” Phys. Rev. X 6, 031038 (2016).
  * Cobo  _et al._ (2017) S. Cobo, A. Godoy-Lorite, M. Sales-Pardo, and R. Guimerà, “Optimal prediction of decisions and model selection in social dilemmas using block models with metadata,” (2017).


# On the Impact of Fixed Point Hardware for Optical Fiber Nonlinearity
Compensation Algorithms

Tom Sherborne*, Benjamin Banks*, Daniel Semrau, , Robert I. Killey , Polina
Bayvel, , and Domaniç Lavery Manuscript received January, 2018. This work was
funded by United Kingdom (UK) Engineering and Physical Sciences Research
Council (EPSRC) Programme Grant UNLOC (UNLocking the capacity of Optical
Communications), EP/J017582/1. D. Lavery is supported by the Royal Academy of
Engineering under the Research Fellowships scheme.*These authors contributed
equally to this work. T. Sherborne and B. Banks were students within the
Department of Electronic and Electrical Engineering, University College
London, London WC1E 7JE, U.K. (e-mail<EMAIL_ADDRESS>zceebjb@ucl.ac.uk).D. Lavery, D. Semrau and P. Bayvel are with the Optical
Networks Group, Department of Electronic and Electrical Engineering,
University College London, London WC1E 7JE, U.K. (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>p.bayvel@ucl.ac.uk).

###### Abstract

Nonlinearity mitigation using digital signal processing has been shown to
increase the achievable data rates of optical fiber transmission links. One
especially effective technique is digital back propagation (DBP), an algorithm
capable of simultaneously compensating for linear and nonlinear channel
distortions. The most significant barrier to implementing this technique,
however, is its high computational complexity. In recent years, there have
been several proposed alternatives to DBP with reduced computational
complexity, although such techniques have not demonstrated performance
benefits commensurate with the complexity of implementation. In order to fully
characterize the computational requirements of DBP, there is a need to model
the algorithm behavior when constrained to the logic used in a digital
coherent receiver. Such a model allows for the analysis of any signal recovery
algorithm in terms of true hardware complexity which, crucially, includes the
bit-depth of the multiplication operation. With a limited bit depth, there is
quantization noise, introduced with each arithmetic operation, and it can no
longer be assumed that the conventional DBP algorithm will outperform its low
complexity alternatives. In this work, DBP and a single nonlinear step DBP
implementation, the Enhanced Split Step Fourier method (ESSFM), were compared
with linear equalization using a generic software model of fixed point
hardware. The requirements of bit depth and fast Fourier transform (FFT) size
are discussed to examine the optimal operating regimes for these two schemes
of digital nonlinearity compensation. For a 1000 km transmission system, it
was found that (assuming an optimized FFT size), in terms of SNR, the ESSFM
algorithm outperformed the conventional DBP for all hardware resolutions up to
13 bits.

###### Index Terms:

Optical Fiber Communication, Digital Signal Processing, Nonlinearity
Compensation

##  I Introduction

Signal processing techniques to overcome optical fiber nonlinearity have
evolved considerably in recent years. With the aim of increasing the capacity
and reach of coherent detection systems, the use of digital signal processing
(DSP) algorithms has been shown to significantly improve system performance
[1] and reduce complexity [2] when compared to earlier systems [3]. The DBP
algorithm is conventionally applied at the receiver, however it can also be
applied at the transmitter (i.e., digital pre-compensation [4]), or as some
combination of both link endpoints, as investigated in [5]. In all cases, the
aim of DBP is to simulate a reverse ‘virtual’ link (see Fig. 1) and,
consequently, reverse the deterministic nonlinear and linear effects of fiber
propagation[3].

Figure 1: Transmission model used to investigate the performance of various
DBP schemes. This is implemented as a series of linear and nonlinear operators
applied to the buffered signal vector. Ideal DBP describes the most intensive
implementation of DBP with $k$ steps equal to the number of steps applied over
$N_{s}$ spans in simulated forward propagation.

However, the promising performance gains achieved with DBP typically incur a
computational expense beyond the capabilities of present digital signal
processing hardware[6], due to the use of the split step Fourier method to
solve the reverse Nonlinear Schrödinger equation (NLSE) as a Manakov
system[7]. In order to achieve real-time processing, several implementations
of DBP have been proposed with differing perspectives on the trade-off between
performance and complexity. These included increasing the nonlinear step size
of the split step Fourier method in the DBP virtual link to minimize the
number of calculations required for signal recovery. 1 step-per-span DBP
(1sps-DBP) [3] has been proposed as a natural trade-off between performance
and complexity for single channel DBP. More advanced proposals include the low
pass filter DBP (LPF-DBP) - an implementation of DBP that introduces a phase
noise filter in the nonlinear step of the algorithm[8]. High frequency phase
noise components are present in the phase spectrum when measured across
multiple spans, and the use of LPF-DBP acts to suppress such components to
achieve reasonable performance with $<1$ step per span. LPF-DBP has been shown
to perform comparably to 1sps-DBP with 1 step every 4 spans [8], corresponding
to an approximately 75% reduction in the number of calculations required for
NLC. Note that the overall numerical complexity of the aforementioned NLC
schemes is dependent on the link length.

DBP is typically implemented as a Wiener-Hammerstein system[6] on a ‘virtual’
link expressing a reverse journey to the physical forward propagation. The
conventional Manakov system for modeling propagation concatenates operators of
fixed distance to describe the full link, and so the complexity of standard
DBP, which uses this approach, is proportional to the link length. A recent
development in the field has been the design of NLC methods with complexity
that is decoupled from the physical link length[9].

For realizable NLC in a coherent receiver, any solution would ideally
demonstrate performance improvement over electronic dispersion compensation
(EDC) while avoiding a significant increase in computational complexity. This
has spurred interest in single step-per-link implementations of DBP, which
reduce the computational complexity compared to multiple step (e.g 1 step-per-
span) schemes. The Enhanced Split Step Fourier Method (ESSFM) is a general
implementation of filtered DBP, which applies the reversed channel with the
step size equal to the link length [9]. Similarly to LPF-DBP, the ESSFM uses
an additional operation stage in the nonlinear algorithm sub-step. The ESSFM
algorithm has demonstrated a performance increase of 0.7 dB [10] over a linear
Feed-Forward Equalizer (FFE) system with a comparable demand of complex
multiplications[10]. Both methods of augmented DBP use nonlinear optimization
routines for the determination of a nonlinear filter[11, 10]. This routine
solves a numerical problem and not the underlying problem of solving an
indeterminate equation describing a physical system, and therefore it is
uncertain if either method demonstrates optimality in their solutions.

These implementations of DBP have demonstrated the potential performance
benefits of using NLC over EDC in offline implementations. However, to make a
case for the practical deployment of such systems there remains a question of
DBP performance when constrained to the fixed point arithmetic required in a
high throughput coherent receiver.

This has previously been investigated in [12], however the study was
constrained by the use of time domain dispersion filters, with a manually
optimized dispersion compensation ratio between the linear filters, and a
specific target ASIC implementation. To extend the general understanding of
quantization effects in NLC methods, here we investigate frequency domain
dispersion compensation, and make no assumptions on the target hardware. A
further advantage of this hardware- and system-agnostic approach is that it
facilitates a clearer comparison between different NLC methods, where a change
in dispersive block length can have a significant impact on performance, while
removing the consideration of distortions introduced by finite-length time
domain filters[13].

In this paper, we develop a model of finite precision arithmetic and simulate
the operation of several implementations of DBP. We assess and compare the
performance of constrained DBP using the linear EDC algorithm as a benchmark.
Finally we determine operating regimes for which a specific implementation of
NLC outperforms EDC in terms of SNR.

This paper is organized as follows. In Section II-A, different NLC algorithms
are introduced. In Section II-B, the model of fixed-point arithmetic and
hardware simulation is described in detail. The results of the transmission
simulations are given in Section III-A for the DBP fixed point algorithm and
Section III-B for the ESSFM fixed point algorithm. Conclusions are drawn in
Section IV.

##  II Methodology

###  II-A Algorithm Design

In modeling the behavior of digital NLC we used a simulation model of an
optical link (Fig. 1), assuming ideal noise-free transceivers, with the link
parameters given in Table I. The signals under test were single channel 32 GBd
dual polarization (DP) quadrature phase shift keying (QPSK) and 16-ary
quadrature amplitude modulation (DP-16QAM) signals. The signal was sampled at
4 samples/symbol and shaped using a root-raise cosine (RRC) filter. In this
approach, we numerically solved the Manakov system for forward pulse
propagation to simulate 1000 km (25$\times$40 km span) of standard single-mode
fiber (SSMF), as a representative long-haul transmission system, and following
the methodology in [10]. An EDFA was included after each span to fully
compensate for the power loss due to signal attenuation over the span.

TABLE I: Summary of System Parameters Parameter | Value | Units
---|---|---
Fiber attenuation | 0.2 | dB/km
Dispersion parameter | 17 | ps/(nm $\cdot$ km)
Fiber nonlinear coefficient | 1.2 | 1/(W $\cdot$ km)
Span length | 40 | km
Simulation step size | 100 | m
1sps-DBP step size | 40,000 | m
1spl-DBP/ESSFM step size | Varies | m
Symbol rate | 32 | GBd
EDFA noise figure | 5 | dB
Pulse shape | RRC, 1% | rolloff
DBP Wiener-Hammerstein split | 0.85 | -
ESSFM Wiener-Hammerstein split | 0.4 | -

Fig. 2 shows the receiver model used in this work. The received data is
detected, resampled to 2 samples/symbol and normalized to unit average power.
The signal is then processed using either EDC (implemented as a single step
frequency domain filter), or the NLC algorithm under test. The approaches to
NLC in this work are the ESSFM[9], 1 step-per-span DBP, and 1 step-per-link
DBP. In the implementation of each NLC algorithm, fixed-point arithmetic is
used to model the finite precision available in practical hardware. Details of
the fixed point modules are discussed in Section II-B. For both DBP and EDC,
the Fast Fourier Transform (FFT) is based on the radix-2 Cooley-Tukey
algorithm [14] using an ‘overlap and save’ procedure to model a buffered
application of dispersion compensation. Exponentials were approximated using
the coordinate rotation digital computer (CORDIC) algorithm to reduce the
computational demand of approximation using a long Taylor series. The filter
coefficients for the ESSFM algorithm were computed offline using a double
floating-point precision non-linear optimization algorithm, but quantized to
the appropriate bit depth at runtime. Following NLC or EDC, matched filtering
was applied to the signal, before normalization and downsampling to 1
sample/symbol. SNR estimation was performed on the central $2^{14}$ symbols
according to the method detailed in [15] and is reported for the optimum
launch power unless otherwise specified.

Figure 2: The signal processing chain used in this work. The complex valued
pulse train is captured using a model of balanced detection, downsampled to 2
samples/symbol, NLC applied and match filtered before unitary scaling and an
SNR for the symbols estimated. This routine is also the function, $f(x)$, for
the optimisation of the ESSFM filter coefficients.

DBP was applied using the split-step Fourier method (SSFM) to solve the
Manakov equation [16] with the inverse loss, dispersion and nonlinearity fiber
coefficients according to the parameters in Table I. The nonlinear interaction
between accumulated ASE noise and the co-propagating signal limits the channel
performance even with application of ideal DBP. This can be attributed to
nonlinear phase noise accumulation (Gordon-Mollonauer effect)[6]. The single
channel transmission at 32 GBd with 40 km spans was chosen for equitable
comparison of algorithm performance between DBP and ESSFM, again, following
the methodology of [10]. The Manakov inverse channel model[16], over some
incremental step, $h$, approximates the deterministic effects of the channel
as

$\begin{split}\mathbf{E}(z+h,T)&\simeq\textrm{exp}(h\hat{D})\textrm{exp}(h\hat{N})\mathbf{E}(z,T)\\\
&\simeq\textrm{exp}\left(\frac{h}{2}\hat{D}_{1}\right)\textrm{exp}(L_{\textrm{eff}}\hat{N})\textrm{exp}\left(\frac{h}{2}\hat{D}_{2}\right)\mathbf{E}(z,T),\end{split}$
(1)

where $\hat{D}$ is the linear operator, given by

$\hat{D}=\frac{i\beta_{2}}{2}\frac{\partial^{2}}{\partial t^{2}},$ (2)

and $\hat{N}$ is the nonlinear operator, defined as

$\hat{N}=-i\gamma\frac{8}{9}\mathbf{E^{H}E}-\frac{\alpha}{2}.$ (3)

Here, $\mathbf{E}=[E_{X},E_{Y}]$ is the optical field in the $X$ and $Y$
polarization states, $\alpha$ is the fiber loss parameter, $\beta_{2}$ is the
group velocity dispersion, $\gamma$ is the fiber nonlinearity coefficient, and
$L_{\textrm{eff}}$ is the effective length of $h$.

Following the approach described in [6], a 3- block Wiener-Hammerstein model
for each step was used for the DBP. Linear compensation sub-steps were
performed in the frequency domain with a circular convolution using 4 FFTs per
step. For DBP, adjacent linear blocks were combined over the total link for
$N_{s}+1$ total FFT blocks. Nonlinear substeps were applied in the time domain
as a phase shift proportional to the instantaneous power of the back
propagating signal. The responsibility of application of the linear components
between each $\hat{D}_{i}$ block in (1) is dictated by the Wiener-Hammerstein
(WH) split. This parameter has the effect of controlling the position of the
nonlinear $\hat{N}$ term across step distance, $h$.

To implement the ESSFM algorithm proposed in [9], we use a modified nonlinear
operator

$\hat{N}_{k}=-i\gamma
L_{eff}\left(\sum_{i=0}^{N_{c}}c_{i}\left(|\mathbf{E_{k-i}}|^{2}+|\mathbf{E_{k+i}}|^{2}\right)\right)$
(4)

for each sample, $k$, with the application of the ESSFM filter coefficients,
$c_{i}$, for filter size $N_{c}+1$. The filter coefficients are determined
through a nonlinear multivariate equation solver. The DSP chain for signal
processing, including the ESSFM, is formulated as some function, $f(x)$, which
returns a final SNR. In this work, a Quasi-Newton cubic line search procedure
was used to minimize a negated SNR output of $f(x)$. This blind optimization
approach uses gradient descent with a relative tolerance of $10^{-3}$ to
converge on a set of filter coefficients. Similarly to [10], we applied
nonlinear optimization to maximize the SNR of a sample sequence using the
coefficients, $c_{i}$ as the degrees of freedom. The number of coefficients in
the filter was restricted to a power of two, again for direct comparison with
[10]. The number of multiplications per transmitted symbol for the ESSFM
$\hat{N}$ operator can be approximated by $N_{s}(2N_{c}+1)$ when the
multiplications for instantaneous power are disregarded111Note that the ESSFM
filter is a symmetric FIR filter and, thus, this filter could be efficiently
implemented using frequency domain convolution. Given the relatively short
filter length, the computational complexity in this case would be comparable,
but this observation is noted here for completeness.. This linear relationship
suggests the possibility for an optimality trade-off between ESSFM filter size
and performance in SNR. Diminishing returns for the system in [10] were
observed for filters of 32 tap weights and above. Here, this investigation was
repeated for our system parameters, as it can be seen that additional
multiplications in the $\hat{N}$ operator have implications for the
performance limitations from the quantization overhead in the fixed point
model; there is a trade off between nonlinear performance gain and
quantization noise introduced from the additional filter taps. The theoretical
basis for the operation of the ESSFM and the filter coefficients is further
discussed in [17, 10].

As the ESSFM operates at the 1 step-per-link level, the Wiener-Hammerstein
(WH) split now determines the position of the only nonlinearity compensation
operator across all spans. The optimal WH split value, previously determined
to be 0.85 for conventional DBP, assumes a step, $h$, of maximum size of a
single span[6]. We found this value for the split to be sub optimal in our
system for DBP, however the SNR benefits from additional optimization only
reached a maximum of $0.01$ dB and, therefore, a WH-split value of 0.85 was
retained when simulating DBP as the NLC under test. However for the ESSFM, the
WH split was fully re-optimized. The intuition for this is that since the
distribution of power over distance for a step size of 1 span length is only a
subset of the distribution over the full link, the location of the single
$\hat{N}$ would likely change. In [10] the ESSFM is proposed using the Wiener
Model of the $\hat{N}$ term preceding a singular $\hat{D}$ term. We found this
arrangement to be sub-optimal in terms of nonlinear compensation performance,
despite the additional quantization from two FFT pairs per link, and we
subsequently used an optimized WH split value of 0.4.

###  II-B Hardware Model

The primary aim of this work is to quantify the performance of different NLC
algorithms when constrained to fixed-point (FXP) arithmetic. In digital
circuits, signals are stored using a bit-wise representation. For $b$ bits it
is possible to represent any integer value in the range
$\left[0,2^{b}-1\right]$. The precision of computations improves with bit
depth, however high precision comes at a cost of more associated digital
circuitry, and therefore complexity, which introduces greater power
consumption, circuit area and, ultimately, cost. At the same time there is a
need, especially in the context of a receiver, to maintain enough numerical
precision so that the received signal meets a given SNR requirement. This
trade-off is investigated in this work where solution precision will have a
direct influence the signal-to-noise ratio (SNR) of a signal which has gone
through the arithmetic process of NLC. For the different NLC schemes, it is
possible to compare the performance and complexity merits of such algorithms
in a realistic scenario. To do this, it is necessary to model the bit-level
arithmetic occurring in realtime receiver hardware.

A software model of FXP arithmetic for an arbitrary quantization level was
developed. This model remains agnostic to the hardware platform, with no
intended target system. The basis for this model are functional blocks for FXP
addition and multiplication, which receive some quantization level as an input
parameter and return the arithmetic result. The constraints investigated here
are limited to bit depth (the number of logical bits permitted in FXP
arithmetic) as a primary factor in the degree of accumulated quantization and
the size of the FFT window. Intuitively, a decrease in bit depth leads to
greater quantization noise as the bit stream is increasingly compressed. This
limits the ability of the FXP model to approximate the best case (64-bit
floating point) simulation model, resulting in some level of performance
degradation. Arithmetic precision is, therefore, a key design parameter which
has a strong influence on the performance of any NLC scheme.

These atomic, arithmetic functions were used to emulate the operation of more
complex operations, including an FFT. A fixed bit depth was assumed throughout
the system and bit depth expansion was not permitted across multiple
functions, as is typical in real-time signal processing hardware. However bit
depth expansion is permitted for intermediate operations only222We refer to an
intermediate operation as one inside a single function block but between the
atomic function units. The primary usage of this intermediate bit expansion is
during the FFT block, wherein conditional scaling of bit depth is employed to
minimize the addition of quantization noise from successive FXP operations.
This problem could also have been approached by choosing a different
quantization interval, however in this work the bounds $x\in[-1,1)$ are
strictly enforced.. Each number is represented in simulation using the
fractional two’s complement fixed point format, with each representable value
in the range $-1\leq x<1$. This format was chosen as there is a guarantee that
the product of every multiplication will have some magnitude equal or less
than that of either input. The sampling of a continuous signal to some
discrete amplitude introduces quantization noise, resulting in error
propagation as a result of the finite precision of the system. As a sample
propagates through a cascade of multiplication stages (e.g., in the FFT) then
each multiplication introduces some quantization noise. The signal-to-
quantization noise ratio of a signal quantized to $B$ bits [18] can be
described as

$SQNR=\frac{1}{\sigma^{2}},$ (5)

where the quantization noise variance is given by

$\sigma^{2}=\frac{2^{-2B}}{12}$ (6)

The logarithm of Eq. (5) reveals a gain of 6 dB in SQNR for the provision of
an extra bit. The crux of Eq. (5) in the context of digital NLC is the trade-
off between maximizing the number of DBP steps in order to ensure solution
accuracy and minimizing the number of steps to prevent the accumulation of
quantization noise.

When an input signal is first quantized to a set bit depth, each amplitude is
rounded towards the nearest available discrete amplitude. In the event of an
arithmetic overflow333Arithmetic overflow refers to some discrete number
outside of the $[-1,1)$ bounds. in this system, the erroneous value is clipped
at maximum or minimum quantized amplitude. This method also ensures the
quantization noise distribution remains zero-centered.

As previously noted, linear filtering stages of NLC are performed in the
frequency domain. Transformations between the time and frequency domain are
performed by a radix-2 Cooley-Tukey FFT [14]. This format of FFT is chosen for
the ability to implement a wide range of FFT sizes and because the complexity
of this FFT is well studied[10]. We did not consider split radix FFTs in this
work, which may have changed the overall NLC performance. However, because the
same FFT implementation is common to all algorithms considered, any
performance change would affect all NLC implementations. The quantization
noise introduced from one butterfly operation in an FFT is $4\sigma^{2}$ with
a quantization noise variance introduced from an FFT as
$4(N-1)\sigma^{2}$[18]. This corresponds to the optimal SNR achievable for a
signal processed using an FFT and, as such, places a ceiling on total system
performance.

For an accurate representation of online processing of a continuous stream of
input samples, the incoming signal is buffered and processed in batches equal
to the size of the FFT under test. The Overlap and Save algorithm was used to
be able to model the continuous filtering required in hardware. The size of
the buffer to process must correspond to a radix-2 FFT size, or $2^{n}$, where
$n\in{}\\{5,\ldots 15\\}$. An overlap size of $N/4$ was chosen, where $N$ is
the length of a signal buffer, in correspondence with the work carried out in
[10]. The Overlap and Save algorithm used here also applies the dispersion map
as the linear substeps of DBP and the full CDC operation. Nonlinearity
compensation is provided in the time domain by multiplying the signal with a
vector of complex exponentials. Using Euler’s identity, it is possible to
approximate the complex exponential using the CORDIC algorithm[19]. In this
algorithm, the sine and cosine of a provided angle is approximated
iteratively, with the approximate solution accuracy improving with each
iteration. Note that, in contrast to the Taylor expansion model used in [12],
CORDIC provides a greater degree of accuracy with an acceptable computational
complexity.

##  III Results

###  III-A Digital Back Propagation

The initial simulations of the transmission of $2^{18}$ DP-QPSK modulated
samples at 32 GBd over a 1000 km fiber link. DBP was implemented with the
parameters described in Section II-A with the finite precision logic described
in Section II-B. For comparison, we also implemented CDC using the same FXP
modeling constraints. The FFT size was optimized in each case and is therefore
not included as a parameter. In the following, we examine the behavior of DBP
in the FXP environment, and discuss the additional complexity of this form of
NLC compared to CDC.

Fig. 3 shows the SNR improvement relative to EDC at the same FXP bit depth,
across a range of numbers of nonlinear steps per link. For an increased number
of steps in DBP, the virtual link model improves in efficacy, as expected,
however this is at the cost of an increased FXP bit depth, which is required
to reduce the impact of quantization noise when multiple steps are used.
Additionally, Fig. 3 identifies that there exists some quantization limit in
performance, at a given number of steps-per-link, which can be reached with a
sufficient bit depth in the digital logic. In the single step per link case,
with the fewest multiplications, there is a negligible SNR gain of 0.1 dB over
CDC for bit depths greater than 13 bits. We observe no gain in SNR over CDC
for bit depths below this, for any number of steps-per-link. The greatest
performance difference over CDC is a 3.5 dB gain, observed using 50 steps per
link, equivalent to 2 steps per span, was at a bit depth of 15 bits. However,
when compared to CDC it can be seen that any gain requires a minimum of 13
bits.

$8$$10$$12$$14$$16$$-4$$-2$$0$$2$$4$Bit depth$\Delta$SNR over CDC (dB)1 SPL10
SPL25 SPL50 SPL Figure 3:  SNR improvement for DBP compared to CDC at 1000 km
for DP-QPSK transmission. Both DBP and CDC are simulated with FXP arithmetic
and optimized FFT size. DBP is simulated with a range of steps-per-link (SPL)
to investigate the influence of the number of steps on the quantization error
on the signal.

We expand on this analysis in Fig. 4 by examining the absolute SNR performance
of DBP against the variation in steps-per-link between 7 and 16 bits. From
this result it can be inferred that for $\leq 12$ bits, there is insufficient
arithmetic precision to tolerate the increased accumulation of quantization
noise arising from increasing the number of DBP steps per link. While
increasing the steps per link increases the accuracy of theoretical DBP our
results show that increasing the number of FFTs and multiplications is
monotonically detrimental for low bit depths due to quantization error.
However for $\geq 13$ bits, there is a noted transition to an operating regime
for which an increase in the steps per link results in some SNR gain. This
confirms that at 1000 km, a minimum of 13 bits are required to overcome the
quantization overhead present in our implementation of DBP. In particular for
13 bits, there is an increase in SNR as a result of increased steps up to 25
SPL, where a peak SNR of $23.6$ dB is observed. Beyond this a subsequent
decrease in SNR is observed, likely as the quantization at 50 SPL becomes
dominant even for 13 bits resolution. Another way to view this result is that
there exists an optimum number of steps per link for a given FXP hardware bit
depth. Previously reported results in this area show that DBP performance
increases monotonically with number of steps per link and, thus, this result
highlights the importance of considering fixed point arithmetic in performance
analyses.

$1$$10$$-5$$0$$5$$10$$15$$20$$25$Steps per linkSNR (dB)7 bits8 bits9 bits10
bits11 bits12 bits13 bits14 bits15 bits16 bits Figure 4:  SNR against the
steps per link for DBP at 1000 km for 32 GBd DP-QPSK transmission. DBP is
simulated with FXP arithmetic with the number of steps in the virtual fiber
model and bit depth varied.

###  III-B Enhanced Split Step Fourier Method

$-4$$-3$$-2$$-1$$0$$20.5$$21$$21.5$$22$$22.5$$23$Launch Power, $P_{lch}$
(dBm)SNR (dB)CDCESSFM $N_{c}$=16ESSFM $N_{c}$=32ESSFM $N_{c}$=64ESSFM
$N_{c}$=128ESSFM $N_{c}$=256 Figure 5: SNR against launch power for the ESSFM
applied as NLC to a 32 GBd DP-QPSK signal. The CDC baseline is also shown.
This simulation was produced using double-precision floating point arithmetic,
representing the ceiling of achievable performance for a comparable fixed
point simulation which will have a non-negligible quantization penalty.

The ESSFM algorithm described in Section II-A was simulated in the system
detailed by Fig. 2. For the DP-QPSK scenario, Fig. 5 identifies the variation
in SNR against launch power. Each power and coefficient filter size was
separately optimized up to a filter of size $N_{c}=256$.

Improvements in SNR compared to CDC were observed for all input powers near
the optimum. The maximum SNR improvement was 0.7 dB for DP-QPSK and 0.8 dB for
DP-16QAM with a 256-tap filter. It should be noted that the optimum launch
power increased to -2.5 dBm for both the DP-QPSK and DP-16QAM cases. These two
results together give confidence in the effective operation of the ESSFM, and
are comparable to the original results for this algorithm as presented in
[10].

Similarly to previous work, marginal gains in SNR were observed for
significant increases in the filter size. Fig. 6 highlights the aforementioned
trade-off between the filter size and the efficacy of the NLC operator. This
result indicates that the algorithm saturates for the given system parameters
near $N_{c}=128$, demonstrating only 0.1 dB increase in SNR for a doubling of
the filter size and increased hardware complexity. For the remainder of this
work, the ESSFM reported in this work uses a filter of size $128$ because of
this observed saturation.

$50$$100$$150$$200$$250$$0.7$$0.72$$0.74$$0.76$$0.78$$0.8$Number of ESSFM tap
weights $N_{c}$$\Delta$SNR (dB)QPSK16QAM Figure 6: $\Delta$SNR over CDC for
the ESSFM at optimum $P_{lch}$ for increasing filter size.

Previous research into modifying the nonlinear operator in DBP[10, 20, 11]
used a black box approach generating an optimal filter. We follow this
methodology, however we chose to inspect further the behavior of the filter
returned through the optimization routine detailed in Section II-A. In Fig. 7,
the 128-tap filter for processing a DP-QPSK signal at the optimum launch power
is analyzed in the frequency domain. The profile of this filter is effectively
a weighting of a neighboring samples’ power to each sample. Fig. 7
demonstrates that this filter exhibits a low pass filtering behavior that was
not a constraint intended as part of the design. We note that this highlights
more similarity between the methods of the LPF-DBP and the ESSFM, as the
coefficients in the ESSFM were not originally proposed as a filter of any
specific kind. With the understanding that the $\check{N}$ operator in the
ESSFM is performing a linear phase444The linear phase response of the ESSFM
128-tap filter is not shown in this work for brevity. Also not shown is
similar low pass behavior for the 128-tap filter designed for DP-16QAM. low
pass filtering with one NLC term per link. We suggest that in future work, the
approaches of [20] and [10] could be integrated, e.g., through direct
optimization of the filter bandwidth within the ESSFM.

$0.0$$0.25$$0.50$$0.75$$1$$-20$$0$$20$Normalised frequency
($\times\pi$/sample)Magnitude (dB) Figure 7: Frequency response of the 128-tap
ESSFM $\check{N}$ operator for 32 GBd DP-QPSK transmission.

Finally, the behavior of the ESSFM as NLC in a finite precision arithmetic
simulation is explored. Fig. 8 shows the maximum SNR improvement over FXP CDC
with a varying number of bits and FFT size optimized for each bit depth. It
can be seen that, as anticipated, there exists some quantization penalty
limiting the finite precision algorithm performance when compared to double
precision. The difference between the overheads of 0.024 dB for DP-QPSK and
0.18 dB for DP-16QAM highlights a potential level of increased sensitivity to
quantization for multi-level modulation formats. As the level of quantization
noise is constant for any distance with the ESSFM, there exists some link
length at which channel noise dominates over quantization noise and this
overhead becomes negligible. The exact location of this crossover point was
not explored, as it is an inherently system-dependent parameter, however it is
estimated from preliminary results to be near a link length of 2000 km for the
system considered herein.

$6$$8$$10$$12$$14$$16$$-3$$-2$$-1$$0$$1$Number of bits$\Delta$SNR
(dB)QPSK16QAM Figure 8: $\Delta$SNR over CDC for the ESSFM against bit depth
at 1000 km for 32 GBd DP-QPSK and DP-16QAM transmission.

The results in Fig. 8 also highlight the limitations of even low complexity
DBP when compared to CDC, as a minimum of 10 bits are required to improve over
linear equalization. To achieve gains above 0.7 dB over CDC in the DP-QPSK
case, at least 13 bits are required, with 14 bits required for the DP-16QAM
algorithm to exhibit a benefit of 0.6 dB in SNR. These requirements for bit
depth are greater than the operating size of typical ASICs used in high speed
coherent receivers, and present a barrier to achieving the full SNR gains
promised by this NLC algorithm. Without further investigation into nonlinear
signal recovery, our work appears to show the continued strength of linear
algorithms for hardware implementation.

Comparing Fig. 3 and Fig. 8, we can observe that, at a 1 step-per-link level,
the low complexity ESSFM demonstrates improvements in SNR performance when
compared to both DBP and EDC. However, this gain in SNR remains less than $1$
dB and is sub-optimal when compared to a DBP using a greater number of steps-
per-link. Provided there is a sufficient number of bits to saturate the system
performance, Fig. 3 demonstrates that more steps-per-link will yield greater
SNR improvement compared to the addition of the nonlinear phase filter in the
ESSFM. Nevertheless, it should be noted that the ESSFM outperforms DBP for all
bit depths $\leq$13 bits, irrespective of the number of steps used for DBP.

##  IV Conclusion

The behavior of algorithms for optical fiber nonlinearity compensation were
investigated in detail under the constraint of fixed-point arithmetic. Using a
simulation including an ideal transmitter and coherent receiver, limitations
on performance for signal recovery in DSP, have been shown, and are primarily
a result of quantization noise. It was found that a minimum of 13 bits is
required to observe any additional SNR gain from $\geq{}1$ steps-per-link DBP.
A similar bit depth is required to reach the maximum achievable performance
benefit over CDC for the ESSFM algorithm, however we note that some
improvement over CDC can be seen with the reduced constraint of 10 bits.
Although other algorithms for simplified DBP exist [21, 22], the work
presented herein has shown that it is crucial to include an analysis of bit
depth in order to accurately compare the performance of different NLC
algorithms, and that the number of mathematical operations, alone, cannot
serve as a proxy for complexity.

We also highlight, for the first time, the similarities in behavior of the
ESSFM and LPF-DBP through the observation that the frequency response of the
ESSFM phase filter exhibits low pass filtering behavior. Future work on this
topic could include the investigation into the behavior of the aforementioned
NLC in a multi-channel or super-channel system, as well as expanding on the
impact of modulation format on the quantization penalty of the finite-
precision ESSFM.

More generally, we note that the fixed point models developed herein could be
used as a basis for estimating the true performance of any NLC algorithm, and
should, therefore, be included as inherent consideration in the NLC algorithm
design process.

## Acknowledgment

The authors wish to extend their thanks to Drs Rachid Bouziane and Philip
Watts, for enlightening discussions on the topic of hardware design, and
assistance with developing the models used in this work.

## References

  * [1] P. Bayvel, C. Behrens, and D. S. Millar, “Chapter 5 - Digital Signal Processing (DSP) and its application in optical communication systems,” in _Optical Fiber Telecommunications (Sixth Edition)_ , 6th ed., ser. Optics and Photonics. Boston: Academic Press, 2013, pp. 163–219.
  * [2] A. Napoli, Z. Maalej, V. A. J. M. Sleiffer, M. Kuschnerov, D. Rafique, E. Timmers, B. Spinnler, T. Rahman, L. D. Coelho, and N. Hanik, “Reduced complexity digital back-propagation methods for optical communication systems,” _Journal of Lightwave Technology_ , vol. 32, no. 7, pp. 1351–1362, 2014.
  * [3] E. Ip and J. M. Kahn, “Compensation of dispersion and nonlinear impairments using digital backpropagation,” _Journal of Lightwave Technology_ , vol. 26, no. 20, pp. 3416–3425, 2008.
  * [4] E. Temprana, E. Myslivets, L. Liu, V. Ataie, A. Wiberg, B. Kuo, N. Alic, and S. Radic, “Two-fold transmission reach enhancement enabled by transmitter-side digital backpropagation and optical frequency comb-derived information carriers,” _Opt. Express_ , vol. 23, no. 16, pp. 20 774–20 783, Aug 2015.
  * [5] D. Lavery, D. Ives, G. Liga, A. Alvarado, S. J. Savory, and P. Bayvel, “The Benefit of Split Nonlinearity Compensation for Single-Channel Optical Fiber Communications,” _IEEE Photonics Technology Letters_ , vol. 28, no. 17, pp. 1803–1806, 2016.
  * [6] D. S. Millar, S. Makovejs, C. Behrens, S. Hellerbrand, R. I. Killey, P. Bayvel, and S. J. Savory, “Mitigation of fiber nonlinearity using a digital coherent receiver,” _IEEE Journal on Selected Topics in Quantum Electronics_ , vol. 16, no. 5, pp. 1217–1226, 2010.
  * [7] G. P. Agrawal, _Applications of Nonlinear Fiber Optics (2nd Edition)_ , 2nd ed. Elsevier, 2008.
  * [8] L. B. Du and A. J. Lowery, “Improved single channel backpropagation for intra-channel fiber nonlinearity compensation in long-haul optical communication systems.” _Optics express_ , vol. 18, no. 16, pp. 17 075–17 088, 2010.
  * [9] M. Secondini, D. Marsella, and E. Forestieri, “Enhanced split-step Fourier method for digital backpropagation,” _European Conference on Optical Communication, ECOC_ , no. 1, pp. 3–5, 2014.
  * [10] M. Secondini, S. Rommel, G. Meloni, F. Fresi, E. Forestieri, and L. Potì, “Single-step digital backpropagation for nonlinearity mitigation,” _Photonic Network Communications_ , vol. 31, no. 3, pp. 493–502, 2016.
  * [11] Y. Gao, J. H. Ke, J. C. Cartledge, and S. Yam, “Low-Pass Filter Assisted Digital Back Propagation Algorithm for 112 Gb/s DP 16-QAM,” _Paper SPT5D.5 in Proc. Signal Processing in Photonic Communications (SPPCom) 2013_.
  * [12] C. Fougstedt, M. Mazur, L. Svensson, H. eliasson, M. Karlsson, and P. Larsson-Edefors, “Time-domain digital back propagation: Algorithm and finite-precision implementation aspects,” in _Optical Fiber Communication Conference_. Optical Society of America, 2017, p. W1G.4.
  * [13] S. J. Savory, “Digital filters for coherent optical receivers,” _Optics Express_ , vol. 16, no. 2, 2008.
  * [14] J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex fourier series,” _Mathematics of Computation_ , vol. 19, no. 90, 1965\.
  * [15] A. Alvarado, E. Agrell, D. Lavery, R. Maher, and P. Bayvel, “Replacing the soft-decision FEC limit paradigm in the design of optical communication systems,” _Journal of Lightwave Technology_ , vol. 33, no. 20, pp. 4338–4352, Oct 2015.
  * [16] D. Marcuse, C. R. Menyuk, and P. K. A. Wai, “Application of the Manakov-PMD equation to studies of signal propagation in optical fibers with randomly varying birefringence,” _Journal of Lightwave Technology_ , vol. 15, no. 9, pp. 1735–1745, 1997.
  * [17] M. Secondini and E. Forestieri, “On XPM mitigation in WDM fiber-optic systems,” _IEEE Photonics Technology Letters_ , vol. 26, no. 22, pp. 2252–2255, 2014.
  * [18] J. G. Proakis and D. G. Manolakis, _Digital Signal Processing_ , 4th ed. Pearson Education International, 2007.
  * [19] S. A. Khan, _Digital Design of Signal Processing Systems: A Practical Approach_. John Wiley & Sons, Ltd, 2011\.
  * [20] Y. Gao, J. H. Ke, J. C. Cartledge, and S. S. H. Yam, “Method for determining the low-pass filter bandwidth for the low-pass filter assisted digital back propagation algorithm,” in _39th European Conference and Exhibition on Optical Communication (ECOC 2013)_ , Sept 2013, pp. 1–3.
  * [21] R. Asif, C. Y. Lin, M. Holtmannspoetter, and B. Schmauss, “Logarithmic step-size based digital backward propagation in N-channel 112Gbit/s/ch DP-QPSK transmission,” in _2011 13th International Conference on Transparent Optical Networks_ , June 2011, pp. 1–4.
  * [22] C. Häger and H. D. Pfister, “Nonlinear interference mitigation via deep neural networks,” in _Optical Fiber Communication Conference_. Optical Society of America, 2018, p. W3A.4.