11institutetext: Stevens Institute of Technology, Hoboken NJ 07030, USA,
22institutetext: University of North Carolina at Chapel Hill, Chapel Hill NC
27599, USA

# Augmented Degree Correction for Bipartite Networks with Applications to
Recommender Systems

Benjamin Leinwand  11 Vladas Pipiras  22

###### Abstract

In recommender systems, users rate items, and are subsequently served other
product recommendations based on these ratings. Even though users usually rate
a tiny percentage of the available items, the system tries to estimate
unobserved preferences by finding similarities across users and across items.
In this work, we treat the observed ratings data as partially observed, dense,
weighted, bipartite networks. For a class of systems without outside
information, we adapt an approach developed for dense, weighted networks to
account for unobserved edges and the bipartite nature of the problem. This
approach allows for community structure, and for local estimation of flexible
patterns of ratings across different pairs of communities. We compare the
performance of our proposed approach to existing methods on a simulated data
set, as well as on a data set of joke ratings, examining model performance in
both cases at differing levels of sparsity. Keywords: Degree Correction, Link
Prediction, Bipartite Networks, Weighted Networks, Recommender Systems

##  1 Introduction

Recommender systems are widely used in industry to serve millions of users a
variety of products and services. The goal is to uncover a similarity
structure within and between sets of users and items that can be used to
identify specific items for those users who are expected to be the target
audience for those items. In general, user demographic features and item
features like genre are available. In this work, however, we focus on the
situation where the only information is users’ ratings of items, where each
user may rate a different subset of items. Ratings can be structured in a
matrix with users on one axis, items on the other axis, and each cell
containing the relevant user’s rating of the relevant item. For services with
large corpora and subscriber bases, most of the entries in this matrix will be
missing. A dataset can be treated as a partially observed, dense, weighted,
bipartite network.

The current state of the art in recommender systems uses variations of Graph
Neural Networks (GNNs) [1], and usually extensions of Graph Convolutional
Networks (GCNs) [2, 3]. These techniques incorporate features of both users
and items outside of the graph structure, like demographic or location
information, as well as the edges on the network itself. While GNNs appear to
succeed on a wide variety of tasks and graph structures, they also rely on
message passing, which takes advantage of the topology of the observed
network. In real-world systems, the existing edges may represent some signal,
as the items with which users have chosen to interact may reflect some
preference. For the specific task discussed in this paper, we focus on GCNs as
a feasible alternative to our model. GCNs learn an embedding for each user and
item, and a function mapping these embeddings to estimated ratings. User
defined hyperparameters include the embedding dimension of both users and
items, the number of layers, the number of nodes in each layer, which nodes
are connected, and the activation function of each edge, and the error metric.
A good model is one which, given a user and an item, can approximately predict
the value of the edge weight between the relevant user and item. Though we do
not account for it here, a possible addition to the GCN framework is to
include information about the users and items that is not contained in the
network. For example, genre of an item may be known, and can be incorporated
into the item’s embedding.

The task of estimating unobserved edge weights is reminiscent of matrix
completion, made famous by the Netflix Prize competition [4]. This topic has
spawned its own literature, evidenced by [5, 6] and the references therein. In
general, these works assume the matrix is low rank, and try to find the matrix
with the lowest rank that matches the observed entries, though [7] and [8]
attempt to minimize the nuclear norm in place of the rank. As such, we also
consider collaborative filtering techniques as potential alternatives to our
approach. However, in this setting, as the functions used to generate the
networks are nonlinear, this low rank assumption may not hold.

The method proposed in this paper can be seen as a nonlinear alternative
approach to matrix completion. By segmenting users and items into communities,
we can locally fit models to individual subnetworks, instead of the whole data
set. We adapt a model for the creation and estimation of dense weighted
networks that incorporates community structure and flexible connectivity
patterns, known as “augmented degree correction,” to this setting.
Additionally, we introduce an algorithm that accounts for missing values and
utilizes the bipartite structure of the network, performing community
detection on the user nodes and item nodes simultaneously. We apply our
methodology to a subset of the Jester Dataset [9] and to simulated data.

The rest of this paper is structured as follows. Section 2 describes the
proposed model, a clustering measure, a community detection algorithm, and an
estimation procedure suited to this setting. Section 3 details alternative
modeling algorithms. Section 4 defines the approach for measuring model
performance. A simulation study is conducted in Section 5, and the dataset and
results are discussed in Section 6. Appendix 0.A provides a deeper explanation
of our proposed community detection algorithm.

##  2 Proposed Model, Measure, and Algorithm

We approach the problem of modeling recommender systems by extending the work
in [10] to bipartite networks, to account for missing values, and sometimes
very few distinct values (e.g. when ratings only take on values in
$\\{1,2,3,4,5\\}$). We shall discuss the relevant portions of the approach,
but readers are referred to [10] for a more theoretical treatment of that
model’s motivation. The rationale for applying this approach mirrors the
proposed benefits of that work, in particular to allow for community
structure, nearly arbitrary distinct distributions of edge weights within and
between communities, and “augmented degree correction” characterized by
flexible “sociability” patterns (that is, the tendency for edges joining
particular nodes to have persistently higher edge weights resulting in
possibly non-convex contours demonstrating degree correction effects within
the network) through the use of “$H$-functions.”

The function $H_{ij}:(0,1)\times(0,1)\rightarrow(0,1)$ is an $H$-function if:

  * •

it is monotonic in both arguments

  * •

$\iint_{H(x,y)\leq z}dxdy=z$, for all $z\in(0,1).$

That is, if inputting two i.i.d. U$(0,1)$ random variables into an
$H$-function, the output is distributed as one U$(0,1)$ random variable.
$H$-functions have “positive association” if $H$ is non-decreasing in both
arguments. [10] suggests a recipe for devising a subset of $H$-functions by
taking cumulative distribution functions $F_{1}$ and $F_{2}$, as well as their
convolution $F_{1,2}$, and setting

$H(x,y)=F_{1,2}(F_{1}^{-1}(x)+F_{2}^{-1}(y)).$ (1)

By construction, any function using this recipe satisfies the conditions of
being a positive $H$-function The choice of $F_{1}$ and $F_{2}$ can define the
concavity of the shape and tilt of contours of $H$ on $(0,1)$, allowing for a
a variety of different maps. In the simulation study in Section 5, we generate
a network using $H$-functions with positive association, as well as those with
negative association. $H(x,y)$ is an $H$-function with negative association if
$H(1-x,1-y)$ is an $H$-function with positive association.

In the $H$-Normal Nonlinear Sociability Model (NSM) in [10], nodes $u$ are
assigned communities $i=1,\dots,K$, and “sociability” parameters $\Psi_{u}(k)$
with respect to all communities $k=1,\dots,K$, including their own, where all
$\Psi$ values are assumed to be i.i.d. U$(0,1)$ random variables. The original
model is not tailored to bipartite networks, so in the present situation,
nodes $u$ will be considered users (or row nodes), while $v$ will be
considered items (or column nodes). As a bipartite network, any edge must
connect exactly one user with exactly one item, so we need only be concerned
with modeling those kinds of edges. Consequently, $u$ may only belong to a
user (or row) community $i$, and $v$ may only belong to an item (or column)
community $j$. Each user $u$ is assigned “sociability” parameters
$\Psi_{u}(j)$ with respect to each item (column) communities $\\{j\\}$, while
each item $v$ is assigned “sociability” parameters $\Psi_{v}(i)$ with respect
to each user (row) communities $\\{i\\}$, Each user/item pair of communities
has their own possibly unique edge weight distribution $G_{ij}$, and a
possibly unique $H$-function that maps the two incident nodes’ $\Psi$ values
to a single value, as well as a parameter $\sigma_{ij}$ which modulates how
much the underlying sociability model is corrupted by random fluctuation. Edge
weights $W_{uv}$ are generated by mapping both nodal parameters and the edge
weight distribution into “standard normal” space. In particular, edge weights
are drawn according to the following formula:

$\Phi_{1}^{-1}(G_{ij}(W_{uv}))=\frac{1}{\sqrt{1+\sigma_{ij}^{2}}}\Phi_{1}^{-1}(H_{ij}(\Psi_{u}(j),\Psi_{v}(i)))+\frac{\sigma_{ij}}{\sqrt{1+\sigma_{ij}^{2}}}\epsilon_{uv},$
(2)

where $\epsilon_{uv}$ is an idiosyncratic $\mathcal{N}(0,1)$ random variable,
and $\Phi_{1}$ is the CDF of a $\mathcal{N}(0,1)$ distribution. By
construction, the right-hand side of (2) is a $\mathcal{N}(0,1)$ random
variable, so $W_{uv}$ has distribution $G_{ij}$. Letting $\sigma_{ij}=0$ would
mean repeated sample networks from the same model would all be identical,
while letting $\sigma_{ij}\rightarrow\infty$ would imply that each edge is a
purely random draw from the relevant distribution $G_{ij}$, where the
sociabilities of the incident nodes are immaterial.

###  2.1 Community Detection

Due to the network structure, instead of estimating just one set of
communities, users and items must each be partitioned into communities.
Communities in this context are defined differently than in other work such as
variants of the stochastic block model. Nodes in the same community $i$ all
share common $G_{ij}$, $H_{ij}$ and $\sigma_{ij}$ for each particular
corresponding community $j$. Instead of looking for homophily, patterns in
edge weights can be used to cluster nodes into communities: users in the same
row community $i$ share an “order of preferences” over items in a particular
column community $j$, as reflected by persistently greater edge weights.
Intuitively, a certain subset of users (e.g. those who love horror movies) may
generally share an order of preferences over a certain subset of items (e.g.
horror movies), while a different subset of users (e.g. those who hate horror
movies), may share a completely different order of preferences over the same
subset of items. This is different than describing item quality per se, as the
different subsets of users may demonstrate diametrically opposite preferences,
but it may broadly reflect latent features of the subsets of users and items
(e.g. whether each user enjoys being scared, and the scariness of each movie).
Perhaps less intuitively, though items do not have preferences, in this model,
items in the same column community $j$ also share an “order of preferences”
over users in a particular column community $i$, as reflected by persistently
greater edge weights. In this sense, the user and item communities are
interdependent, as fixing different partitions of users might induce different
partitions of items, and vice versa. Let

$d_{i}(v)=\sum_{u^{\prime}:\,u^{\prime}\in i}W_{u^{\prime}v}\,,\,v\in j,$

be the local degree of node $v$ in community $j$ with respect to community
$i$. Define a node-community correlation as

$C_{ij}(u)=corr\\{(d_{i}(v),W_{uv}):v\in j\\},$

where missing values are ignored. The goal of a metric for community detection
is to reward community structures where $C_{ij}$ values are consistently high
and positive, while also preferring larger communities where possible so as to
avoid overfitting too locally to only a few available edges.

To account for the bipartite structure and missing values, we must generalize
the measure $L$ introduced in [10] to the context of recommender systems. In
this case, let

$\displaystyle
L(\\{\,\widehat{i}\,\\},\\{\,\widehat{j}\,\\},\widehat{K_{r}},\widehat{K_{c}})=$
(3)
$\displaystyle\sum_{\widehat{i}=1}^{\widehat{K_{r}}}\sum_{\widehat{j}=1}^{\widehat{K_{c}}}\left(\overline{C_{\widehat{i}\,\widehat{j}}(u)}\times\left(1-\sqrt{SD(C_{\widehat{i}\,\widehat{j}}(u))}\right)+\overline{C_{\widehat{j}\,\widehat{i}}(v)}\times\left(1-\sqrt{SD(C_{\widehat{j}\,\widehat{i}}(v))}\right)\right)\times((n_{\widehat{i}}-2)(n_{\widehat{j}}-2))_{+}\times
density(\widehat{i},\widehat{j}),$

where $\widehat{K_{r}}$ and $\widehat{K_{c}}$ respectively represent the
numbers of estimated row (user) communities and column (item) communities,
$n_{\widehat{i}}$ is the number of nodes in community $\widehat{i}$, and
$density(\widehat{i},\widehat{j})$ represents the proportion of edges present
in the subnetwork containing only row nodes in estimated row community
$\widehat{i}$ and column nodes in estimated row community $\widehat{j}$.
$\overline{C_{\widehat{i}\,\widehat{j}}(u)}$ and
$SD(C_{\widehat{i}\,\widehat{j}}(u))$ refer to the average node-community
correlation over nodes in community $i$ with respect to community $j$, and the
standard deviation of those same values, respectively. In this form, $L$
values the row community $i$’s preference consistency over the column
community $j$ equally to $j$’s preference consistency over $i$. The reasons to
update the measure in [10] are that the matrix is no longer symmetric nor
square, and that the missing values imply we should not merely account for the
size of each subnetwork, but also for the subnetwork density; otherwise nodes
with consistent observed preferences over only a few observed edges might have
undue influence on the recovered community structures relative to those nodes
with marginally less consistent preferences observed over a much greater set
of incident edges. In a case where edges are not missing completely at random
(MCAR), the very fact of greater (or lesser) interaction between users and
items may itself be indicative of community behavior. This measure is the same
as the measure in the original $H$-Normal NSM context from [10], except with
respect to within community subnetworks. In this case, as the network is
bipartite, there are no within community subnetworks.

The algorithm we propose to maximize this measure $L$ is described in detail
in Appendix 0.A. It begins with a greedy algorithm to cluster singleton nodes
into communities, which is a more complex version of Algorithm 2 in the
supplement to [10]. Instead of clustering both rows and columns
simultaneously, this algorithm must alternate between rows and columns,
agglomerating sets of singletons into larger communities. This initial
agglomeration is especially difficult as $L$ depends on the interplay of row
community detection and column community detection, but at the outset, there
is no known community structure, and there is path dependency whereby the
initial choices of which nodes to merge into communities can reverberate
through the rest of the agglomeration step. That is to say, there is the least
useful information available to the greedy algorithm at the beginning, but
these earliest choices may have the greatest impact on the final recovered
communities. As a result, the second step in the algorithm must be more
complicated. It looks for nodes which have potentially been misclustered, and
tries several different approaches to putting those nodes into alternative
communities, and possibly combining communities. However, as community
structure is dependent on other nodes in all communities, the algorithm
alternates removing nodes from communities and recombining communities,
continuing to iterate as long as the measure $L$ continues to increase after
each cycle. At each cycle, the algorithm greedily seeks to maximize the
increase in $L$, but as the possible proposed partitions are limited by the
current configuration, the algorithm may attempt several possible partitions,
with the aim of avoiding very local minima near the current estimated
configuration.

###  2.2 Model Parameter Estimation

Once the communities have been estimated, the methodology proposed in [10] can
be directly applied to learn model parameters. The estimation procedure can
treat all subnetworks based on pairs of communities in this setting as
“between community” subnetworks. In this section we will review that process,
but as estimation is done at the subnetwork level without leakage from one
subnetwork to another, all $i$ and $j$ subscripts are dropped for convenience.
In short, $G$ is estimated using a modified version of a nonparametric
empirical CDF estimator, by simply ordering the edge weights and mapping this
order to a uniform distribution. $\Psi$ values are estimated by ordering the
values $\sum_{v\in j}\Phi^{-1}(\widehat{G}(W_{uv}))$ within each community
$i$. that is, by mapping the edge weights from their original distribution to
a standard normal distribution, and ordering the nodes by their local degrees
with respect to these transformed edge weights. The estimate of user $u$’s
$\Psi$ value with respect to community $j$ is denoted by
$\widehat{\Psi}_{u}^{(j)}$, while the estimate of item $v$’s $\Psi$ value with
respect to community $i$ is denoted by $\widehat{\Psi}_{v}^{(i)}$. $H$ and
$\sigma$ for each pair of communities is estimated by optimization using the
following equation:

$\displaystyle\\{\widehat{H},\widehat{\sigma}\\}=$
$\displaystyle\underset{H\in\mathcal{H},\sigma\geq
0}{\rm{argmin}}\sum\limits_{u\in i,v\in
j}\left(\Phi_{1}^{-1}(\widehat{G}(W_{uv}))-\frac{1}{\sqrt{1+\sigma^{2}}}\Phi_{1}^{-1}(H(\widehat{\Psi}_{u}^{(j)},\widehat{\Psi}_{v}^{(i)}))\right)^{2}.$
(4)

where $\mathcal{H}$ is a set of $H$-functions. Finally, to estimate the edge
weights in the network,

$\widehat{W}_{uv}=\underset{w}{\rm{argmin}}\left|\Phi_{1}^{-1}(\widehat{G}(w))-\frac{1}{\sqrt{1+\widehat{\sigma}^{2}}}\Phi_{1}^{-1}(\widehat{H}(\widehat{\Psi}_{u}^{(j)},\widehat{\Psi}_{v}^{(i)}))\right|.$
(5)

When not all edges are observed in the network, the network can be estimated
by ignoring missing edges, and subsequently replacing the missing values with
their estimate based on the incident nodes. The network can be iteratively re-
estimated until the estimated edge weights converge. In this paper, we do not
re-estimate until convergence, but instead use 10 iterations.

##  3 Alternative Approaches

Though Graph Neural Networks (GNNs) are a popular approach in recommender
systems, they are not well suited to this task, as they implicitly rely on
network topology as a source of information. That assumption is not
unreasonable, as the missingness of an edge may reflect underlying
preferences. For example some users may consciously avoid certain otherwise
popular items, as they may not expect the item to suit their taste. However,
the simulations in this paper operate under the assumption that edges are
missing completely at random (MCAR). Under the MCAR assumption, incorporating
network topology is spurious, as the presence or absence of an edge has no
correspondence to item or user characteristics, even if the weight of that
edge does. This also appears to be the case with the Jester dataset. As the
underlying network is assumed to be dense, but merely partially observed, the
topology of the network with all edges observed would only be a function of
the edge weights, not adjacencies. Even so, we attempted to fit a model using
message passing, but its poor performance did not warrant further
investigation.

Graph Convolutional Networks, on the other hand, appeared to show greater
promise in data rich simulations. On the simulations below, we embedded both
the users and items into 8-dimensional space, with 2 fully connected hidden
layers made up of 64 and 32 nodes, where all activation functions are RELU,
using MSE as the metric of interest. After some experimentation, these
settings were chosen because they successfully captured the dynamics on our
toy network when a large proportion of the edges are observed. Flexibility
appears to be a strength of GCNs, as this model was not purpose built for
bipartite, dense, weighted networks. It is also not explicitly identifying
communities at all, though those may be recoverable in the estimated
embeddings.

We also tested several algorithms from the recommenderlab R package [11] on
our simulations and dataset. Specifically, we include: Item-based
collaborative filtering (IBCF), User-based collaborative filtering (UBCF),
Latent Factor Models using Singular Value Decomposition (SVD) with column-mean
imputation, and the so-called Funk SVD using gradient descent. For more
information on each of these algorithms, refer to [12] and the references
therein. We do not include results from UBCF and SVD since they are always
competitive with but inferior to IBCF results. In all cases, we train the
algorithm on the entire training set using leave one out cross validation, and
all results for these algorithms are based on the average of the resulting
metric over for that algorithm over 10 runs. Other algorithms available in
recommenderlab were ignored due to runtime, performance, or appropriateness
for the task.

##  4 Model Evaluation

Using our proposed pipeline from Section 2 in the setting where only a single
network is observed, comparing different models’ performances is not
straightforward. Instead, there are several options.

Hold out edges: Treat a subset of the observed edges as missing values by
withholding them during training. Communities can be estimated based on the
observed edges, and the model can estimate the held out edges, which can be
compared to the true, unseen ratings.

Hold out users and items: Treat a subset of both users and items as missing,
and train the model on the remaining data. After training, assign each held
out user and item to the community that maximizes the measure in (3). Estimate
the observed edges based on the model fit on the training set.

Hold out users, items, and edges: Treat a subset of both users and items as
missing, and train the model on the remaining data. After training, among the
held out users and items, continue holding out a subset of edges, and assign
each to the community that maximizes the measure in (3). Finally, estimate the
missing edges using the model trained on only the training set.

Each of these options may be appropriate depending on the context. For
example, when the sets of users and items are known and fixed, holding out
edges alone may suffice. The second case may be helpful for understanding
whether the community detection results generalize to new observations, as
when new users join or items are added to the corpus. When new items and users
are added to the dataset, but one can ask existing users about new items and
new users about existing items, it may be best to hold out users, items, and
edges. In all cases, assuming enough observations are present, the data can be
split into training, validation, and test sets by holding out additional
edges, users, and items as necessary. This may be useful for finding the best
model, including choosing the best data transformation.

Unfortunately, only holding out edges alone works with the simplest versions
of GCNs, as the other options must still use newly observed edges to generate
embeddings for previously unobserved users and items, thereby contaminating
the estimates of the edge weights using the edges in the test set. This could
be avoided if the embeddings were generated based on features of the users and
items, as is the case for GraphSAGE [13], rather than strictly based on
ratings. This would also avoid the so-called “cold start” problem, when there
are no edges observed for a new user or item, but outside information is
available. Neither our model nor vanilla GCNs are built to handle this
problem, but conversely, GraphSAGE is not appropriate when the only available
information is contained in the network itself. In this work, we hold out
edges to compare models, but wish to explore the other methods in the future.

##  5 Simulation Study

The simulated dataset used for examining performance of the proposed
algorithms is shown in the top left plot of Figure 2. It consists of 3 item
communities and 4 user communities. All communities have 73 nodes, and have
nodes with $\Psi$ values equally spaced between .05 and .95, and all use the
$H$-function as in (1) where $F_{1}$ and $F_{2}$ have the gamma density

$f(x)=\begin{cases}\frac{-x^{-.5}e^{x}}{\Gamma(.5)},&x<0,\\\
0,&\text{otherwise},\end{cases}$

and $F_{1,2}$ is given by

$F_{1,2}(z)=\begin{cases}e^{z},&z<0,\\\ 1,&z\geq 0.\end{cases}$

When the subnetwork of interest is on the diagonal, there is positive
association, and the edge weights are distributed according to a U$(0,200)$
distribution. When the subnetwork of interest is on the off-diagonal, there is
negative association, and the edge weights are distributed according to a
U$(0,100)$ distribution. For all other subnetworks, there is negative
association, and the edge weights are distributed according to a U$(0,50)$
distribution. When comparing different models, the whole dataset will not be
observed during the training period. Instead, some subset of the dataset will
be hidden, which can then be used as a test set on which to measure model
performance. Of particular interest is the scenario when very few data points
are available, as in many real world settings for recommender systems.

One shortcoming of this simulation is determining which edges will be
observed, and which will be hidden. For our simulations, we conveniently
assume edges are MCAR, which is unlikely to occur in real datasets,
potentially due to different forms of selection bias. First, users know their
own tastes and can often read reviews before trying an item, so higher valued
edges may be more likely to be observed. At the same time, users may be more
likely to leave a rating for items which elicited strong feelings, whether
they were positive or negative. There may also be some interaction with a
recommender system algorithm, where users believe their ratings will shape
future recommendations. The combination of these behaviors would expect to
result in a relatively large number of extreme valued edges for many users,
but relatively few edges with moderate values. It may be desirable to sample
edges in a fashion consistent with these behaviors in future work, as well as
to use real datasets that are known to feature these patterns of omissions.

“Our proposed method” refers to the pipeline described in Section 2, including
both community detection and estimation (though we also include the case where
communities are known in advance). When the estimation procedure observes and
estimates a network with missing values, the unobserved values are replaced
with their estimates, then the network is re-estimated. Each network (except
the fully observed network) is estimated a total of 10 times, though not
necessarily until convergence, and the tenth estimate is treated as the final
estimate. It may be of interest to feed the network with estimated missing
edges back into the community detection algorithm, though this was not done.
For all GCN results, we fit the model described in Section 3 for 1000 epochs
and a batch size of 2048 in all epochs.

###  5.1 Simulation Performance and Comparison

A total of 21 networks were generated. The first network is the full network,
with no edges missing, so the only meaningful performance reported is that of
the GCN on the training set, if only to demonstrate that the choice of
hyperparameters can suitably capture the network dynamics. The rest of the
networks are 4 sets of 5 networks that have, respectively, 50%, 70%, 80%, and
90% of the edges deleted at random. This model has $\sigma=0$, though the
interaction between missing values and the error magnitude may be worthy of
investigation.

Our proposed method has two modules: community detection and estimation. The
estimated communities will impact the estimation performance, so we seek to
examine the performance of both modules, as well as the overall performance of
the combined methodology as the number of missing edges increases.The left
plot in Figure 1 shows the Normalized Mutual Information (NMI) of the
community memberships and their estimates. The community detection algorithm
performs perfectly in all cases when at most 70% of the edges are missing. The
model also performs well when 80% of the edges are missing. However, community
detection fails when 90% of edges are missing. We suspect the issue is not due
to the proportion of missing values, but rather the lack of observed edges
incident to each node. If we kept the same sparsity but increased the size of
the network to observe some threshold number of total edges incident to each
node, we would expect the community detection to succeed. We have begun to
test this hypothesis, first by removing every other user and item node from a
network missing 50% of its edges. In this case, community detection failed.
Next, we took a network missing 80% of its edges where community detection
failed to recover the true communities, and duplicated all the nodes, thus
creating a network with 4 times as many edges, but no new information. In this
case, true communities were recovered.

Figure 1: Community detection performance of our method, and comparing
estimation performance with GCN.

Figure 2: By column from left to right: observed networks, our model’s
estimates, the average of 5 model estimates, GCN estimates, and the average of
5 GCN estimates.

The left column of Figure 2 shows the network the model observes, where white
spaces represent unobserved edges. The network in the first row contains all
edges, while the next 4 rows are missing 50%, 70%, 80%, and 90% of the edges,
respectively. In the second column from the left, the results of the
estimation procedure of our proposed method. In the middle column, the average
of 5 estimated networks, where each estimate is based on a different observed
network which contains the same proportion of missing edges as the observed
network in that row. In cases with many missing edges, the model may expect
too much from nodes that randomly are not missing their largest edges.
Furthermore, the model may not have a large number of different edge weights
to choose from within each subnetwork, potentially resulting in resolution
loss. Even with much of the data missing, the model recovers the broader
pattern, though edges incident to specific nodes may be overestimated or
underestimated, resulting in disordered or “plaid” looking patterns. However,
when community detection fails, so too does estimation, as the pattern of
preferences is more muted due to the mixing of positive and negative
associations for nodes across the network.

The GCN estimates appear to degrade even more sharply than using our proposed
method. The right plot in Figure 1 shows the RMSE of the estimates for the 21
different networks using the estimates from GCNs, our proposed method
including community detection, and our proposed method with known communities.
While all fit the fully observed network nearly perfectly, a divergence begins
when 70% of the data is missing, and the 3 methods are clearly in separate
clusters when 90% of the data is missing, with GCN performing the worst. This
plot shows the RMSE of edges in the test set, but when the whole network is
observed, there is no test set, in which case the RMSE on the training set is
reported. The right two columns of Figure 2 depict the estimates and averages
from the fitted GCNs for different levels of sparsity, which look somewhat
worse than our method’s estimates.

Figure 3: Estimates of 90% missing network when using the true community
memberships.

The analysis indicates that community detection is pivotal. If the true
communities are recovered, the model can broadly reconstruct connectivity
patterns, even when most edges are missing, as illustrated in Figure 3. The
left plot of that figure displays the estimated 90% missing network shown in
Figure 2 using the true communities. The right plot of Figure 3 represents the
average of 5 estimates using the true communities, rather than the estimated
communities. When there are few nodes and many edges are missing, though,
community detection may remain a difficult problem, simply because there isn’t
enough information in the dataset. Recycling nodes to increase absolute
coverage appears to be a promising option, though when error is included (that
is, when $\sigma>0$), that error will also be amplified.

Missing Proportion |  NMI of Row Community Estimates |  NMI Column Community Estimates |  Proposed Method Test Set MSE |  Proposed Method with Known Communities Test Set MSE |  GCN Test Set MSE |  IBCF Test Tet MSE |  Funk SVD Test Set MSE
---|---|---|---|---|---|---|---
0 |  1 |  1 |  1.87 (training set) |  1.87 (training set) |  3.0 (training set) |  NA |  NA
0.5 |  1 |  1 |  6.44 |  6.44 |  9.91 |  55.42 |  5107.77
0.5 |  1 |  1 |  6.62 |  6.62 |  9.45 |  54.94 |  5114.18
0.5 |  1 |  1 |  5.96 |  5.96 |  15.23 |  53.70 |  5127.62
0.5 |  1 |  1 |  6.30 |  6.30 |  4.88 |  53.68 |  5099.26
0.5 |  1 |  1 |  6.37 |  6.37 |  12.14 |  57.24 |  5117.81
0.7 |  1 |  1 |  20.44 |  20.44 |  60.26 |  72.32 |  5110.85
0.7 |  1 |  1 |  18.12 |  18.12 |  48.36 |  69.59 |  5143.52
0.7 |  1 |  1 |  22.40 |  22.40 |  29.57 |  73.81 |  5124.30
0.7 |  1 |  1 |  21.92 |  21.92 |  33.42 |  73.54 |  5091.43
0.7 |  1 |  1 |  20.70 |  20.70 |  19.15 |  71.55 |  5162.01
0.8 |  0.98 |  1 |  52.36 |  51.16 |  278.42 |  125.56 |  5137.09
0.8 |  1 |  1 |  47.56 |  47.56 |  206.19 |  116.28 |  5103.40
0.8 |  0.96 |  1 |  58.87 |  52.61 |  160.97 |  127.31 |  5122.81
0.8 |  0.99 |  1 |  50.73 |  50.12 |  192.34 |  132.46 |  5086.41
0.8 |  0.99 |  1 |  50.43 |  49.55 |  83.82 |  117.93 |  5111.58
0.9 |  0.49 |  0 |  1667.46 |  152.13 |  2262.41 |  1532.01 |  336.04
0.9 |  0.50 |  0 |  1736.39 |  145.21 |  2004.16 |  1473.11 |  251.05
0.9 |  0.44 |  0.03 |  1129.04 |  169.02 |  1906.90 |  1567.24 |  311.70
0.9 |  0 |  0.45 |  1806.97 |  167.90 |  1891.07 |  1513.37 |  390.95
0.9 |  0 |  0.46 |  1851.84 |  159.07 |  2057.29 |  1466.07 |  267.55
Table 1: Comparing the performance of our model against alternatives on the simulated data. Varying the proportion of hidden ratings, hiding 50%, 70%, 80%, and 90% of the data. Community detection performance is measured using the NMI between the recovered communities and the true communities. The best result for each run (based on test set MSE) is bolded. The gray column represents the test set performance of our estimation procedure if the true communities are known in advance. Approximate missing proportion |  Transformation |  Number of Row Communities detected |  Number of Column Communities detected |  NMAE for proposed model |  NMAE for GCN |  NMAE for GCN training set |  NMAE for IBCF (average of 10)
---|---|---|---|---|---|---|---
.445 |  none |  15 |  1 |  .1760 |  .1860 |  .1312 |  .1680
.445 |  none |  4 |  1 |  .1737 |  .1831 |  .1298 |  .1655
.445 |  Center Columns |  14 |  1 |  .1739 |  .1831 |  .1296 |  .1659
.63 |  none |  6 |  1 |  .1766 |  .1952 |  .1254 |  .1693
.63 |  none |  7 |  2 |  .1769 |  .1912 |  .1268 |  .1689
.63 |  none |  3 |  2 |  .1758 |  .1940 |  .1233 |  .1703
.815 |  none |  5 |  1 |  .1804 |  .2082 |  .09602 |  .1850
.815 |  none |  4 |  3 |  .1812 |  .2082 |  .0963 |  .1838
.815 |  none |  5 |  2 |  .1796 |  .2103 |  .0961 |  .1850
|  |  |  |  |  |  |
.8125 (75% dense with further 75% hidden) |  no transformation, warm start for the column communities |  3 |  4 |  .1762 |  .1960 |  .1028 |  0.1889
Table 2: Comparing the performance of our model against alternatives on the
Jester Dataset for users who rated exactly 74 jokes. Varying the proportion of
ratings which are hidden, hiding $\sim$ 25%, $\sim$ 50%, and $\sim$ 75% of the
74% dense network. As true communities are unknown, the number of recovered
row and column communities are reported by the best performing transformation
of the data according to training set cross validated NMAE. The best result
for each run based on test set NMAE is bolded. The gray column represents
training set performance of the GCN, indicating alternative hyperparameters
are unlikely to drastically improve GCN test set performance. The final row
contains results for a different dataset, users who rated exactly 75 jokes.

Table 1 includes detailed results from the simulation. Our proposed method
outperforms other algorithms when at most 80% of the edges are missing, and
the GCN is competitive when most 70% of the edges are missing. A few
interesting things occur when the network gets very sparse, though. First,
poor community detection means other approaches are preferred when 90% of the
edges missing. Even more surprisingly, while the Funk SVD algorithm performs
poorly when at least 20% of the network is observed, it is by far the best
performer with only 10% of the data, the only algorithm to improve with less
information.

##  6 Analysis of Jester Dataset

Using the Jester Data set from [14] of joke ratings from 24,983 users who have
rated at least 36 out of 100 jokes, we include only those users who rated
exactly 74 jokes. This was chosen to get a feasibly large dataset, and one
with a reasonable amount of missing data (the network is chosen to be 26%
missing). This leaves 626 users and 100 jokes. The jokes do not appear to be
MCAR, but instead many jokes are far more frequently rated than others, likely
due to ordering effects. Ratings run from $-10$ to 10. The distribution of
ratings for either users or items shouldn’t be expected to be i.i.d. Some
jokes may be broadly perceived to be funnier than others, and some users may
be more critical than others. It may be desirable to transform the data to
account for these effects, as it may be more informative to identify nodes
with similar relative order of preferences. For this reason, when performing
community detection only, we do not just use the observed matrix of ratings,
but also the matrix after column centering, after row centering, after column
normalizing, and after row normalizing. Column (row) centering means each
column (row) will have a mean of 0, while column (row) normalizing means each
column (row) will be demeaned and have a standard deviation of 1.

From the 74% dense dataset described above, we create 9 (non-independent)
experiments. In the first 3, we keep 75% of the edges (both observed and
missing) in the training set, and put the remaining 25% in the test set.
Accounting for the initial network’s density, the training sets have about
55.5% density. In the next 3, the train-test split is 50/50, so the training
set is about 37% dense. The final 3 examples, the train test split is 25/75,
with training set density about 18.5%. As the test sets are chosen randomly,
edges are expected to appear in multiple test sets. For each of these 9
training sets, we conduct 3 fold cross validation. For every fold, we estimate
communities based on the remaining 2 folds of the training set. Then, using
the estimated communities and the original values (not including the test set
or validation fold), we estimate the network, running 10 iterations to impute
and stabilize the missing values. These estimates are compared to the edge
weights in the validation fold. The transformation that yields the minimal
average Normalized Mean Absolute Error (NMAE) is used for community detection
for this experiment. Although MSE was the metric for the simulated data, we
use NMAE here, following [14]. As the range of the data is $-10$ to 10, NMAE
is MAE divided by 20. Finally, communities are estimated using the entire
training set and the best performing transformation. Those communities are
used to estimate the edge weights in the test set, iterating 10 times.

Figure 4: By column from left to right: Observed networks, our model’s
estimates, the average of 3 model estimates, GCN estimates, and the average of
3 GCN estimates. For legibility purposes, all networks are reordered according
to the estimated communities as detected by the eighth row in Table 2, with 4
row communities and 3 column communities, which is depicted in the bottom plot
in the second column from the left.

Table 2 contains detailed results for these experiments. Unlike the simulation
dataset, true communities are unknown here, so NMI values can’t be computed.
Instead, we include the GCN performance on the training set, to demonstrate
that the GCN approach is likely overfitting to the training set. All models
perform worse with sparser training sets. Our model is best in the sparsest
case, while approaches from recommenderlab achieve better results when more
data is available. For this dataset, IBCF dominates all other algorithms
besides our proposed model on these experiments. These results indicate that
the proposed approach may not be well-suited to detecting subtler signals
restricted to a smaller subset of users and items, leaving room for
enhancement. Nonetheless, as this model can robustly recover large scale,
coarse patterns in the data, there appears to be a sparsity regime where the
proposed approach provides better than state-of-the-art performance. As real
world datasets are often extremely sparse, the sparsest setting may most
resemble realistic conditions, albeit without including additional features
about users and items.

Figure 4 shows the results of our approach and GCNs, with all networks and
predictions reordered by the communities detected in the run depicted in the
bottom row of Figure 4, which is why that plot contains smooth contours. Each
run detected different community structures, but a consistent ordering was
chosen in order to depict broader trends throughout the different results. The
ordering of the original dataset is not informative beyond the fact that
certain jokes have many more ratings than other (the less rated jokes look
like white strips in the top left plot in Figure 4). This doesn’t appear to be
related to the content of the jokes, and splitting the data into training and
test sets obviates some of this disparity. Though there is perhaps more fine
grained prediction with more data, our method’s predictions do not appear to
change too much even in sparse settings. The right two columns of Figure 4
show some GCN estimates appear to be more extreme. In contrast to our model
and GCN, recommenderlab only makes predictions on the test set, so plots would
be misleading as they would frequently include observations directly from the
training set.

Based on the resulting performance of the recovered communities on different
folds of the training set, we can determine which transformation appears to
find meaningful sets of communities. In practice, using the original dataset
led to better performance than any transformation in 8 of 9 experiments. As
our community detection algorithm is slow, cross validation took far longer
than any other aspect of the modeling pipeline. Our recommendation is
therefore to skip cross validation and simply perform community detection and
estimation on the original training set.

Finally, to ensure that this performance is not limited to the dataset of
users who rated exactly 74 jokes, we attempted a similar analysis with the
dataset of users who rated exactly 75 jokes. Though the jokes are the same,
this is a different dataset, as all users are previously unobserved. We
included only 25% of data in the training set, leaving an approximately 18.75%
dense network as the training set. Starting with a previously recovered set of
item (joke) communities from the eighth run in Table 2, but initializing all
users as singletons, we re-ran community detection without cross validation,
and ran our estimation procedure. In this case of “transfer learning,” our
model again outperformed other approaches, indicating the model’s success in
this sparsity regime is not simply a quirk of the previous dataset.

##  7 Future work

While our proposed methodology provides promising results, there is room for
improvement. In real world, online recommender systems, additional information
is often known about users or items which may be leveraged to improve
recommendations, especially for new users or items. One option would be to
generalize the model discussed here, but it might be simpler to adapt the
present notion of communities into a GCN, as was done with [15], which has
communities based on modularity. That option may combine the functional
fitting strengths of neural networks while restricting the fit within
estimated subnetworks. However, the community detection process is too slow
for industrial scale. Faster approaches are necessary, likely requiring
developments in both theory and algorithms. It is difficult to upgrade from a
“good” clustering to an optimal or near optimal one because the measure $L$ is
based on node-community correlations, and because poor user clusterings can
cascade into poor item clusterings or vice versa.

Though cross validation was used to try to find the best transformation of the
data, it ultimately proved unhelpful. However, cross validation could be used
to choose model complexity. One interpretation of the slower rate of
performance decay for our model is that our model is able to find coarse,
globally consistent patterns with few data points. However, given more data,
algorithms like IBCF may be able to pick up cross-cutting patterns that are
too fine for our method to detect. Some distinct communities may be
incorrectly merged based on similar large scale dynamics, displaying a kind of
resolution limit. Given our results, it may be useful to run the estimation
pipeline, and subsequently run it on the matrix consisting of residuals of the
initial estimates, allowing for more complicated behaviors. Cross validation
can be used to determine the appropriate model complexity. Ideally, this would
encourage a model to capture only simpler, global patterns when data is
sparse, while using more locally appropriate information when that signal can
be detected using a richer dataset, achieving excellent performance in both
regimes.

## References

  * [1] Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C. & Sun, M. Graph neural networks: A review of methods and applications. AI Open. 1 pp. 57-81 (2020)
  * [2] Hamilton, W., Ying, R. & Leskovec, J. Representation learning on graphs: methods and applications. IEEE Data Engineering Bulletin. 40, 52-74 (2017)
  * [3] Schlichtkrull, M., Kipf, T., Bloem, P., Berg, R., Titov, I. & Welling, M. Modeling relational data with graph convolutional networks. European Semantic Web Conference. pp. 593-607 (2018)
  * [4] Bennett, J., Lanning, S. & Others The Netflix prize. Proceedings Of KDD Cup And Workshop. 2007 pp. 35 (2007)
  * [5] Nguyen, L., Kim, J. & Shim, B. Low-rank matrix completion: A contemporary survey. IEEE Access. 7 pp. 94215-94237 (2019)
  * [6] Chen, Z. & Wang, S. A review on matrix completion for recommender systems. Knowledge And Information Systems. 64, 1-34 (2022,1)
  * [7] Candès, E. & Recht, B. Exact matrix completion via convex optimization. Foundations Of Computational Mathematics. 9, 717-772 (2009)
  * [8] Candès, E. & Tao, T. The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions On Information Theory. 56, 2053-2080 (2010)
  * [9] Jester Dataset, https://goldberg.berkeley.edu/jester-data
  * [10] Leinwand, B. & Pipiras, V. Block dense weighted networks with augmented degree correction. Network Science. 10, 301-321 (2022)
  * [11] Hahsler, M. recommenderlab: Lab for Developing and Testing Recommender Algorithms. (2023), https://CRAN.R-project.org/package=recommenderlab, R package version 1.0.3
  * [12] Hahsler, M. recommenderlab: An R Framework for Developing and Testing Recommendation Algorithms. (arXiv:2205.12371 [cs.IR],2022,5)
  * [13] Hamilton, W., Ying, Z. & Leskovec, J. Inductive representation learning on large graphs. Advances In Neural Information Processing Systems. 30 (2017)
  * [14] Goldberg, K., Roeder, T., Gupta, D. & Perkins, C. Eigentaste: A constant time collaborative filtering algorithm. Information Retrieval. 4 pp. 133-151 (2001)
  * [15] Tsitsulin, A., Palowitch, J., Perozzi, B. & Müller, E. Graph clustering with graph neural networks. Journal Of Machine Learning Research. 24, 1-21 (2023)

##  Appendix 0.A Algorithms

This appendix contains algorithms for community detection for partially
observed dense weighted bipartite networks. The overarching algorithm
displayed as Algorithm 1 is relatively simple, though the implementation is
technically complex.

All row (user) and column (item) nodes begin in their own singleton
communities, so the first stage is to alternate combining row communities and
column communities, as long as the suggested change improves the measure $L$
in Equation (3). This is quite similar to Algorithm 2 from the supplemental
technical appendix to [10], but has the extra wrinkle of having to account for
both row nodes and column nodes, and includes adaptations of the attemptMerge
and sweep procedures in Algorithms 4 and 5. When no more communities can be
combined, this first part of the algorithm, agglomeration, shown in Algorithm
LABEL:alg:agglom, terminates. While this is intended to estimate a reasonable
set of communities essentially from scratch, the issue is that row (or column)
communities that are combined early on in agglomeration are not accounting for
the column (or row) community structure that develops later. This is a
feature, though, not a bug. Community detection is a partitioning problem, and
NP-hard. As a means of getting any traction, agglomeration is designed only to
combine communities, not to consider how to break them down more efficiently.
Especially with a lot of missing values and no known community structure, it
is to be expected that agglomeration will cluster many nodes into the wrong
communities.

Result: $\\{\,\widehat{i}\,\\},\\{\,\widehat{j}\,\\}$ // estimated community
memberships for row (user) and column (item) nodes

Input: $W$ // the observed network

1 [rowLabels, columnLabels, rowCommunityAggregate, columnCommunityAggregate,
L] = agglom($W$, $\\{1,...,n_{r}\\},\\{1,...,n_{c}\\}$) // $n_{r}$ and $n_{c}$
are the number of rows (users) and columns (items), respectively

2 [$\\{\,\widehat{i}\,\\}$, $\\{\,\widehat{j}\,\\}$] = fixCommunities($W$,
rowLabels, columnLabels, L)

Algorithm 1 Greedy algorithm using $L$

Procedure  _agglom(_W, rowLabels, columnLabels_)_

1 $q$=1

2 rowLabels${}^{(0)}=\vec{0}_{n_{r}}$; rowLabels(1) = rowLabels

3 columnLabels${}^{(0)}=\vec{0}_{n_{c}}$; columnLabels(1) = columnLabels

4 rowCommunityAggregate = $W^{\prime}$

5 colCommunityAggregate = $W$ // $rowCommunityAggregate$ is a
$\widehat{K_{r}}\times n_{c}$ matrix that treats each row community as a
‘‘supernode,’’ containing aggregate edge weights to each column node.
$colCommunityAggregate$ does the same for column communities.

6 while  _$rowLabels^{(q)}\neq rowLabels^{(q-1)}\,||\,columnLabels^{(q)}\neq
columnLabels^{(q-1)}$_ do

7 $q$++

8 rowLabels(q) = rowLabels(q-1)

9 columnLabels(q) = columnLabels(q-1)

10 rowLabelOrder = sort labels by decreasing column sum of
rowCommunityAggregate

11 colLabelOrder = sort labels by decreasing column sum of
colCommunityAggregate

12 while  _$\max(length(rowLabelOrder),length(colLabelOrder)) >0$_ do

13 rowLabels(q), rowCommunityAggregate, rowLabelOrder =
attemptMergeBipartite(_$W$ , rowCommunityAggregate, rowLabels(q),
rowLabelOrder, columnLabels(q)_)

14 dequeue(rowLabelOrder(1)) // remove first entry

15 columnLabels(q), columnCommunityAggregate, columnLabelOrder =
attemptMergeBipartite(_$W^{\prime}$ , columnCommunityAggregate,
columnLabels(q), columnLabelOrder, rowLabels(q)_)

16 dequeue(colLabelOrder(1)) // remove first entry

end while

17if  _$rowLabels^{(q)}=rowLabels^{(q-1)}\,\
&\&\,columnLabels^{(q)}=columnLabels^{(q-1)}$_ then

// if no communities are combined using $attemptMergeBipartite$, try to
combine using $sweep$

18 rowLabels(q), rowCommunityAggregate = sweepBipartite(_$W$ ,
rowCommunityAggregate, rowLabels(q), columnLabels(q)_)

19 columnLabels(q), columnCommunityAggregate = sweepBipartite(_$W^{\prime}$ ,
columnCommunityAggregate, columnLabels(q), rowLabels(q)_)

end if

end while

20L = calculateMeasureL(W, rowLabels(q), columnLabels${}^{(q)})$

21 return rowLabels(q), columnLabels(q), rowCommunityAggregate,
columnCommunityAggregate, L

Algorithm 2 agglom

The goal of the second step, shown in Algorithm 3, is to identify and correct
any clustering errors. The thorny issue that remains is one of
interdependence. Even in perfect conditions with easily defined communities
and $\sigma=0$, if the column community assignments are wrong, a wrong
partition of row nodes may outperform the true communities in the measure $L$.
So while it’s difficult to immediately reclassify misclassified nodes, there
are heuristic approaches available, which can be used iteratively to get to
better and better partitions. Those heuristics are used in Algorithms 7, 8, 9,
and 10, which will be discussed, but at a higher level, Algorithm 3 starts
with the results of agglomeration, and checks many alternative clusterings by
using these heuristics, calling Algorithm 6. If the best performing
alternative clustering achieves a higher $L$ than that of agglomeration, this
new clustering is adopted and the process is repeated. If no heuristically
induced clustering outperforms the current estimated communities, the
algorithm takes the best performing result of the heuristics where all
communities are of size $>2$, and looks for alternative clusterings based on
that structure. If any outperform the current clustering, the best is adopted
and algorithm returns to the beginning. If not, as a last ditch effort, the
model starts with the current clustering and attempts to combine each pair of
row communities, then each pair of column communities. If the best performing
change results in a higher $L$ than the current clustering, it is adopted and
algorithm returns to the beginning. The whole purpose of all this machinery is
to deal with the interdependence, where an early error in agglomeration may
require several iterations to fix, yet to try to avoid being bogged down
attempting every possible change.

Procedure  _fixCommunities(_W, rowLabels, columnLabels, L_)_

1 oldMeasure = 0; newMeasure = L

while  _oldMeasure $<$ newMeasure_ do

2 oldMeasure = newMeasure

3 newMeasure = getLargestOneStepL(_W, rowLabels, columnLabels_)

4 if  _newMeasure(L) $>$ oldMeasure_ then

5 update rowLabels and columnLabels to those that produced newMeasure

end if

else

6 largerGroups = newMeasure(largerGroups)

7newMeasure = getLargestOneStepL(_W, largerGroups(rowLabels),
largerGroups(columnLabels)_)

8 if  _newMeasure(L) $>$ oldMeasure_ then

9 update rowLabels and columnLabels to those that produced newMeasure

end if

if  _newMeasure $\leq$ oldMeasure_ then

10 Combine each pair of rowCommunities, and each pair of columnCommunities. If
the best performing amalgamation produces a better $L$ value than oldMeasure,
update rowLabels, columnLabels, and newMeasure to those that produced this
higher $L$ value // As this requires ${\widehat{K_{r}}\choose
2}+{\widehat{K_{c}}\choose 2}$ comparisons, we only run this when all other,
more efficient heuristics for adjusting community assignments fail

end if

end if

end while

11return rowLabels, columnLabels

Algorithm 3 fixCommunities

Procedure  _attemptMergeBipartite(_W , communityAggregate, relevantLabels,
relevantLabelOrder, otherLabels_)_

1 newLabels = relevantLabels

2 newOrder = relevantLabelOrder

3 A = newOrder(1)

4 B = $\underset{\text{B}^{\prime}\neq\text{A}}{\rm{argmax}}$
$corr([\text{communityAggregate(A)}],[\text{communityAggregate(B')}])$

// [communityAggregate(A)] represents the Ath column in the communityAggregate
matrix (either rowCommunityAggregate or columnCommunityAggregate). Correlation
is the "pairwise complete" correlation here and everywhere else in this
document. This algorithm seeks to combine community ‘‘A’’ with the community
whose edge weights are most correlated with A’s

5 mergedLabels = newLabels

6 mergedLabels(mergedLabels = B) = A

7 mergedOrder = newOrder

8 mergedOrder(mergedOrder = B) = A

// The above four lines will tentatively merge the communities represented by
columns A and B to see if it improves the measure L below

9 if  _calculateMeasureL( W, mergedLabels, otherLabels) $\geq$
calculateMeasureL(W, newLabels, otherLabels)_ then

10 newLabels = mergedLabels

11 newOrder = mergedOrder

12 [communityAggregate(A)] += [communityAggregate(B)]

13 Remove [communityAggregate(B)] column from communityAggregate

// If the measure L is improvedby combining the communities, update both in
the community labels and in the relevant communityAggregate matrix

end if

14return newLabels, communityAggregate, newOrder

Algorithm 4 attemptMerge for bipartite networks

Procedure  _sweepBipartite(_W, communityAggregate, relevantLabels,
otherLabels_)_

1 correlationOrder = sort non-matching label pairs by decreasing correlations
of columns in communityAggregate with each row containing 2 communities

2 while  _$length(correlationOrder) >0$_ do

3 newLabels = relevantLabels

4 A = correlationOrder(1, 1) // First community in first row of
correlationOrder

5 B = correlationOrder(1, 2) // Second community in first row of
correlationOrder

6 mergedLabels = newLabels

7 mergedLabels(mergedLabels = B) = A

8 if  _calculateMeasureL(W, mergedLabels, otherLabels) $\geq$
calculateMeasureL(W, newLabels, otherLabels)_ then

9 newLabels = mergedLabels

10 [communityAggregate(A)] += [communityAggregate(B)]

11 Remove [communityAggregate(B)] column from communityAggregate

12 dequeueAll(correlationOrder)

else

13 dequeue(correlationOrder(1))

end if

end while

14return newLabels, communityAggregate

Algorithm 5 Sweep for bipartite networks

Procedure  _getLargestOneStepL(_W, rowLabels, columnLabels_)_

1 // Below, we write the assignment informally. Each assignment inherits the
resulting rowLabels, columnLabels, rowCommunityAggregate,
columnCommunityAggregate, and L, and possibly other results from the procedure
on the right hand side. When feeding in results from one procedure to another,
any necessary additional arguments are provided, but omitted here for brevity

a = removeNodesA(_W, rowLabels, columnLabels_)

2 b = regroupNodesA(_a, banRemain = 0_)

3 c = regroupNodesA(_a, banRemain = 1_)

4 d = agglom(_a_); e = agglom(_b_); f = agglom(_c_)

5 g = removeNodesB(_W, rowLabels, columnLabels_)

6 h = regroupNodesB(_g, banRemain = 0_)

7 i = regroupNodesB(_g, banRemain = 1_)

8 j = agglom(_g_); k = agglom(_h_); l = agglom(_i_)

9 m = regroupNodesB(_a, banRemain = 0_); n = agglom(_m_)

10return $\underset{a,b,\dots,n}{\rm{argmax}}(L)$,
$\underset{a,b,\dots,n}{\rm{argmax}}\\{L:min(\text{community size})>2\\}$ //
The first list of returned variables contains all features of the best
clustering according to the measure $L$. The second list of returned variables
contains all features of the best clustering according to measure $L$ where
there are no communities of size 1 or 2, and is thus referred to as
$largerGroups$. Both returned lists may be identical

Algorithm 6 getLargestOneStepL

As for these heuristics, the first step is identifying nodes which may be
misclassified. Going back to Section 2.1, we expect a good estimated
clustering to have large positive $C_{\widehat{i}\widehat{j}}$ values
everywhere. Algorithms 7 and 8 are built on the insight that if there is a
negative $C_{\widehat{i}\widehat{j}}$ value, then a row (or column) node’s
edges are anticorrelated with its assigned community’s edges, with respect to
some column (row) community, which may mean the row node’s comunity assignment
is misclassified. However, if all the negative row
$C_{\widehat{i}\widehat{j}}$ values are due to one particular column community
$j$, then perhaps that whole estimated column community is the root of the
problem. Algorithms 7 and 8 handle this issue somewhat differently. If all row
node issues are due to only one column community, Algorithm 7 does not
consider these row nodes to be misclassified, while Algorithm 8 still does.
Nevertheless, the purpose of each is to identify those problematic nodes and
remove their current community labels. Algorithm 7 sends any potentially
misclassified node to its own new singleton community. Algorithm 8 will remove
problematic nodes from their current communities, but will take any row nodes
that came from the same community that also have negative
$C_{\widehat{i}\widehat{j}}$ values with respect to the exact same row
communities, and place them together into a newly formed community. Especially
in the case where many edges are missing, having communities with multiple
nodes may make future steps more robust. However, because having exactly two
nodes in a community can make it very easy for correlations to be very high,
Algorithm 8 also breaks up any community of exactly two nodes into two
different communities. The reason for including both of these algorithms is
clustering may be unstable in myriad ways, so these approaches take different
paths, giving more options to find marginally better, or even much better
clusterings. Aggomeration is run on both of these results to see if the newly
“pure” communities can be recombined, but there are also two more algorithms
for recombination.

Procedure  _removeNodesA(_W, rowLabels, columnLabels_)_

1 rowCorrMat = $n_{r}\times\widehat{K_{c}}$ matrix where entry ($u,j$)
contains $C_{\hat{i}\hat{j}}(u)$ where $u\in\widehat{i}$

2 columnCorrMat = $n_{c}\times\widehat{K_{r}}$ matrix where entry ($v,i$)
contains $C_{\hat{j}\hat{i}}(v)$ where $v\in\widehat{j}$

// $n_{r}$ and $n_{c}$ represent the number of rows and columns in the
network, respectively. $rowCorrMat_{u,j}$ shows correlation between the edge
weights connecting row node $u$ to estimated column community $j$ and the
total degree of column nodes in estimated column community $j$ with respect to
estimated row community $i$. $columnCorrMat_{v,i}$ shows correlation between
the edge weights connecting column node $v$ to estimated row community $i$ and
the total degree of row nodes in estimated row community $i$ with respect to
estimated column community $j$. Any negative correlations may indicate a
clustering issue. This algorithm moves each problematic node to a new
singleton community

3 rowWrongNodes = $\\{\\}$; columnWrongNodes= $\\{\\}$

54 correspondingColumns = $\\{\\}$; correspondingRows = $\\{\\}$ // Every top
level for loop or if statement done for row nodes will be done correspondingly
immediately thereafter for column nodes. From here on, for simplicity, only
the row node steps are shown. for  _u in 1 to $n_{r}$_ do

if  _min( $rowCorrMat_{i}$)$<$0 or max($rowCorrMat_{i}$)= 0_ then

queue(rowWrongNodes, u)

queue(correspondingColumns, which($rowCorrMat_{i}$)$<$0)

// Mark problematic nodes and where they cause issues

end if

end for

6if  _max(correspondingColumns) = min(correspondingColumns) and
length(correspondingColumns)$\geq$ $|columnLabels=min(correspondingColumns)|$
_ then

rowWrongNodes = $\\{\\}$; queue(columnWrongNodes,
$\\{v:columnLabels(v)=min(correspondingColumns)\\}$ )

end if

// if all problematic row nodes are due to one column community, and there are
more marked row nodes than members of that column community, mark that whole
column community as problematic, and unmark all row nodes

7 for  _u in rowWrongNodes_ do

$\widehat{K_{r}}$++; rowLabels(u) = $\widehat{K_{r}}$ // move each marked node
to its own new singleton communities

end for

8L = calculateMeasureL(W, rowlabels, columnLabels)

9 return rowLabels, columnLabels, rowWrongNodes, rowColumnNodes, L

Algorithm 7 removeNodesA

Procedure  _removeNodesB(_W, rowLabels, columnLabels_)_

1 // $n_{r}$, $n_{c}$, $rowCorrMat$, and $colCorrMat$ are as in the previous
algorithm. Negative values may indicate a misclustering. This algorithm moves
problematic nodes to new communities based on their pattern of negative
correlations with other communities

2 rowWrongNodes = $\\{\\}$; columnWrongNodes= $\\{\\}$

3 correspondingColumns = $\\{\\}$; correspondingRows = $\\{\\}$; // Every top
level for loop or if statement done for row nodes will be done correspondingly
immediately thereafter for column nodes. From here on, for simplicity, only
the row node steps are shown.

4 for  _u in 1 to $n_{r}$_ do

if  _min( $rowCorrMat_{i}$)$<$0 or max($rowCorrMat_{i}$)= 0  or
$|rowLabels=rowLabels(u)|=2$_ then

queue(rowWrongNodes, u) // same step as in above algorithm

queue(correspondingColumns, which$(rowCorrMat_{i})<0$)

end if

end for

5if  _max(correspondingColumns) = min(correspondingColumns) and
length(correspondingColumns)$\geq$ $|columnLabels=min(correspondingColumns)|$
_ then

queue(columnWrongNodes, $\\{v:columnLabels(v)=min(correspondingColumns)\\}$ )
// if all problematic row nodes are due to one column community, and there are
more marked row nodes than members of that column community, mark that whole
column community as problematic

end if

additionalRowCommunities = 0

6for  _i in 1 to $\widehat{K_{r}}$_ do

binaryString = $\\{\\}$

for  _$\\{u:\\{rowLabels(u)=i\\}\cap rowWrongNodes\\}$_ do

binaryString(u) = binaryToInteger($rowCorrMat_{u}<0$) // nodes coming from
same community that have same pattern of negatives are moved to same new
community

end for

for  _val in unique(binaryString)_ do

additionalRowCommunities++

rowLabels($\\{$u: binaryString(u) = val$\\}$) =
$\widehat{K_{r}}$+additionalRowCommunities

end for

end for

7L = calculateMeasureL(W, rowlabels, columnLabels)

return rowLabels, columnLabels, rowWrongNodes, rowColumnNodes, L

Algorithm 8 removeNodesB

Algorithm 9 is meant only to be run on the result of Algorithm 7, seeking to
place those newly removed singleton communities into pre-existing communities.
Algorithm 10 is a bit more radical, in that it uses the same criteria to move
communities into other communities, but it will even combine pre-existing
communities. Furthermore, both algorithms have a toggle called banRemain,
which, if turned on, prevents nodes that have just been removed from their
previous community from rejoining that community. Algorithm 9 is run both with
and without banRemain on the reuslt of Algorithm 7, and Algorithm 8 is run
both with and without banRemain on the result of 8, as well as without
banRemain on the result of 7. Agglomeration is subsequently run on all these
results. All told, that will create 14 clustering options to be compared with
the current result.

Procedure  _regoupNodesA(_W, rowLabels, columnLabels, rowWrongNodes,
rowColumnNodes, banRemain, previousRowLabels, previousColumnLabels_)_

1 for  _i in max(previousRowLabels)+1 to max(rowLabels)_ do

2 for  _k in 1 to max(previousRowLabels)_ do

rowCommunityvalues(i, k) =
$\sum\limits_{j=1}^{\widehat{K_{c}}}max(0,n_{k}-2)\times corr\big{(}W(i,{v\in
j}),rowCommunityAggregate(k,{v\in j})\big{)}$

end for

end for

// repeat above for loop for columns

3 for  _i in max(previousRowLabels)+1 to max(rowLabels)_ do

if  _banRemain = 0_ then

rowLabels$(\\{u:rowLabels=i\\})=\underset{i^{\prime}\in\\{1,...,max(previousRowLabels)\\}}{\rm{argmax}}(rowCommunityvalues_{i})$

end if

else

rowLabels$(\\{u:rowLabels=i\\}\cap\\{u\not\in
rowWrongNodes\\})=\underset{i^{\prime}\in\\{1,...,max(previousRowLabels)\\}}{\rm{argmax}}(rowCommunityvalues_{i})$

rowLabels$(\\{u:rowLabels=i\\}\cap\\{u\in
rowWrongNodes\\})=\underset{i^{\prime}\in\\{1,...,max(previousRowLabels)\\}\backslash\\{previousNodeLabels(u)\\}}{\rm{argmax}}(rowCommunityvalues_{i})$

end if

end for

// repeat above for loop for columns

4 L = calculateMeasureL(W, rowlabels, columnLabels)

5 return rowLabels, columnLabels, L

Algorithm 9 regroupNodesA

Procedure  _regoupNodesB(_W, rowLabels, columnLabels, rowWrongNodes,
rowColumnNodes, banRemain, previousRowLabels, previousColumnLabels_)_

1 for  _i in 1 to max(rowLabels)_ do

2 for  _k in 1 to max(rowLabels)_ do

rowCommunityvalues(i, k) =
$\sum\limits_{j=1}^{\widehat{K_{c}}}max(0,n_{k}-2)\times
corr\big{(}rowCommunityAggregate(i,v\in j),rowCommunityAggregate(k,{v\in
j})\big{)}$

end for

end for

// repeat above for loop for columns

3 for  _i in 1 to max(rowLabels)_ do

if  _banRemain = 0_ then

rowLabels$(\\{u:rowLabels=i\\})=\underset{i^{\prime}\in\\{1,...,max(rowLabels)\\}}{\rm{argmax}}(rowCommunityvalues_{i})$

end if

else

rowLabels$(\\{u:rowLabels=i\\}\cap\\{u\not\in
rowWrongNodes\\})=\underset{i^{\prime}\in\\{1,...,max(rowLabels)\\}}{\rm{argmax}}(rowCommunityvalues_{i})$

rowLabels$(\\{u:rowLabels=i\\}\cap\\{u\in
rowWrongNodes\\})=\underset{i^{\prime}\in\\{1,...,max(rowLabels)\\}\backslash\\{previousNodeLabels(u)\\}}{\rm{argmax}}(rowCommunityvalues_{i})$

end if

end for

// repeat above for loop for columns

4 L = calculateMeasureL(W, rowlabels, columnLabels)

5 return rowLabels, columnLabels, L

Algorithm 10 regroupNodesB


# Emergence of extreme events in a quasi-periodic oscillator

Premraj Durairaj1, Sathiyadevi Kanagaraj1, Suresh Kumarasamy1, Karthikeyan
Rajagopal1,2 1Centre for Nonlinear Systems, Chennai Institute of Technology,
Chennai-600 069, Tamilnadu, India.
2Department of Electronics and Communications Engineering, University Centre
for
Research and Development, Chandigarh University, Mohali, 140 413, Punjab,
India.

(:)

###### Abstract

Extreme events are unusual and rare large-amplitude fluctuations that occur
can unexpectedly in nonlinear dynamical systems. Events above the extreme
event threshold of the probability distribution of a nonlinear process
characterize extreme events. Different mechanisms for the generation of
extreme events and their prediction measures have been reported in the
literature. Based on the properties of extreme events, such as rare in the
frequency of occurrence and extreme in amplitude, various studies have shown
that extreme events are both linear and nonlinear in nature. Interestingly, in
this work, we report on a special class of extreme events which are nonchaotic
and nonperiodic. These nonchaotic extreme events appear in between the quasi-
periodic and chaotic dynamics of the system. We report the existence of such
extreme events with various statistical measures and characterization
techniques.

Extreme events, Morse Oscillator, nonchaotic extreme events

###### pacs:

05.45.-a

Extreme events are unanticipated, rare events that occur in many natural and
engineering systems. Extreme events (EE) can exist in various forms, including
floods, cyclones, droughts, pandemics, power outages, material ruptures,
explosions, chemical contamination, and stock market crashes, among others
albeverio06 . Such events have a severe impact on real-world situations. Thus,
it is necessary to understand the relevant mechanism and its generic
characteristics for the occurrence of EE in order to prevent such EE. As a
result, the researchers focused on exploring the EE in diverse nonlinear
oscillators hohmann10 ; metzger14 ; mathis15 ; birkholz16 , maps prem1 , and
neural networks saha . Further, the extreme events have also been identified
in a super-fluid helium ganshin08 , plasmabailung11 , optical fibers solli07 ,
lasers zamora13 , and capillary wave Onorato  etc.

However, depending on the characteristics of a dynamical system, the
occurrence of EE has been discovered under a variety of mechanisms, including
internal crises, on-off intermittency, blowout bifurcations, stick-slip
bifurcations, and so on prem1 ; Venk ; Grebogi ; zamora13 ; Reinoso . For
instance, prior studies reveal that EE can arise as a result of the abrupt
expansion and destruction of chaotic attractors produced by internal or
external crises zamora13 ; Grebogi . Further, interior crises are found to be
a critical mechanism for the occurrence of EE, when the trajectory of chaotic
attractors reaches the stable manifold of a saddle or unstable periodic orbit,
which increases the size of the chaotic attractors. Such a sudden expansion of
the chaotic attractor may result in EE. In addition, Pomeau-Manneville
intermittency is identified as another mechanism for the existence of EE. Such
intermittency can occur when the periodic oscillations are interspersed by
chaotic bursts, which further results in very large amplitude events. EEs can
also exist through the following other mechanisms. The sliding bifurcation
near the discontinuous boundary can cause EE. The trajectory of the attractors
might hop between coexisting attractors due to noise in multi-stable systems,
which can cause unusually large events. This is referred to as noise-induced
intermittency. The trajectory of the attractors in coupled systems departs
from the synchronization manifold to the transverse direction of the manifold.
During such a transition, a synchronization error of dynamics can show zero or
nonzero and is referred to as on-off intermittency review .

Moreover, previous studies discovered that extreme or rare events can occur as
a result of chaotic or stochastic processes review . In particular, the
appearance of EE has been reported in micro-electromechanical cantilevers with
discontinuous boundaries and diode lasers with phase-conjugate feedback
Mercier ; Kumarasamy . By applying the harmonic pump modulation to the fiber
laser the emergence of Rogue waves has been identified Pisarchik1 ; Pisarchik2
; rw2 ; hohmann10 . The EE in stochastic transport on networks has been
demonstrated using multiple random walks on complex networks Santhanam1 ;
Santhanam2 . Now the interesting question is whether extreme events can be
induced by nonchaotic signals. In literature, a study has shown nonchaotic and
nonperiodic have been well studied in the name of strange nonchaotic dynamics,
which arises during the attractor transition from quasi-periodicity to chaos
Pikovsky . One can find the generation mechanisms of these strange nonchaotic
attractors in literature Pikovsky ; AP ; prem2 . The results in the present
work show that similar to the strange nonchaotic dynamics, the nonperiodic and
nonchaotic dynamics show large-amplitude extreme events. The present study
opens a new area of study where the nonchaotic nonlinear process can also lead
to extreme events and the same has not been found reported.

To show the nonchaotic extreme events, we consider the Morse Oscillator (MO)
which is used to describe the motion of diatomic molecules. Importantly, the
MO has made substantial contributions in the fields of classical, semi-
classical, and quantum mechanics Lie ; Knop . The MO was used for photo-
dissociation molecules without any damping. In the presence of driving and
damping, the MO was exploited for multi-photon excitation of molecules,
pumping the local mode of polyatomic molecules Zdravkovic . We consider the
quasi-periodically forced MO and its dynamical equation can be written as

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle y$
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle
fsin(\omega_{1}t)+gsin(\omega_{2}t)+e^{-2x}-e^{-x}-\gamma y$ (1)

where $x$, and $y$ are the state variables of the system and $\gamma$ is a
damping parameter. The amplitudes of the first and second force are
represented by $f$ and $g$ and the corresponding frequencies are denoted by
$\omega_{1}$ and $\omega_{2}$, respectively.

Figure 1: (a) Time evolution of $x_{n}$ for the nonchaotic dynamics with
forcing amplitudes as (a) $f=g=0.255$, and (b) $f=g=0.278$. The $x_{n}$ is the
$n^{th}$ local peaks of the variable $x$. The horizontal black dot-dashed and
red dashed lines are the critical threshold lines defining the extreme events
(refer to the text for the meaning of $N$ and $A$). (c) The probability
distribution function corresponds to the extreme events and (d) return
interval (R) (inter-event interval) with respect to the probability of
recurrence times (PR) of the EE for (b). The filled circles and solid lines in
(d) represent the numerical data and the corresponding power-law fit. We fixed
the other parameter values as $\gamma=0.35$, $\omega_{1}=0.3$, and
$\omega_{2}=(\frac{\sqrt{5}-1}{2}$).

To manifest the existence of extreme events, we first depicted the time
evolution of the $x$-variable in Fig. 1(a) and Fig. 1(b) by fixing the
amplitude of the first and second forcing as $f=g=0.255$ and $f=g=0.278$. We
observe from Fig. 1(a) that some of the oscillation(event) has larger
amplitudes, while the rest of them take lower amplitudes. To check the larger
amplitude oscillations satisfy the extreme events criteria defined in the
literature, we use the following relation:

$\displaystyle x_{{}_{EE}}=<x_{n}>+N\sigma_{x{{}_{n}}},$ (2)

where $x_{{}_{EE}}$ is the critical amplitude threshold and $N$ is a
multiplication factor. The mean and standard deviation of the variable $x$ is
represented by $<x_{n}>$ and $\sigma_{x{{}_{n}}}$, respectively. Here, the
$x_{n}$ (an event) are the local peaks of the variable $x$. An event or a
local peak can satisfy extreme event criteria if it has a value higher than
the critical threshold defined by Eq. (2) with $N\geq 4$. To confirm the
presence of EE, we plotted the critical threshold on the time series for $N=5$
and $N=4$ in Figs. 1(a) and 1(b). We used two different $N$ values depending
on the time series. Though the choice of $N$ is arbitrary, we set the minimum
$N$ value as 4 in the present study. We also find the critical value of
$N_{max}$ for a range of each $f$ value – the details will be discussed below.
In both cases, we can see that some of the large amplitude events cross the
threshold line, confirming the presence of EE. Since the choice of $N$ is
arbitrary in the previous criterion, we use another criterion defined by the
abnormality index; $A_{n}=\frac{Hf_{n}}{H_{1/3}}$ Mercier , where $Hf_{n}$ is
the difference between the maximum height of the event $n$ and the mean height
of its population, $Hf_{n}=x_{n}-\langle x_{n}\rangle_{n}$ and $H_{1/3}$ is
the average value among the highest one-third values of $Hf_{n}$. If an event
$x_{n}$ has abnormality $A_{n}$ greater than 2 then the event is termed an
extreme event. We find that both cases in Figs. 1(a) and 1(b) satisfy the
above criterion with abnormality index $A=3$ denoted by a dashed horizontal
line in the plots. It is evident that a few rare large amplitude events cross
the abnormality index line. We computed the probability distribution function
(PDF) in Fig. 1(c) for the time series shown in Fig. 1(b). The EE critical
threshold at $N=4$ is plotted as a vertical dashed line on the PDF diagram. In
the plot, the events with a finite probability above the critical threshold
line characterize the extreme events. We can plot similar probability
distribution for Fig. 1(a), however, for simplicity, we have plotted the PDF
corresponding to Fig. 1(b).

The above analysis shows that the observed behavior satisfies the extreme
events criterion in the amplitudes. Another important characteristic of
extreme events is an inter-event interval. The inter-event interval defines
the frequency occurrence of the events and should not have discrete values
(discrete values mean the periodic occurrence of events), rather it should
have a distribution over a range. In order to examine the distribution of
events in the observed time series, we find inter-event intervals (R) between
successive extreme events. Subsequently, we find the probability of such
inter-event intervals (PR) as shown in Fig. 1(d). Inter-event interval and its
probability obey power-law relations as given by
$log_{10}(PR)=a~{}log_{10}(R)^{b}$, where $a$ and $b$ are constants with
values $a=-0.006$ and $b=2.96$, respectively. The obtained numerical values
are depicted in a filled circle, and a continuous line shows the corresponding
power-law fit. The route for the emergence of EE and its transitions is
further estimated below using Lyapunov exponents (LE), amplitude maxima
$X_{max}$, critical factor $N_{max}$, and two-parameter analysis.

Figure 2:  (a) The two-parameter bifurcation diagram in $(f,g)$ space. Using
the range of Lyapunov exponents ($\lambda$) (denoted by the color bar) the
dynamical regions are marked. (b) The maximum Lyapunov exponents as a function
of forcing amplitude $f(=g)$, (c) maximum amplitude of the events $x_{max}$
(red) and the corresponding $N_{max}$ (Eq. 3) of the event (blue) by varying
the magnitude of $f(=g)$. The black line represents the extreme events
critical threshold ($x_{EE}$) drawn from Eq. (2) for N=4. The other parameter
values are fixed as the same as in Fig. 1.

To illustrate the global dynamical transition of the attractors and route of
the EE, the two-parameter diagram is drawn in $(f,g)$ space using the maximum
LE as shown in Fig. 2(a). The range of LE (shown in the color bar) denotes the
emergence of quasi-periodic, nonchaotic, and chaotic attractors in the
respective parameters of $f$ and $g$. If the forcing amplitudes $f$ and $g$
are small, attractors have a maximum negative LE, indicating the presence of a
quasi-periodic (QP) attractor region. To better comprehend QP attractors, we
plotted their time-evolution and phase portrait trajectories in Supplementary
Material Fig. S1 a(i, ii) for $f=g=0.23$, which show their bounded nature.
Thus, the EE critical threshold for this attractor is greater than the
amplitude of QP attractors. By increasing $f$ and $g$ values, the QP attractor
transits to a chaotic (CH) attractor via strange and nonchaotic dynamics in
which the LE takes the values from negative (near zero) to positive. To
distinguish between the strange nonchaotic and chaotic attractors, the time-
evolution and phase portrait trajectories are shown in Figs. S1 b(i,ii) and
Figs. S1 c(i,ii) in the supplementary materials by fixing $f=g=0.278$ and
$f=g=0.33$, respectively. Also, the frequency spectra can be used to
distinguish quasiperiodic, SNA and chaos. We have the frequency spectrum
analysis in the Supplementary material in Fig. S4 (a-c). When compared to the
chaotic attractor (which has a greater number of large amplitude
oscillations), we found the SNA shows fewer large amplitude oscillations. The
supplemental material’s Fig. S1 can be consulted for more information.
Furthermore, to show the dynamical transitions clearly, we displayed maximum
Lyapunov exponents in Fig. 2(b) by keeping the parameter ($f=g$) and varying
it along the diagonal dashed line shown in Fig. 2(a). In Fig. 2(b), the
maximum LE is illustrated as a function of forcing amplitudes $f$ and $g$
($f=g$) in the range $(0.23<f(=g)<0.32)$. We observe that when the forcing
amplitudes are minimum in the mentioned range, LE takes negative values,
indicating quasi-periodic dynamics. While increasing the parameter, the
transition of LE from negative to positive values indicates the dynamical
transition of quasi-periodic behavior to chaotic behavior. Furthermore, we
found that the negative values of LE near-zero exhibit strange nonchaotic
behavior; extreme events are seen in this region. The literature has shown
that the EEs occur under chaotic dynamics review  through distinct routes and
stochastic processes like stochastic transport on networks has been
demonstrated using multiple random walks on complex networks Santhanam1 ;
Santhanam2 . Among the various routes, the occurrence of EEs in nonchaotic
dynamics is new and it has not been reported to the best of our knowledge.

To validate the occurrence of EEs in the SNA region, we find the maximum
amplitude $x_{max}$, extreme event threshold $x_{EE}$, and maximum value of
$N$ ($N_{max}$) of a given time series. In Fig. 2(c), we have plotted the
above quantities by varying the magnitude of $f=g$. The plot explains the
regime of extreme events in the following way. During the non-extreme regime,
the critical threshold $x_{EE}$ is larger than the $x_{max}$. It means that
the threshold is larger than the large amplitude oscillations and does not
satisfy the extreme events criterion. While in the EE regime, the $x_{max}$ is
larger than the EE critical threshold $x_{EE}$ (shaded EE region). This
explains that extreme events have a larger amplitude than the extreme event
criterion. Note that the SNA regime in the parameter range $f\in{0.28}$ to
${0.2912}$ shows no extreme events. As we discussed above, we fixed $N=4$ as
an arbitrary constant from the literature Bonatto . However, the maximum value
of the $N$ can be determined by rewriting Eq. (2), as

$\displaystyle N_{max}=\frac{max(x_{EE})-\langle
x_{n}\rangle}{\sigma_{x_{n}}}.$ (3)

Figure 3: Basin of attraction for $f=g=0.278$. $QP1$, $QP2$, and $QP3$ are the
quasi-periodic attractor-1, quasi-periodic attractor-2, and quasi-periodic
attractor-3, respectively. SNA represents the strange nonchaotic attractor. We
fixed the other parameter values the same as in Fig. 1.

In the SNA region shown in Fig. 2(c), we found that the multiplication factor
taking values between $4\leq N_{max}\leq 5.611$ when the forcing amplitudes in
the range from $0.256$ to $0.28$ denoted by shaded transparent pattern. The
plot of $N_{max}$ shows that depending on the parameter choice, the arbitrary
value can be chosen N$\in\\{4,5.611\\}$. Thus above results satisfy all the
criteria proposed for the extreme events and justify the existence of EEs in
the SNA regime.

As we discussed earlier, the observed EEs are nonchaotic and nonperiodic. At
the same time, the parameters corresponding to the strange nonchaotic EEs show
multiple stable behaviors. The multi-stable behavior can be seen from the
basins of attraction drawn for a range of initial conditions. Figure 3 is
drawn by varying the initial states $x_{0}$ and $y_{0}$ of the system for the
parameters given in Fig. 1 caption. We can see that basin of nonchaotic and
nonperiodic behavior or SNA is embedded within the basin of quasi-periodic
dynamics. Outside the SNA basin, we have found three different basins which
contain quasi-periodic attractors. All the three different quasi-periodic
attractor basins and the SNA basin, denoted by QP1, QP2, QP3, and SNA
respectively in Fig. 3. In supplemental material Fig. (S2), each of the quasi-
periodic attractors is depicted. Figure 3 shows that extreme events occur for
specific values of initial conditions. The size of these basins changes as we
vary the parameter within the EEs regime marked in Fig. 2.

Figure 4: Two parameter phase diagram in $(f,g)$ space (plotted using Eq. 2
for fixed initial condition $(x_{0},y_{0})=(0.3,0.2)$), to distinguish the
existence of extreme events (EE) and non-extreme events (NEE), respectively.
We fixed the other parameter values as the same as in Fig. 1.

Similarly, to determine the regime of the extreme event in the parametric
space between $f$ and $g$, a two-parameter diagram is drawn as shown in Fig.
4. The white regime in the plot shows the extreme events for the combinations
of parameter $(f,g)$ separated with the help of Eq. 2 from the non-extreme
events (NEE– denoted by blue color). By comparing Fig. 2(a) with Fig. 4 we can
say that EEs occur in the SNA region (however some of the SNA parameter regime
may not contain EEs).

Figure 5: Singular continuous spectrum for fixing the forcing amplitudes $f,g$
$=$ $0.278$. (a) The logarithmic plot of $|x(\alpha,N)|^{2}$ against $N$. The
red and black lines denote the numerical values and the corresponding power-
law fit. (b) Fractal path in the complex plane of $x$. The other parameter
values are defined as $\gamma=0.35$, $\omega_{1}=0.3$,
$\omega_{2}=(\frac{\sqrt{5}-1}{2}$).

To show the generality of the existence of EEs in the SNA regime, we present
the regime of EEs for $\gamma=0.4$ in the supplementary material Figs. S3
(a),(b). This result validates the presence of strange nonchaotic extreme
events in the selected parameter regime. In the following section, we
characterize the observed behavior as strange and nonchaotic in nature. For
this purpose, we perform singular continuous spectrum analysis and
distribution of finite-time Lyapunov exponents.

To validate the strange nonchaotic dynamics, we plot singular continuous
spectrum Pikovsky  in Fig. 5 using partial Fourier sum of the signal $x$ given
by $X(\alpha,N)=\sum_{m=1}^{N}x_{m}e^{2\pi im\alpha}$, where $\alpha$ is
proportional to the external frequency ($\omega_{1}$) and $N$ is the length of
the time series. The red and black lines show the singular continuous spectrum
and the corresponding power-law fit. When N is considered as time,
$|X(\alpha,N)|^{2}$ grows with N, that is $|X(\alpha,N)|^{2}\sim N^{\beta}$,
where $\beta$ is the slope. When the signal possesses the properties of
strange nonchaotic dynamics, the corresponding slope values lie between
$1<\beta<2$. For this case, the slope value $\beta=1.576$ confirms the
existence of strange nonchaotic dynamics shown in Fig. 5(a). The corresponding
path of Brownian motion with fractal structure in complex [$Re(x),Im(x)$]
plane also confirms the strange nonchaotic dynamics in Fig. 5(b).

Figure 6: Finite time Lyapunov exponent with respect to probability
distribution function (PDF) for SNAs by fixing the three distinct finite time
periods $T=500$ (red line), $T=1000$ (blue dashed line), and $T=1500$ (black
dotted line) with $f,g=0.278$.

The strange nonchaotic dynamics are also validated using another statistical
characterization known as the distribution of finite-time Lyapunov exponents.
The distribution takes both positive and negative values, but the area under
the curve is maximum in the negative regime for strange nonchaotic dynamics.
Figure 6 plotted for three different finite time intervals $T=500,~{}1000$,
and $1500$, the distribution has a large negative region compared to the
positive region showing nonchaotic dynamics. From these analyses, the observed
dynamics are strange (nonperiodic) as well as nonchaotic, which also shows the
large amplitude and rare events.

The present letter shows a mechanism of the emergence of extreme events in a
quasi-periodically forced Morse oscillator. As a function of forcing
amplitude, we found the transition from quasi-periodic (QP) to chaotic (CH)
attractor via strange nonchaotic extreme events. During such extreme event
dynamics, we found a long excursion of trajectories that are away from the
bounded attractor, while the chaotic attractors show many higher amplitude
peaks. To confirm the existence of EEs, we estimated the critical threshold,
and it is observed that the higher amplitude peaks in the EE cross the
critical threshold while the peaks in the CH and QP attractor do not. The
dynamical transitions of the attractors and the occurrence of nonchaotic EE
dynamics are manifested through maximum Lyapunov exponents. The observed
extreme events are further validated using the probability distribution and
return interval (inter-event interval) with respect to the probability of
recurrence times of the EE. Extreme events are abnormal and unexpected events
that occur in many natural and man-made systems. Understanding the mechanism
or route can help to anticipate the onset of EEs. Early works on extreme
events show the chaotic nature of the extreme events because of the rare and
extreme amplitude properties of extreme events. The present study shows an
unknown emergence of extreme events that are nonchaotic and nonperiodic
extreme events. This finding shed light on the new direction where extreme
events can happen as a nonchaotic process.

We gratefully acknowledge this work is funded by the Center for Nonlinear
Systems, Chennai Institute of Technology (CIT), India, vide funding number
CIT/CNS/2022/RP-016.

## References

  * (1) S. Albeverio, V. Jentsch and H. Kantz, Extreme Events in Nature and Society, The Frontiers Collection (Springer, Berlin, 2006).
  * (2) R. Höhmann, U. Kuhl, H. J. Stöckmann, L. Kaplan and E. J. Heller, Phys. Rev. Lett., 104, (2010) 093901.
  * (3) J. J. Metzger, R. Fleischmann and T. Geisel, Phys. Rev. Lett., 112, (2014) 203903.
  * (4) A. Mathis, L. Froehly, S. Toenger, F. Dias, G. Genty and J. M. Dudley, Sci. Rep., 5, (2015) 12822.
  * (5) S. Birkholz, C. Brée, I. Veselić, A. Demircan and G. Steinmeyer, Sci. Rep., 6, (2016) 35207.
  * (6) D. Premraj, K. Suresh, S. A. Pawar, L. Kabiraj, A. Prasad, and R. I. Sujith, EPL, 134, (2021) 34006.
  * (7) A. Saha, and U. Feudel, Phys. Rev. E, 5, (2017) 062219.
  * (8) A. N. Ganshin, V. B. Efimov, G. V. Kolmakov, L. P. Mezhov-Deglin and P. V. E. McClintock, Phys. Rev. Lett., 101, (2008) 065303.
  * (9) H. Bailung, S. K. Sharma and Y. Nakamura, Phys. Rev. Lett., 107, (2011) 255005.
  * (10) D. R. Solli, C. Ropers, P. Koonath and B. Jalali, Nature, 450, (2007) 1054.
  * (11) J. Zamora-Munt, B. Garbin, S. Barland, M. Giudici, J. R. Rios Leite, C. Masoller and J. R. Tredicce, Phys. Rev. A, 87, (2013) 035802.
  * (12) M. Onorato, T. Waseda, A. Toffoli, L. Cavaleri, O. Gramstad, P. A. E. M. Janssen, T. Kinoshita, J. Monbaliu, N. Mori, A. R. Osborne and M. Serio, Phys. Rev. Lett. 102, (2009) 114502.
  * (13) S. C. Venkataramani, B. R. Hunt, E. Ott, D. J. Gauthier and J. C. Bienfang, Phys. Rev. Lett., 77, (1996) 5361.
  * (14) C. Grebogi, E. Ott, F. Romeiras and J. A. Yorke, Phys. Rev. A, 36, (1987) 5365.
  * (15) J. A. Reinoso, J. Zamora-Munt, and C. Masoller, Phys. Rev. E, 87, (2013) 062913.
  * (16) S. N. Chowdhury, A. Ray, S. K. Dana and D. Ghosh, Phys. Rep., 966, (2022) 1-52.
  * (17) E. Mercier, A. Even, E. Mirisola, D. Wolfersberger, and M. Sciamanna, Phys. Rev. E, 91, 042914 (2015).
  * (18) K. Suresh, and A. N. Pisarchik, Phys. Rev. E, 98, (2018) 032203.
  * (19) A. N. Pisarchik, R. Jaimes-Reátegui, R. Sevilla-Escoboza, G. Huerta-Cuellar, and M. Taki, Phys. Rev. Let., 107, 274101 (2011).
  * (20) A.N. Pisarchik, A.E. Hramov, Multistability in Physical and Living Systems: Characterization and Applications (Springer, 2022).
  * (21) F. T. Arecchi, U. Bortolozzo, A. Montina, and S. Residori, Phys. Rev. Lett., 106, 153901 (2011).
  * (22) A. S. Pikovsky, U. Feudel and S. P. Kuznetsov, Strange nonchaotic attractors: Dynamics between order and chaos in quasi periodically forced systems (World Scientific (2006)).
  * (23) V. Kishore, M. S. Santhanam, and R. E. Amritkar, Phys. Rev. let., 106, 188701 (2011).
  * (24) A. Kumar, S. Kulkarni and M. S. Santhanam, Chaos, 30, 043111, (2020).
  * (25) A. Prasad, S. S. Negi and R. Ramaswamy, Inter. J. Bifur. and Chaos, 11, (2001) 291-309.
  * (26) D. Premraj, K. Suresh, J. Palanivel and K. Thamilmaran, Communica. Non. Scie. Numer. Simula., 50, (2017) 103.; D. Premraj, K. Suresh, K. Thamilmaran and K. Rajagopal, The European Phys. J. Spe. Top., (2022) 1-7.
  * (27) G. C. Lie and J. M. Yuan, The J. chem. phys., 84, (1986) 5486.
  * (28) W. Knop and W. Lauterborn, The J. chem. phys., 93, (1990) 3950.
  * (29) S. Zdravkovic, A. N. Bugay and A. Y. Parkhomenko, Nonlinear Dyn., 90, (2017) 2841.
  * (30) C. Bonatto and A. Endler, Phys. Rev. E, 96, (2017) 012216.
  * (31) A. S. Pikovsky and U. Feudel, J. Phys. A, 27, (1994) 5209.