# Grass-roots optimization of coupled oscillator networks

Pranick R. Chamlagai  Department of Mathematics, Trinity College, Hartford, CT
06106, USA  Dane Taylor  Department of Mathematics, University at Buffalo,
State University of New York, Buffalo, NY 14260, USA  Per Sebastian Skardal
<EMAIL_ADDRESS>Department of Mathematics, Trinity College,
Hartford, CT 06106, USA

###### Abstract

Synchronization is critical for system function in applications ranging from
cardiac pacemakers to power grids. Existing optimization techniques rely
largely on global information, and while they induce certain local properties,
those alone do not yield optimal systems. Therefore, while useful for
designing man-made systems, existing theory provides limited insight into
self-optimization of naturally-occurring systems that rely on local
information and offer limited potential for decentralized optimization. Here
we present a method for “grass-roots” optimization of synchronization, which
is a multiscale mechanism involving local optimizations of smaller subsystems
that are coordinated to collectively optimize an entire system, and the
dynamics of such systems are particularly robust to islanding or targeted
attacks. In addition to shedding light on self-optimization in natural
systems, grass-roots optimization can also support the parallelizable and
scalable engineering of man-made systems.

The ability for large systems of dynamical units to self-organize and produce
robust collective behavior continues to drive a large body of research
Pikovsky2003 ; Strogatz2004 . Applications include cardiac dynamics
Bychkov2020JACC , brain dynamics Kopell2000PNAS , cell signaling
Prindle2012Nature , and power grids Rohden2012PRL . Weak synchronization and
desynchronization events often lead to pathological behavior, e.g., spiral
wave breakup in cardiac tissue Fenton2002Chaos ; Panfilov2007PNAS  and black
outs in power grids Dorfler2013PNAS , thereby motivating optimized systems for
strong, robust synchronization. While man-made systems such as power grids can
be designed and calibrated using global structural and dynamical information
Pecora1998PRL ; Nishikawa2006PRE , such information is likely unavailable to
naturally occurring systems. Notably, a great deal is known about how
biological systems function, however comparatively little is understood about
the self-optimization processes that are tasked with constructing and
maintaining/repairing such systems. Prominent examples include cardiac
pacemakers that initialize strongly synchronized pulses that propagate through
tissue Mangoni2008  and coordination of chromosomal activity through cell
differentiation Rajapakse2009PNAS . For synchronizing systems that rely on
collective behavior, it is reasonable to assume that the related optimization
mechanisms are themselves a collective, coordinated behavior. A stronger
theoretical understanding of such mechanisms for collective (self)
optimization will deepen our understanding of diverse types of biological (and
other) systems and has the potential to revolutionize the way we engineer
systems—or rather, design systems to engineer themselves. To this end,
collective optimizations constitute an under-explored family of collective
behavior, and there is a lack of multiscale optimization theory to provide
insight into how local optimizations might coordinate to optimize
globally–both in the context of synchronization and more broadly.

In this paper, we explore grass-roots optimization for coupled oscillator
networks, whereby the parallel optimization of smaller subsystems can be
coordinated to collectively optimize the global synchronization properties of
the entire system. In general, subsystems of a network can be defined in a
variety of ways: communities Girvan2002PNAS , spatially distinct regions in a
geometric network Barthelemy2011PhysRep , or other partitions of a network
after embedding in a more general metric space Coiffman2005PNAS . Our main
finding is an intuitive multiscale mechanism for grass-roots optimization of
synchronization that involves two steps: local subsystem optimization, whereby
subsystems are optimized in parallel; and global subsystem balancing, whereby
the subsystems are balanced with one another. We derive this mechanism from
first principles using the Synchrony Alignment Function (SAF) framework, which
provides an objective measure of a system’s synchronization properties and has
been used in a number of synchronization optimization tasks Skardal2014PRL ;
Skardal2016Chaos ; Taylor2016SIAM ; Skardal2017Chaos ; Skardal2019SIAM ;
Arola2021Chaos . We demonstrate the utility of grass-roots optimization across
a range of networks where structural subsystems arise naturally: networks with
community structure, a power grid, and noisy geometric networks that systems
with spatial constraints for connections, such as calcium release sites in
cardiac pacemaker cells Bychkov2020JACC  and self-coordinating chromosomes
Rajapakse2009PNAS . We show that the global synchronization properties of
grass-roots optimized systems are nearly identical to those of globally-
optimized systems, and importantly, these properties are also more robust to
subsystem dismantling, e.g., due to targeted attack or intentional
‘islanding’. Grass-roots optimization provides a viable mechanism by which
biological systems can robustly self-optimize and provides engineering
strategies that are decentralized, parallelizable, and scalable.

Figure 1: Grass-roots synchronization. Illustrations of (a) a random network
with two communities, (b) the IEEE RTS 96 power grid, and (c) a random
geometric network. (d)–(f) The degree of synchronization $r$ and (g)–(i)
synchronization error $1-r$ as a function of coupling strength $K$ for the
three respective network types with either randomly allocated frequencies
(green triangles), globally-optimized frequencies (blue circles), or grass-
roots optimized frequencies (red crosses).

We consider networks of coupled, heterogeneous phase oscillators whose
dynamics are given by

$\displaystyle\dot{\theta}_{i}=\omega_{i}+K\sum_{j=1}^{N}A_{ij}H(\theta_{j}-\theta_{i}),$
(1)

where $\theta_{i}$ and $\omega_{i}$ are the phase and natural frequency of
oscillator $i=1,\dots,N$, parameter $K$ is the global coupling strength,
network structure is encoded in an adjacency matrix $A$, and $H(\cdot)$ is a
$2\pi$-periodic coupling function. Here, we focus on the case of unweighted,
undirected networks with $A_{ij}=1$ if oscillators $i$ and $j$ are connected
and $0$ otherwise, although these properties may be relaxed without much
trouble. We also use classical Kuramoto coupling Kuramoto , given by
$H(\cdot)=\sin(\cdot)$, but emphasize that one may choose other functions $H$
provided that $H^{\prime}(0)>0$ and $H(\Delta\theta)=0$ for some
$\Delta\theta$ near zero. Notably, phase oscillator models such as Eq. (1)
have been found to be suitable models for naturally-occuring phenomena such as
chromosomal coordination Rajapakse2009PNAS  and integrate and fire dynamics of
cardiac pacemakers Politi2015PRE , as well as mechanical systems such as power
grids Porco2013 ; Skardal2015SciAdv . The degree of synchronization is
measured by the magnitude $r\in[0,1]$ of the Kuramoto order parameter
$re^{i\psi}=N^{-1}\sum_{j=1}^{N}e^{i\theta_{j}}$. By linearizing around the
synchronized state one obtains

$\displaystyle r\approx
1-\frac{J(\bm{\omega},L)}{2K^{2}},~{}\text{where}~{}J(\bm{\omega},L)=\frac{1}{N}\sum_{j=2}^{N}\frac{\langle\bm{v}^{j},\bm{\omega}\rangle^{2}}{\lambda_{j}^{2}}$
(2)

is the Synchrony Alignment Function (SAF) Skardal2014PRL . The SAF utilizes
the alignment of the natural frequencies $\bm{\omega}$ with the eigenvalues
$\\{\lambda^{j}\\}_{j=1}^{N}$ and eigenvectors $\\{\bm{v}^{j}\\}_{j=1}^{N}$ of
the combinatorial Laplacian, $L=D-A$, where $D=\text{diag}(k_{1},\dots,k_{N})$
is a diagonal matrix that encodes the nodal degrees,
$k_{i}=\sum_{j=1}^{N}A_{ij}$. Synchronization is optimized (i.e., $r$ is
maximized) by minimizing $J(\bm{\omega},L)$, which may be done by aligning
$\bm{\omega}$ with the eigenvectors of $L$ that are associated with larger
eigenvalues. Minimizing the SAF by setting $\omega=v^{N}$ also reveals
intuitive key properties of synchrony optimized systems including degree-
frequency correlations and anti-correlations between the frequencies of
neighboring oscillators Skardal2014PRL . While such local properties are
associated with synchronization, they alone do not guarantee it, nor do they
offer insight toward mesoscale/multiscale properties and mechanisms enabling
collective optimization.

We now present a method for grass-roots optimization of synchronization,
including a multiscale mechanism in which subsystems coordinate local
optimizations to optimize a system’s global synchronization properties. We
will present a detailed derivation later, and first summarize our main
findings and implications. We consider networks that can be partitioned into
$C$ subsystems such that the adjacency matrix $A$ may be rewritten in a block
form $A=A_{D}+B$, where $A_{D}=\text{diag}(A^{(1)},\dots,A^{(C)})$ is a block-
diagonal matrix containing the subsystems’ separate adjacency matrices, and
the off-diagonal blocks of $B$ encode edges between subsystems. We assume that
the blocks in $B$ are sparser than the diagonal blocks in $A_{D}$. For each
subsystem $s$, we define its associated combinatorial Laplacian matrix
$L^{(s)}$ and its associated vector $\bm{\omega}^{(s)}$ of frequencies. As we
will show below, under the condition where the subsystems’ mean oscillator
frequencies are equal, then the SAF for the full system may be approximated by
a linear combination of the subsystem-specific SAFs,

$\displaystyle
J(\bm{\omega},L)\approx\eta_{1}J(\bm{\omega}^{(1)},L^{(1)})+\cdots+\eta_{C}J(\bm{\omega}^{(C)},L^{(C)}),$
(3)

where $\eta_{s}$ is the fraction of nodes in subsystem $s$. This result leads
to the following multiscale mechanism for grass-roots optimization: (i)
_Global balancing of subsystems_ : achieve a balanced set of local mean
frequencies across all $C$ subsystems, i.e., minimize
$\text{max}_{s,s^{\prime}}|\langle\bm{\omega}^{(s)}\rangle-\langle\bm{\omega}^{(s^{\prime})}\rangle|$;
(ii) _Local optimization of subsystems_ : optimize the local SAFs, i.e.,
minimize $J(\bm{\omega}^{(s)},L^{(s)})$ for each $s$. This framework is
flexible and may be used under a wide range of application-specific
constraints. Notably, these two intuitive steps are a plausible mechanism that
can be utilized by biological (and other natural) systems to self-optimize
using local/global mechanisms, and it helps fill the theoretical gap between
existing (global) optimization theory and known (local) heuristic properties
that promote synchrony (e.g., degree-frequency correlations).

We now illustrate the effectiveness of grass-roots optimization across three
classes of networks: (i) networks with community structure (which are
generated by the stochastic block model Holland1983 , contain two communities
of sizes $N^{(1;2)}=100$, and have mean intra-degree $\langle
k^{(1;2)}\rangle=5$ and mean inter-degree $\langle k^{(12)}\rangle=1$); (ii)
the RTS 96 power grid Grigg1999IEEE ; (iii) and noisy geometric networks
Taylor2015NatComms  (which consist of $N=200$ nodes placed uniformly within a
$4\times 1$ box with $95\%$ of links placed between the closest possible nodes
pairs and the other $5\%$ of links placed randomly, resulting in a mean degree
of $\langle k\rangle=8$). As shown in Figs. 1(a)–(c), we partition the three
classes of networks into two, three, and four subsystems, respectively. (The
four subsystems of the geometric networks are defined by the $\pm$ sign
combinations in the first two non-trivial eigenvectors of $L$.)

For each network, we assume that natural frequencies are given and cannot be
modified, but may be rearranged. Thus, a global balance between subsystems
[step (i)] may be obtained by shuffling frequencies between subsystems, while
the subsystems may be locally optimized [step (ii)] by then shuffling
frequencies within each subsystem. To optimize each network, we use an accept-
reject algorithm, proposing $5\times 10^{4}$ switches between randomly chosen
pairs of frequencies and accepting switches that decrease the SAF. In Figs.
1(d)–(f), we plot the synchronization profile, $r$ vs $K$, for systems with
randomly allocated frequencies (green triangles), globally optimized
frequencies (blue circles) and grass-roots optimized frequencies (red crosses)
for the three classes of networks. All data points are averaged across $50$
random networks and natural frequency realizations (drawn from the standard
normal distribution) except for the power grid, where the same network is used
throughout. Note the comparably strong synchronization properties for both the
global and grass-roots optimized cases. To differentiate the two cases we plot
the synchronization error $1-r$ vs $K$ in a log-log scale in Figs. 1(g)–(i),
revealing that grass-roots optimization is very effective across a wide range
of network structures.

Figure 2: Robustness to islanding and target attacks. (a) Example of local
(subsystem) order parameters for the RTS 96 power grid before ($t<0$) and
after ($t\geq 0$) islanding for global (solid blue) and grass-roots (dashed
red) optimization. (b) Density of local (subsystem) SAFs after islanding for
global (solid blue) and grass-roots (dashed red) optimization. (c)
Illustration of the islanded subsystems in the RTS 96 power grid.

Next we show that grass-roots optimized networks outperform globally optimized
networks when subsystems are islanded from one another or otherwise dismantled
by a targeted attack. For instance, modern power grids feature microgrids that
are smaller subsystems that may be separated, i.e., “islanded”, from the
larger grid Porco2013 . We predict such a feature to be advantageous to
synchronizing biological processes, which is a main motivator for our work. As
an example, we consider the RTS 96 power grid before and after the islanding
of three subsystems [as indicated in Fig. 2(c)]. In Fig. 1(a) we plot time
series of three local order parameters for system designed using global (solid
blue) and grass-root (dashed red) optimization. We use $K=1$ and normally-
distributed frequencies. Edges between subsystems are removed at time $t=0$.
Before islanding ($t<0$) both cases display strong synchronization properties.
After islanding ($t\geq 0$) the globally-optimized system displays
significantly weaker synchronization properties and a desynchronization event
(indicated by oscillations). On the other hand, the grass-roots optimized
system maintains its strong synchronization properties. This is further
demonstrated in Fig. 2(b), where we plot the density of local, i.e.,
subsystem-specific, SAFs for globally (solid blue) and grass-roots (dashed
red) optimized systems obtained from $10^{4}$ realizations. We indicate the
respective means $\overline{J(\bm{\omega},L)}=0.1427$ and $0.0629$ of the
local SAFs for the globally and grass-roots optimized cases with vertical
lines.

We conclude by finally presenting our local approximation of the SAF, which
has allowed us to identify the multiscale mechanism (i.e., steps i-ii)
underlying grass-roots optimization. For simplicity, we first consider the
case of two subsystems, leaving further generalization to the Supplemental
Material (SM). Writing the adjacency matrix as
$A=\begin{bmatrix}A^{(1)}&B^{(12)}\\\ B^{(12)T}&A^{(2)}\end{bmatrix}$, where
$A^{(1)}\in\mathbb{R}^{N_{1}\times N_{1}}$, $A^{(2)}\in\mathbb{R}^{N_{2}\times
N_{2}}$, $B^{(12)}\in\mathbb{R}^{N_{1}\times N_{2}}$, and $N_{1}$ and $N_{2}$
are the sizes of the respective subsystems, the Laplacian is given by
$L=L_{0}+L_{B}$, where $L_{0}=\begin{bmatrix}L^{(1)}&0\\\
0&L^{(2)}\end{bmatrix}$, $L_{B}=\begin{bmatrix}D_{B^{(12)}}&-B^{(12)}\\\
-B^{(12)T}&D_{B^{(12)T}}\end{bmatrix}$, and $L^{(1,2)}=D^{(1,2)}-A^{(1,2)}$
with diagonal matrices $D_{B^{(12)}}$ and $D_{B^{(12)T}}$ whose entries
correspond to row sums of $B^{(12)}$ and $B^{(12)T}$, respectively. We assume
$B^{(12)}$ to be sparser than $A^{(1)}$ and $A^{(2)}$ so that
$\left\|L_{B}\right\|\ll\left\|L_{0}\right\|$ under a suitable matrix norm
(e.g., the Frobenius norm). We then define $\Delta
L=(\|L_{0}\|/\|L_{B}\|)L_{B}$ so that $L(\epsilon)=L_{0}+\epsilon\Delta L$
recovers the original network structure for the choice
$\epsilon=\|L_{B}\|/\|L_{0}\|\ll 1$.

Next, we discuss the spectral properties of $L_{0}$. Since this matrix encodes
the two subsystems in isolation, its eigenvalue spectrum is the union of the
eigenvalue spectrum of $L^{(1)}$ and $L^{(2)}$. Specifically, ordering the
eigenvalues of $L^{(1)}$ and $L^{(2)}$, respectively,
$0=\mu_{1}<\mu_{2}\leq\cdots\leq\mu_{N_{1}}$ and
$0=\nu_{1}<\nu_{2}\leq\cdots\leq\nu_{N_{2}}$ (where we assume that the
subsystems are themselves connected), this implies that $L_{0}$ has two zero
eigenvalues, $\lambda_{1}=\lambda_{2}=0$, and the rest are positive. Since $0$
is a repeated eigenvalue of $L_{0}$, its nullspace requires some care. Rather
than choosing eigenvectors $\bm{v}^{1}\propto[\bm{1},\bm{0}]^{T}$ and
$\bm{v}^{2}\propto[\bm{0},\bm{1}]^{T}$, whose entries are constant within one
subsystem and zero within the other, it is advantageous to instead choose
$\bm{v}^{1}=\frac{1}{\sqrt{N}}[\bm{1},\bm{1}]^{T}$ and
$\bm{v}^{2}=\frac{\sqrt{N_{1}N_{2}}}{N}[\bm{1}/N_{1},-\bm{1}/N_{2}]^{T}$ so
that $\bm{v}^{1}$ is independent of $\epsilon$ and characterizes the nullspace
of $L(\epsilon)$, and $\bm{v}^{2}$ is associated with an eigenvalue that
converges to 0 as $\epsilon\to 0$ but is strictly positive for $\epsilon>0$.
The other $N-2$ eigenvectors of $L_{0}$ are given by
$\\{\bm{v}^{j}\\}_{j=3}^{N}=\left\\{[\bm{u}^{j},\bm{0}]^{T}\right\\}_{j=2}^{N_{1}}\bigcup\left\\{[\bm{0},\bm{x}^{j}]^{T}\right\\}_{j=2}^{N_{2}}$,
where $\\{\bm{u}^{j}\\}_{j=1}^{N_{1}}$ and $\\{\bm{x}^{j}\\}_{j=1}^{N_{2}}$
are the eigenvectors of $L^{(1)}$ and $L^{(2)}$.

Considering $0<\epsilon\ll 1$, each eigenvalue of $L(\epsilon)$ varies
continuously with $\epsilon$ Kato2013 , so we may write
$\lambda_{j}(\epsilon)=\lambda_{j}+\epsilon\delta\lambda_{j}^{(1)}+\epsilon^{2}\delta\lambda_{j}^{(2)}+\mathcal{O}(\epsilon^{3})$.
We similarly assume
$\bm{v}^{j}(\epsilon)=\bm{v}^{j}+\epsilon\delta\bm{v}^{j(1)}+\epsilon^{2}\delta\bm{v}^{j(2)}+\mathcal{O}(\epsilon^{3})$.
Since $\lambda_{2}(\epsilon)\ll 1$ and $\lambda_{j}(\epsilon)\sim 1$ for
$j=3,\dots,N$, the term associated with $j=2$ needs to be treated separately
from the others, and

$\displaystyle
J(\bm{\omega},L(\epsilon))=\frac{1}{N}\left(\frac{\langle\bm{\omega},\bm{v}^{2}(\epsilon)\rangle}{\lambda_{2}(\epsilon)}\right)^{2}+\frac{1}{N}\sum_{j=3}^{N}\left(\frac{\langle\bm{\omega},\bm{v}^{j}(\epsilon)\rangle}{\lambda_{j}(\epsilon)}\right)^{2}.$
(4)

Upon expanding the $N-1$ terms contributing to the SAF in Eq. (4), we find
that they all take a similar form except for a factor of $\epsilon$,

$\displaystyle\left(\frac{\langle\bm{\omega},\bm{v}^{j}(\epsilon)\rangle}{\lambda_{j}(\epsilon)}\right)^{2}=\epsilon^{\alpha_{j}}\left(\frac{\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{2}}\right)+\epsilon^{1+\alpha_{j}}\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{2}}-\frac{2\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{3}}\right)$
$\displaystyle+\epsilon^{2+\alpha_{j}}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle^{2}+2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(2)}\rangle}{(\lambda_{j})^{2}}-\frac{4\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{3}}+\frac{(3(\delta\lambda_{j}^{(1)})^{2}-2\lambda_{j}\delta\lambda_{j}^{(2)})\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{4}}\right)+\mathcal{O}(\epsilon^{3+\alpha_{j}}),$
(5)

where $\alpha_{j}$ is a term that equals -2 when $j=2$ and 0 when $j\geq 3$.
Due the the different scaling with $\epsilon$, the terms associated with $j=2$
are much larger than those for $j\geq 3$. Inserting Eq. (5) into Eq. (4)
yields

$\displaystyle J(\bm{\omega},$ $\displaystyle
L(\epsilon))=N^{-1}\epsilon^{-2}\left(\frac{\langle\bm{\omega},\bm{v}^{2}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{2}}\right)+N^{-1}\epsilon^{-1}\left(\frac{2\langle\bm{\omega},\bm{v}^{2}\rangle\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle}{(\delta\lambda_{2}^{(1)})^{2}}-\frac{2\delta\lambda_{2}^{(2)}\langle\bm{\omega},\bm{v}^{2}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{3}}\right)$
$\displaystyle+N^{-1}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle^{2}+2\langle\bm{\omega},\bm{v}^{2}\rangle\langle\bm{\omega},\delta\bm{v}^{2(2)}\rangle}{(\delta\lambda_{2}^{(1)})^{2}}-\frac{4\delta\lambda_{2}^{(2)}\langle\bm{\omega},\bm{v}^{2}\rangle\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle}{(\delta\lambda_{2}^{(1)})^{3}}+\frac{(3(\delta\lambda_{2}^{(2)})^{2}-2\delta\lambda_{2}^{(1)}\delta\lambda_{2}^{(3)})\langle\bm{\omega},\bm{v}^{2}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{4}}\right)$
$\displaystyle+\eta_{1}J(\bm{\omega}^{1},L_{1})+\eta_{2}J(\bm{\omega}^{2},L_{2})+\epsilon\left[N^{-1}\sum_{j=3}^{N}\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{2}}-\frac{2\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{3}}\right)\right]+\mathcal{O}(N^{-1}\epsilon,\epsilon^{2}),$
(6)

where we have used that
$\frac{1}{N}\sum_{j=3}^{N}\frac{\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{\lambda_{j}^{2}}=\eta_{1}J(\bm{\omega}^{(1)},L^{(1)})+\eta_{2}J(\bm{\omega}^{(2)},L^{(2)})$.
While Eq. (6) may appear daunting, the key insight is the presence of an inner
product $\langle\bm{\omega},\bm{v}^{2}\rangle$ in several leading-order terms.
Recalling the structure of $\bm{v}^{2}$, and writing
$\bm{\omega}=[\bm{\omega}^{(1)},\bm{\omega}^{(2)}]^{T}$, where
$\bm{\omega}^{(1)}$ and $\bm{\omega}^{(2)}$ are the frequency vectors
corresponding to the two subsystems, we have that
$\langle\bm{\omega},\bm{v}^{2}\rangle=\sqrt{\eta_{1}\eta_{2}}(\langle\bm{\omega}^{(1)}\rangle-\langle\bm{\omega}^{(2)}\rangle)$.
Thus, if the subsystems’ mean frequencies can be engineered to match,
$\langle\bm{\omega}^{(1)}\rangle=\langle\bm{\omega}^{(2)}\rangle$, then many
terms vanish to yield

$\displaystyle J$
$\displaystyle(\bm{\omega},L(\epsilon))=\eta_{1}J(\bm{\omega}^{1},L_{1})+\eta_{2}J(\bm{\omega}^{2},L_{2})$
$\displaystyle+\epsilon\left[N^{-1}\sum_{j=3}^{N}\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{2}}-\frac{2\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{3}}\right)\right]$
$\displaystyle+N^{-1}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{2}}\right)+\mathcal{O}(N^{-1}\epsilon,\epsilon^{2}),$
(7)

which recovers Eq. (3) to leading order for the case of $C=2$ subsystems. See
the SM for further generalization insights.

While recent progress has been made in optimizing collective behavior in
complex systems, the resulting techniques and methodologies rely largely on
global network information. These approaches express certain local properties
such as correlations between nodal degrees and natural frequencies
Skardal2014PRL ; Skardal2016Chaos , however such properties alone do not
optimize systems. This leaves open the critical question of how naturally-
occurring systems tune their own structure and dynamics to self-optimize, and
it is reasonable to consider that the optimization itself is a collective
behavior.

Grass-roots optimization is a multiscale mechanism for coordinating and
optimizing the local synchronization properties of a network’s subsystems and
is a plausible mechanism for collective (self) optimization within naturally-
occurring systems that have access to only local information, such as cardiac
pacemakers Bychkov2020JACC  and genetic oscillators Rajapakse2009PNAS . It can
also support the design of decentralized, parallelizable and scalable
algorithms to engineer man-made systems that are robust to network
dismantling. Notably, these very same features may have provided an
evolutionary advantage for biological systems that crucially depend on
synchronization.

###### Acknowledgements.

PRC was supported by the Interdisciplinary Science Program and Summer Research
Program at Trinity College. PSS was supported by NSF grant MCB-2126177. DT was
supported by NSF grant DMS-2052720 and Simons Foundation award #578333.

## References

  * [1] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: a universal concept in nonlinear sciences (Cambridge University Press, 2003).
  * [2] S. H. Strogatz, Sync: The emerging science of spontaneous order (Penguin UK, 2004).
  * [3] R. Bychkov, M. Juhaszova, K. Tsutsui, C. Coletta, M. D. Stern, V. A. Maltsev, and E. G. Lakatta, Synchronized Cardiac Impulses Emerge From Heterogeneous Local Calcium Signals Within and Among Cells of Pacemaker Tissue, JACC Clin. Electrophysiol. 6, 907 (2020).
  * [4] N. Kopell, G. B. Ermentrout, M. A. Whittington, and R. D. Traub, Gamma rhythms and beta rhythms have different synchronization properties, Proc. Natl. Acad. Sci. U.S.A. 97, 1867 (2000).
  * [5] A. Prindle, P. Samayoa, I. Razinkov, T. Danino, L. S. Tsimring, and J. Hasty, A sensing array of radically couples genertic ‘biopixels’, Nature 481, 39 (2012).
  * [6] M. Rohden, A. Sorge, M. Timme, and D. Witthaut, Self-organized synchronization in decentralized power grids, Phys. Rev. Lett. 109, 064101 (2012).
  * [7] F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans, Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity, Chaos 12, 852 (2002).
  * [8] A. V. Panfilov, R. H. Keldermann, and M. P. Nash, Drift and breakup of spiral waves in reaction–diffusion–mechanics systems, Proc. Natl. Acad. Sci. U.S.A 104, 7922 (2007)
  * [9] F. Dörfler, M. Chertkov, and F. Bullo, Synchronization in complex oscillator networks and smart grids, Proc. Natl. Acad. Sci. U.S.A. 110, 1005 (2013).
  * [10] M. E. Mangoni and J. Nargeot, Genesis and regulation of the heart automaticity, Physiol. Rev. 88, 919 (2008).
  * [11] I. Rajapakse, M. D. Perlman, D. Scalzo, C. Kooperberg, M. Groudine, and S. T. Kosak, The emergence of lineage-specific chromosomal topologies from coordinate gene regulation, Proc. Natl. Acad. Sci. U.S.A. 106, 6679 (2009).
  * [12] L. M. Pecora and T. L. Carroll, Master stability function for synchronized coupled systems, Phys. Rev. Lett. 80, 2109 (1998).
  * [13] T. Nishikawa and A. Motter, Synchronization is optimal in nondiagonalizable networks, Phys. Rev. E 73, 065106(R) (2006).
  * [14] M. Girvan and M. E. J. Newman, Community structure in social and biological networks, Proc. Natl. Acad. Sci U.S.A. 99, 7821 (2002).
  * [15] M. Barthélemy, Spatial Networks, Phys. Rep. 499, 1 (2011).
  * [16] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker, Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps, Proc. Natl. Acad. Sci. U.S.A. 102, 7426 (2005).
  * [17] P. S. Skardal, D. Taylor, and J. Sun, Optimal synchronization of complex networks, Phys. Rev. Lett. 113, 144101 (2014).
  * [18] P. S. Skardal, D. Taylor, and J. Sun, Optimal synchronization of directed complex networks, Chaos 26, 094807 (2016).
  * [19] D. Taylor, P. S. Skardal, and J. Sun, Synchronization of heterogeneous oscillators under network modifications: Perturbation and optimization of the synchrony alignment function, SIAM J. Appl. Math. 76, 1984 (2016).
  * [20] L. Arola-Fernández, P. S. Skardal, and A. Arenas, Geometric unfolding of synchronization dynamics on networks, Chaos 31, 061105 (2021).
  * [21] P. S. Skardal, R. Sevilla-Escoboza, V. Vera-Ávila, and J. M. Buldú, Optimal phase synchronization in networks of phase-coherent chaotic oscillators, Chaos 27, 013111 (2017).
  * [22] P. S. Skardal, D. Taylor, and J. Sun, Synchronization of network-coupled oscillators with uncertain dynamics, SIAM J. Appl. Math. 79, 2409 (2019).
  * [23] Y. Kuramoto, Chemical oscillations, waves, and turbulence (Springer, 2012).
  * [24] A. Politi and M. Rosenblum, Equivalence of phase-oscillator and integrate-and-fire models, Phys. Rev. E 91, 042916 (2015).
  * [25] J. W. Simpson-Porco, F. Dörfler, and F. Bullo, Synchronization and power sharing for droop-controlled inverters in islanded microgrids, Automatica 49, 2603 (2013).
  * [26] P. S. Skardal and A. Arenas, Control of coupled oscillator networks with application to microgrid technologies, Sci. Adv. 1, e1500339 (2015).
  * [27] P. W. Holland, K. B. Laskey, and S. Leinhardt, Stochastic blockmodels: First steps, Soc. Networks 5, 109 (1983).
  * [28] C. Grigg et al., The IEEE Reliability Test System—1996. A report prepared by the Reliability Test System Task Force of the Application of Probability Methods Subcommittee, IEEE Trans. Power Syst. 14, 1010 (1999).
  * [29] D. Taylor, F. Klimm, H. A. Harrington, M. Kramar, K. Mischaikow, M. A. Porter, and P.J. Mucha, Topological data analysis of contagion maps for examining spreading processes on networks, Nat. Commun. 6, 7723 (2015).
  * [30] T. Kato, Perturbation theory for linear operators, vol. 132 (Springer Science & Business Media, 2013).

Supplemental Material: Grass-roots optimization of coupled oscillator networks

In this Supplemental Material we generalize the local approximation of the SAF
for more than two subsystems. First we present the full derivation of the
approximation for the case of three subsystems, and then we discuss the
generalization to an arbitrary number of subsystems.

## Local Approximation of the SAF for Networks with Three Subsystems

To provide insight into systems with more than two subsystems, we present here
the case of three subsystems and derive a local approximation to the SAF
analogous to the one which we presented in the main text. In this case the
network adjacency matrix can be written in block form as

$\displaystyle A=\begin{bmatrix}A^{(1)}&B^{(12)}&B^{(13)}\\\
B^{(12)T}&A^{(2)}&B^{(23)}\\\ B^{(13)T}&B^{(23)T}&A^{(3)}\end{bmatrix},$ (1)

where $A^{(1)}$, $A^{(2)}$, and $A^{(3)}$ are the adjacency matrices for the
three subsystems and $B^{(12)}$, $B^{(13)}$, and $B^{(23)}$ captures the
connections between the respective subsystems. We denote the sizes of the
three subsystems by $N_{1}$, $N_{2}$, and $N_{3}$ so that
$A^{(1)}\in\mathbb{R}^{N_{1}\times N_{1}}$, $A^{(2)}\in\mathbb{R}^{N_{2}\times
N_{2}}$, $A^{(3)}\in\mathbb{R}^{N_{3}\times N_{3}}$,
$B^{(12)}\in\mathbb{R}^{N_{1}\times N_{2}}$,
$B^{(13)}\in\mathbb{R}^{N_{1}\times N_{3}}$, and
$B^{(23)}\in\mathbb{R}^{N_{2}\times N_{3}}$. We are interested then in the
perturbed combinatorial Laplacian, given by

$\displaystyle L(\epsilon)=L_{0}+\epsilon\Delta L,$ (2)

where

$\displaystyle L_{0}=\begin{bmatrix}L^{(1)}&0&0\\\ 0&L^{(2)}&0\\\
0&0&L^{(3)}\end{bmatrix},$ (3)

$\Delta L=(\|L_{0}\|/\|L_{B}\|)L_{B}$, and

$\displaystyle
L_{B}=\begin{bmatrix}D_{B^{(12)}+B^{(13)}}&-B^{(12)}&-B^{(13)}\\\
-B^{(12)T}&D_{B^{(12)T}+B^{(23)}}&-B^{(23)}\\\
-B^{(13)T}&-B^{(23)T}&D_{B^{(13)T}+B^{(23)T}}\end{bmatrix}.$ (4)

Once again, the choice $\epsilon=\|L_{B}\|/\|L_{0}\|\ll 1$ recovers the
original Laplacian matrix.

As in the two-subsystem case, it is useful to first discuss the spectral
properties of $L_{0}$. Since it is a block-diagonal matrix, its eigenvalues
are given by the union of the eigenvalues of the respective blocks,

$\displaystyle\\{\lambda_{j}\\}_{j=1}^{N}=\\{\mu_{j}\\}_{j=1}^{N_{1}}\bigcup\\{\nu_{j}\\}_{j=1}^{N_{2}}\bigcup\\{\eta_{j}\\}_{j=1}^{N_{3}},$
(5)

where $\\{\mu_{j}\\}_{j=1}^{N_{1}}$ denotes the eigenvalues of $L^{(1)}$,
$\\{\nu_{j}\\}_{j=1}^{N_{2}}$ denotes the eigenvalues of $L^{(2)}$, and
$\\{\eta_{j}\\}_{j=1}^{N_{3}}$ denotes the eigenvalues of $L^{(3)}$. The
associated eigenvectors are given by

$\displaystyle\\{\bm{v}^{j}\\}_{j=1}^{N}=\left\\{\begin{bmatrix}\bm{u}^{j}\\\
\bm{0}\\\
\bm{0}\end{bmatrix}\right\\}_{j=1}^{N_{1}}\bigcup\left\\{\begin{bmatrix}\bm{0}\\\
\bm{x}^{j}\\\
\bm{0}\end{bmatrix}\right\\}_{j=1}^{N_{2}}\bigcup\left\\{\begin{bmatrix}\bm{0}\\\
\bm{0}\\\ \bm{y}^{j}\end{bmatrix}\right\\}_{j=1}^{N_{3}}.$ (6)

where $\\{\bm{u}^{j}\\}_{j=1}^{N_{1}}$, $\\{\bm{x}^{j}\\}_{j=1}^{N_{2}}$, and
$\\{\bm{y}^{j}\\}_{j=1}^{N_{3}}$ are the associated eigenvectors for
$L^{(1)}$, $L^{(2)}$, and $L^{(3)}$, respectively. The most critical
observation to make is that each diagonal block of $L_{0}$ has a trivial
eigenvalue, namely, $\mu_{1},\nu_{1},\eta_{1}=0$, so the nullspace of $L_{0}$
is three-dimensional since it has a triple eigenvalue degeneracy at
$\lambda_{1,2,3}=0$. It is then convenient to rewrite the basis vectors for
this trivial eigenspace using the following eigenvectors:

$\displaystyle\bm{v}^{1}=\frac{1}{\sqrt{N}}\begin{bmatrix}\bm{1}\\\ \bm{1}\\\
\bm{1}\end{bmatrix},~{}~{}~{}\bm{v}^{2}=\frac{\sqrt{N_{1}N_{2}}}{N_{1}+N_{2}}\begin{bmatrix}\bm{1}/N_{1}\\\
-\bm{1}/N_{2}\\\
\bm{0}\end{bmatrix},~{}~{}~{}\bm{v}^{3}=\frac{\sqrt{N_{2}N_{3}}}{N_{2}+N_{3}}\begin{bmatrix}\bm{0}\\\
\bm{1}/N_{2}\\\ -\bm{1}/N_{3}\end{bmatrix},$ (7)

where, similar to the two subsystem case, $\bm{v}^{1}$ is the constant-valued
eigenvector that is associated with the synchronization manifold and whose
eigenvalue $\lambda_{1}=0$ remains constant as $\epsilon$ increases (i.e.,
$v^{1}(\epsilon)=v^{1}$ regardless of $\epsilon$). On the other hand,
$\bm{v}^{2}$ and $\bm{v}^{3}$ will play important roles in the perturbation
analysis since $\lambda_{2}(\epsilon)$ and $\lambda_{3}(\epsilon)$ must take
positive values for any $\epsilon>0$. We note that the vector
$\sqrt{N_{1}N_{3}}/(N_{1}+N_{3})\begin{bmatrix}\bm{1}/N_{1}\\\ \bm{0}\\\
-\bm{1}^{T}/N_{3}\end{bmatrix}$ may also be used in place of either
$\bm{v}^{2}$ or $\bm{v}^{3}$, but as it is just a linear combination of the
two vectors already chosen, it yields the same results given below.

Given the initial spectral properties of $L_{0}$, we consider the following
perturbative expansions. Specifically, for the eigenvalues of $L(\epsilon)$ we
have

$\displaystyle\lambda_{j}(\epsilon)$
$\displaystyle=\epsilon\delta\lambda_{j}^{(1)}+\epsilon^{2}\delta\lambda_{j}^{(2)}+\mathcal{O}(\epsilon^{3}),$
(8)

for $j=2,3$ and

$\displaystyle\lambda_{j}(\epsilon)$
$\displaystyle=\lambda_{j}+\epsilon\delta\lambda_{j}^{(1)}+\epsilon^{2}\delta\lambda_{j}^{(2)}+\mathcal{O}(\epsilon^{3}),$
(9)

for $j=4,\dots,N$. We again assume that the eigenvectors of $L(\epsilon)$ are
continuously differentiable to approximate

$\displaystyle\bm{v}^{j}(\epsilon)$
$\displaystyle=\bm{v}^{j}+\epsilon\delta\bm{v}^{j(1)}+\epsilon^{2}\delta\bm{v}^{j(2)}+\mathcal{O}(\epsilon^{3}).$
(10)

for $j=2,\dots,N$.

Our primary interest is the SAF of the perturbed network, and as we did in the
two subsystem case with the term associated with $j=2$, here we will treat the
terms associated with $j=2$ and $3$ separately:

$\displaystyle
J(\bm{\omega},L(\epsilon))=\frac{1}{N}\left(\frac{\langle\bm{\omega},\bm{v}^{2}(\epsilon)\rangle}{\lambda_{2}(\epsilon)}\right)^{2}+\frac{1}{N}\left(\frac{\langle\bm{\omega},\bm{v}^{3}(\epsilon)\rangle}{\lambda_{3}(\epsilon)}\right)^{2}+\frac{1}{N}\sum_{j=4}^{N}\left(\frac{\langle\bm{\omega},\bm{v}^{j}(\epsilon)\rangle}{\lambda_{j}(\epsilon)}\right)^{2}.$
(11)

We now consider the contribution of these different terms. Beginning with the
terms associated with $j=2$ and $3$, insert Eqs. (8) and (10) into the
relevant terms in Eq. (11), expand, and collect similar terms to obtain

$\displaystyle\frac{1}{N}\left(\frac{\langle\bm{\omega},\bm{v}^{j}(\epsilon)\rangle}{\lambda_{j}(\epsilon)}\right)^{2}$
$\displaystyle=N^{-1}\epsilon^{-2}\left(\frac{\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\delta\lambda_{j}^{(1)})^{2}}\right)+N^{-1}\epsilon^{-1}\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\delta\lambda_{j}^{(1)})^{2}}-\frac{2\delta\lambda_{j}^{(2)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\delta\lambda_{j}^{(1)})^{3}}\right)$
$\displaystyle+N^{-1}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle^{2}+2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(2)}\rangle}{(\delta\lambda_{j}^{(1)})^{2}}-\frac{4\delta\lambda_{j}^{(2)}\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\delta\lambda_{j}^{(1)})^{3}}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\frac{(3(\delta\lambda_{j}^{(2)})^{2}-2\delta\lambda_{j}^{(1)}\delta\lambda_{j}^{(3)})\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\delta\lambda_{j}^{(1)})^{4}}\right)+\mathcal{O}(N^{-1}\epsilon).$
(12)

On the other hand, for $j=4,\dots,N$, we insert Eqs. (9) and (10) into the
relevant terms in Eq. (11), expand, and collect similar terms to obtain

$\displaystyle\frac{1}{N}\left(\frac{\langle\bm{\omega},\bm{v}^{j}(\epsilon)\rangle}{\lambda_{j}(\epsilon)}\right)^{2}$
$\displaystyle=N^{-1}\left(\frac{\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{2}}\right)+N^{-1}\epsilon\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{2}}-\frac{2\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{3}}\right)$
$\displaystyle+N^{-1}\epsilon^{2}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle^{2}+2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(2)}\rangle}{(\lambda_{j})^{2}}-\frac{4\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{3}}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\frac{(3(\delta\lambda_{j}^{(1)})^{2}-2\lambda_{j}\delta\lambda_{j}^{(2)})\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{4}}\right)+\mathcal{O}(N^{-1}\epsilon^{3}).$
(13)

Inserting Eqs. (12) and (13) into Eq. (11), we then obtain

$\displaystyle J($
$\displaystyle\bm{\omega},L(\epsilon))=N^{-1}\epsilon^{-2}\left(\frac{\langle\bm{\omega},\bm{v}^{2}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{2}}+\frac{\langle\bm{\omega},\bm{v}^{3}\rangle^{2}}{(\delta\lambda_{3}^{(1)})^{2}}\right)$
$\displaystyle+N^{-1}\epsilon^{-1}\left(\frac{2\langle\bm{\omega},\bm{v}^{2}\rangle\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle}{(\delta\lambda_{2}^{(1)})^{2}}-\frac{2\delta\lambda_{2}^{(2)}\langle\bm{\omega},\bm{v}^{2}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{3}}+\frac{2\langle\bm{\omega},\bm{v}^{3}\rangle\langle\bm{\omega},\delta\bm{v}^{3(1)}\rangle}{(\delta\lambda_{3}^{(1)})^{2}}-\frac{2\delta\lambda_{3}^{(2)}\langle\bm{\omega},\bm{v}^{3}\rangle^{2}}{(\delta\lambda_{3}^{(1)})^{3}}\right)$
$\displaystyle+N^{-1}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle^{2}+2\langle\bm{\omega},\bm{v}^{2}\rangle\langle\bm{\omega},\delta\bm{v}^{2(2)}\rangle}{(\delta\lambda_{2}^{(1)})^{2}}-\frac{4\delta\lambda_{2}^{(2)}\langle\bm{\omega},\bm{v}^{2}\rangle\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle}{(\delta\lambda_{2}^{(1)})^{3}}+\frac{(3(\delta\lambda_{2}^{(2)})^{2}-2\delta\lambda_{2}^{(1)}\delta\lambda_{2}^{(3)})\langle\bm{\omega},\bm{v}^{2}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{4}}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left.\frac{\langle\bm{\omega},\delta\bm{v}^{3(1)}\rangle^{2}+2\langle\bm{\omega},\bm{v}^{3}\rangle\langle\bm{\omega},\delta\bm{v}^{3(2)}\rangle}{(\delta\lambda_{3}^{(1)})^{2}}-\frac{4\delta\lambda_{3}^{(2)}\langle\bm{\omega},\bm{v}^{3}\rangle\langle\bm{\omega},\delta\bm{v}^{3(1)}\rangle}{(\delta\lambda_{3}^{(1)})^{3}}+\frac{(3(\delta\lambda_{3}^{(2)})^{2}-2\delta\lambda_{3}^{(1)}\delta\lambda_{3}^{(3)})\langle\bm{\omega},\bm{v}^{3}\rangle^{2}}{(\delta\lambda_{3}^{(1)})^{4}}\right)$
$\displaystyle+\eta_{1}J(\bm{\omega}^{1},L_{1})+\eta_{2}J(\bm{\omega}^{2},L_{2})+\eta_{3}J(\bm{\omega}^{2},L_{3})+\epsilon\left[N^{-1}\sum_{j=4}^{N}\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{2}}-\frac{2\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{3}}\right)\right]+\mathcal{O}(N^{-1}\epsilon,\epsilon^{2}),$
(14)

where we have used that, for the three subsystem case, we have

$\displaystyle\frac{1}{N}\sum_{j=4}^{N}\frac{\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{\lambda_{j}^{2}}=\eta_{1}J(\bm{\omega}^{1},L_{1})+\eta_{2}J(\bm{\omega}^{2},L_{2})+\eta_{3}J(\bm{\omega}^{3},L_{3}).$
(15)

Lastly, to complete the analysis we consider not only the contributions of
$\langle\bm{\omega},\bm{v}^{2}\rangle$, but also
$\langle\bm{\omega},\bm{v}^{3}\rangle$. In particular, we note that

$\displaystyle\langle\bm{\omega},\bm{v}^{2}\rangle=\frac{\sqrt{\eta_{1}\eta_{2}}}{\eta_{12}}(\langle\omega^{1}\rangle-\langle\omega^{2}\rangle),$
(16)

and

$\displaystyle\langle\bm{\omega},\bm{v}^{3}\rangle=\frac{\sqrt{\eta_{2}\eta_{3}}}{\eta_{23}}(\langle\omega^{2}\rangle-\langle\omega^{3}\rangle),$
(17)

where $\eta_{ij}=(N_{i}+N_{j})/N$. Thus, if we may engineer the network such
that
$\langle\omega^{1}\rangle=\langle\omega^{2}\rangle=\langle\omega^{3}\rangle$,
then all terms in Eq. (14) with $\langle\bm{\omega},\bm{v}^{2}\rangle$ or
$\langle\bm{\omega},\bm{v}^{3}\rangle$ vanish, yielding

$\displaystyle J(\bm{\omega},L(\epsilon))$
$\displaystyle=\eta_{1}J(\bm{\omega}^{1},L_{1})+\eta_{2}J(\bm{\omega}^{2},L_{2})+\eta_{3}J(\bm{\omega}^{2},L_{3})+N^{-1}\left(\frac{\langle\bm{\omega},\delta\bm{v}^{2(1)}\rangle^{2}}{(\delta\lambda_{2}^{(1)})^{2}}+\frac{\langle\bm{\omega},\delta\bm{v}^{3(1)}\rangle^{2}}{(\delta\lambda_{3}^{(1)})^{2}}\right)$
$\displaystyle+\epsilon\left[N^{-1}\sum_{j=4}^{N}\left(\frac{2\langle\bm{\omega},\bm{v}^{j}\rangle\langle\bm{\omega},\delta\bm{v}^{j(1)}\rangle}{(\lambda_{j})^{2}}-\frac{2\delta\lambda_{j}^{(1)}\langle\bm{\omega},\bm{v}^{j}\rangle^{2}}{(\lambda_{j})^{3}}\right)\right]+\mathcal{O}(N^{-1}\epsilon,\epsilon^{2}),$
(18)

where the leading-order behavior of the perturbed SAF is simply given by a
weighted average of the subsystem-specific SAFs and the weights come from
their relative sizes, which is our desired result and the analogous version of
Eq. (7) in the main text.

## Local Approximation of the SAF for Networks with an Arbitrary Number of
Subsystems

Before concluding, we emphasize that the three subsystem case above informs
the generalization of the local approximation to an arbitrary number of
subsystems. In particular, for $C$ subsystems, the unperturbed Laplacian
$L_{0}$ will contain $C$ diagonal blocks, each with a trivial eigenvalue.
Thus, a basis for the trivial eigenspace must be chosen so that, in addition
to $\bm{v}^{1}\propto\bm{1}$, there are $C-1$ eigenvectors whose eigenvalues
will becomes positive for positive $\epsilon$. This can be done by choosing,
for instance,

$\displaystyle\bm{v}^{2}=\begin{bmatrix}\bm{1}/N_{1}\\\ -\bm{1}/N_{2}\\\
\bm{0}\\\ \vdots\\\
\bm{0}\end{bmatrix},~{}~{}\bm{v}^{3}=\begin{bmatrix}\bm{0}\\\ \bm{1}/N_{2}\\\
-\bm{1}/N_{3}\\\ \vdots\\\
\bm{0}\end{bmatrix},~{}~{}\cdots~{}~{},~{}~{}\bm{v}^{j}=\begin{bmatrix}\vdots\\\
\bm{1}/N_{j-1}\\\ -\bm{1}/N_{j}\\\ \vdots\\\
\bm{0}\end{bmatrix},~{}~{}\cdots~{}~{},~{}~{}\bm{v}^{C}=\begin{bmatrix}\bm{0}\\\
\vdots\\\ \bm{0}\\\ \bm{1}/N_{C-1}\\\ -\bm{1}/N_{C}\end{bmatrix}.$ (19)

Then, after expansion, setting
$\langle\omega^{1}\rangle=\cdots=\langle\omega^{C}\rangle$ causes the two
lowest order contributions to $J(\bm{\omega},L(\epsilon))$ originating from
the terms associated with $j=2,\dots,C$ to vanish, yielding, to leading order,

$\displaystyle
J(\bm{\omega},L(\epsilon))\approx\eta_{1}J(\bm{\omega}^{1},L^{(1)})+\cdots+\eta_{C}J(\bm{\omega}^{C},L^{(C)}).$
(20)


# KnowMAN: Weakly Supervised Multinomial Adversarial Networks

Luisa März ⋄,†, Ehsaneddin Asgari †, Fabienne Braune †,
Franziska Zimmermann† and Benjamin Roth ⋄
⋄ Digital Philology, Research Group Data Mining and Machine Learning,
University of Vienna, Austria
† NLP Expert Center, Data:Lab, Volkswagen AG, Munich, Germany

###### Abstract

The absence of labeled data for training neural models is often addressed by
leveraging knowledge about the specific task, resulting in heuristic but noisy
labels. The knowledge is captured in labeling functions, which detect certain
regularities or patterns in the training samples and annotate corresponding
labels for training. This process of weakly supervised training may result in
an over-reliance on the signals captured by the labeling functions and hinder
models to exploit other signals or to generalize well. We propose KnowMAN, an
adversarial scheme that enables to control influence of signals associated
with specific labeling functions. KnowMAN forces the network to learn
representations that are invariant to those signals and to pick up other
signals that are more generally associated with an output label. KnowMAN
strongly improves results compared to direct weakly supervised learning with a
pre-trained transformer language model and a feature-based baseline.

##  1 Introduction

Neural approaches rely on labeled data sets for training. For many tasks and
languages, such data is either scarce or not available at all. Knowledge-based
weak supervision tackles this problem by employing _labeling functions (LFs)_.
LFs are manually specified properties, e.g. keywords, that trigger the
automatic annotation of a specific label. However, these annotations contain
noise and biases that need to be handled.

A recent approach for denoising weakly supervised data is Snorkel (Ratner et
al., 2020). Snorkel focuses on estimating the reliability of LFs and of the
resulting heuristic _labels_. However, Snorkel does not address biases on the
_input side_ of weakly supervised data, which might lead to learned
representations that overfit the characteristics of specific LFs, hindering
generalization. We address the problem of overfitting to the LFs in this
paper.

Other approaches tackle such overfitting by deleting the LF signal completely
from the input side of an annotated sample: For example, Go et al. (2009)
strip out emoticons that were used for labeling the sentiment in tweets, and
Alt et al. (2019) mask the entities used for distant supervision of relation
extraction training data Mintz et al. (2009). However, as LFs are often
constructed from the most prototypical and reliable signals (e.g., keywords),
deleting them entirely from the feature space might – while preventing over-
reliance on them – hurt prediction quality considerably. However, we find a
way to blur the signals of the LFs instead of removing them.

In this work we propose KnowMAN (Knowledge-based Weakly Supervised Multinomial
Adversarial Networks), a method for controllable _soft deletion_ of LF
signals, allowing a trade-off between reliance and generalization. Inspired by
adversarial learning for domain adaptation Chen and Cardie (2018a); Ganin and
Lempitsky (2015), we consider LFs as domains and aim to learn a LF-invariant
feature extractor in our model. KnowMAN is composed of three modules: a
feature extractor, a classifier, and a discriminator. Specifically, KnowMAN
employs a classifier that learns the actual task and an adversarial opponent,
the LF- discriminator, that learns to distinguish between the different LFs.
Upstream of both is the shared feature extractor to which the gradient of the
classifier and the reversed gradient of the discriminator are propagated. In
our experiments, the feature extractor for encoding the input is a multi-layer
perceptron on top of either a bag-of-words vector or a transformer
architecture, but KnowMAN is in principle usable with any differentiable
feature extractor.

KnowMAN consistently outperforms our baselines by 2 to 30% depending on the
dataset. By setting a hyperparameter $\lambda$ that controls the influence of
the adversarial part we can control the degree of discarding the information
of LF-specific signals. The optimal $\lambda$ value depends on the dataset and
its properties.

The contributions of this work are i) proposing an adversarial architecture
for controlling the influence of signals associated with specific LFs, ii)
consistent improvements over weakly supervised baselines, iii) release of our
code 111https://github.com/LuisaMaerz/KnowMAN. To our knowledge, we are the
first that apply adversarial learning to overcome the noisiness of labels in
weak supervision.

##  2 Method

Figure 1: KnowMAN architecture. The figure depicts one iteration over a batch
of inputs. The parameters of $\mathcal{C}$ and $\mathcal{F}_{s}$ are updated
together, following the green arrows. The LF discriminator $\mathcal{D}$ is
updated following the red arrows. Solid lines indicate forward, dashed lines
the backward pass.

Our approach is composed of three interacting modules i) the shared feature
extractor $\mathcal{F}_{s}$, ii) the classifier $\mathcal{C}$ and iii) the LF
discriminator $\mathcal{D}$. The loss function of $\mathcal{C}$ rewards the
classifier $\mathcal{C}$ for predicting the correct label for the instance,
and the gradient is used for optimizing the shared feature extractor and
classifier modules towards that goal. At the same time, the loss function for
the LF-discriminator $\mathcal{D}$ rewards predicting which LF was responsible
for labeling an instance. However, in adversarial optimization, KnowMAN
backpropagates the _reversed_ gradient for the LF-discriminator, hence the
information indicative for distinguishing between specific LFs is weakened
throughout the network. The hyperparameter $\lambda$ is used to control the
level of weakening the signals - the higher we choose the value the more
influence is assigned to the discriminator information that goes into
$\mathcal{D}$. The result of the interplay between classifier and LF-
discriminator is a shared feature representation that is good at predicting
the labels while reducing the influence of LF-specific signals, encouraging
the shared feature extractor to take other information (correlated with all
LFs for a class) into account.

In Figure 1, the arrows illustrate the training flow of the three modules. Due
to the adversarial nature of the LF discriminator $\mathcal{D}$, it has to be
trained with a separate optimizer (red arrows), while the rest of the network
is updated with the main optimizer (green arrows). When $\mathcal{D}$ is
trained the parameters of $\mathcal{C}$ and $\mathcal{F}_{s}$ are frozen and
vice versa.

To calculate the losses we utilize canonical negative log-likelihood loss
(NLL) and use it for both, the classifier and the LF discriminator. The
classification NLL can be formalized as:

$\mathcal{L}_{C}(\hat{y_{i}},y_{i})=-\log P(\hat{y_{i}}=y_{i})$ (1)

where $y_{i}$ is the (weakly supervised) annotated label and $\hat{y_{i}}$ is
the prediction of the classifier module $\mathcal{C}$, for a training sample
$i$. Analogously, we can define the NLL for the LF discriminator:

$\mathcal{L_{D}}(\hat{lf}_{i},lf_{i})=-\log P(\hat{lf}_{i}=lf_{i})$ (2)

where $lf_{i}$ is the actual LF used for annotating sample $i$ and
$\hat{lf}_{i}$ is the predicted LF by the discriminator $\mathcal{D}$.
Accordingly, we minimize two different objectives within KnowMAN:

$J_{\mathcal{C}}=\sum_{i=1}^{N}\mathcal{L_{C}}(\mathcal{C}(\mathcal{F}_{s}(x_{i});y_{i}))$
(3)
$J_{\mathcal{D}}=\sum_{i=1}^{N}\mathcal{L_{D}}(\mathcal{D}(\mathcal{F}_{s}(x_{i});lf_{i}))$
(4)

Here the shared feature extractor has two different objectives: i) help
$\mathcal{C}$ to achieve better classification performance and ii) make the
feature distribution invariant to the signals from the LFs. This is captured
by the shared objective:

$J_{\mathcal{F}_{s}}=J_{\mathcal{C}}+\lambda\cdot(-J_{\mathcal{D}})$ (5)

where $\lambda$ is the parameter that controls the adversarial influence i.e.
the degree of LF signal blur. $-J_{\mathcal{D}}$ is the reversed loss of the
LF discriminator $\mathcal{D}$ that represents $\mathcal{C}s$ adversarial
opponent. In general, the exact implementation or architecture of the
individual modules is interchangeable and can be set up as required. This
makes KnowMAN a universally applicable and easily customizable architecture.

##  3 Experiments

###  3.1 Data

For our experiments we use three standard datasets for weak supervision.

Spam. Based on the YouTube comments dataset Alberto et al. (2015) there is a
smaller Spam dataset from Snorkel Ratner et al. (2020) where the task is to
classify if a text is relevant to a certain YouTube video or contains spam.
This dataset is very small and does consist of a train and a test set only.
For the $10$ LFs keywords and regular expressions are used.

Spouse. This dataset for extracting the _spouse_ relation has also been
created by Snorkel, it is based on the Signal Media One-Million News Articles
Dataset Corney et al. (2016). The $9$ LFs use information from a knowledge
base, keywords and patterns. One peculiarity of this dataset is that over 90%
of the instances do not hold a spouse relation.

IMDb. The IMDb dataset contains movie reviews that should be classified in
terms of their sentiment (binary, positive or negative sentiment). The LFs
used for this dataset are occurrences of positive and negative keywords from
Hu and Liu (2004). A particular characteristic of this data set is the large
amount of $6800$ LFs, which constitutes a particular challenge to the Snorkel
denoising framework. As a result Snorkel fails to calculate its generative
model, since its memory consumption exceeds the available limit of 32GB RAM.

###  3.2 Experimental setup

For the experiments we use two different methods for encoding the input: i)
TF-IDF encoding and ii) a DistilBERT transformer. For TF-IDF encoding, we
vectorize222https://scikit-
learn.org/stable/modules/generated/sklearn.feature_extraction.text.TfidfVectorizer.html
the input sentences and feed them to a simple MLP. In the transformer setting,
the sequences of words are encoded using a pretrained DistilBERT. Similar to
BERT Devlin et al. (2019), DistilBERT is a masked transformer language model,
which is a smaller, lighter, and faster version leveraging knowledge
distillation while retaining 97% of BERT’s language understanding capabilities
Sanh et al. (2019).

Our encoder takes the representation of the CLS token from a frozen DistilBERT
and learns a non-linear transformation with a drop-out layer to avoid
overfitting Srivastava et al. (2014):

$h_{i}=DistilBERT(Sentence_{i})_{[CLS]}$ ${F_{s}}_{i}=Dropout(ReLU(f(h_{i})))$

where $DistilBERT(.)_{[CLS]}$ generates the hidden state of the BERT’s
classifier token (CLS) and the function $f$ represents a linear transformation
for the $i^{th}$ sentence.

The classifier and discriminator networks following the feature extractor are
in line with the implementation of Chen and Cardie (2018a) for domain-
adversarial learning. Both are simple sequential models with dropout, batch
normalization, $ReLU$ activation and softmax as the last layer. Please see our
code for implementation details. In the TF-IDF setup we use Adam Kingma and Ba
(2014) for both optimizers. When using transformer encoding the $\mathcal{D}$
optimizer again is Adam and the $\mathcal{C}$ optimizer is AdamW Loshchilov
and Hutter (2018), as this yielded more stable results.

Baselines For each input encoding we implemented several baselines. Weakly
supervised TF-IDF (WS TF-IDF) and Weakly supervised DistilBERT (WS
DistilBERT). Both calculate the labels for each instance in the train set
based on their matching LFs. WS TF-IDF directly applies a logistic regression
classifier to the input and the calculated labels. WS DistilBERT directly uses
the DistilBERT uncased model for English Sanh et al. (2019) as a prediction
model. The second baseline (Feature TF-IDF, Feature DistilBERT) uses feature
extractor and classifier layers of KnowMAN without taking the information of
$\mathcal{D}$ into account (this is equal to setting $\lambda$ to zero). We
also fine-tuned the pure language model (Fine-tuned DistilBERT) without
further transformations and without integrating the KnowMAN architecture.

We also compare with training TF-IDF and DistilBERT models on labels denoised
by Snorke (Snorkel TF-IDF, Snorkel DistilBERT). However, Snorkel denoising
failed for the IMDb data set due to the large amount of LFs.

| Spam |  | Spouse |  | IMDb
---|---|---|---|---|---
| Acc | P | R | F1 | Acc
WS TF-IDF | 0.87 | 0.12 | 0.83 | 0.20* | 0.65*
Feature TF-IDF | 0.91 | 0.12 | 0.76 | 0.21* | 0.75*
Snorkel TF-IDF | 0.81 | 0.18 | 0.63 | 0.28* | 0.50*
KnowMAN TF-IDF | 0.94 | 0.16 | 0.72 | 0.35 | 0.77
Fine-tuned DistilBERT | 0.92 | 0.14 | 0.78 | 0.24 | 0.70
WS DistilBERT | 0.87 | 0.09 | 0.90 | 0.17* | 0.67*
Feature DistilBERT | 0.86 | 0.18 | 0.80 | 0.29* | 0.74
Snorkel DistilBERT | 0.88 | 0.13 | 0.70 | 0.23* | 0.49*
KnowMAN DistilBERT | 0.90 | 0.27 | 0.67 | 0.39 | 0.76
Table 1: Results on the test sets. The * indicates that KnowMAN performs
significantly better than the marked model. For the Spouse data set we do
report significance for the F1 scores only.

KnowMAN We refer to the KnowMAN architecture as TF-IDF KnowMAN and DistilBERT
KnowMAN. Depending on the dataset we choose different $\lambda$ values. We
also implemented two ways of evaluation and best model saving during training:
i) evaluate after each batch and save the best model, ii) evaluate after a
certain number of steps in between the batches and save the best model.

Hyperparameters We perform hyperparameter tuning using Bayesian optimization
(Snoek et al., 2012) for the IMDb and Spouse datasets. For Spam,
hyperparameters are not optimized, as no validation set is available. Sampling
history and resulting hyperparameters are reported in the Appendix, Figures 2,
3 as well as hyperparameters chosen for the Spam data set.

Evaluation For the evaluation of the IMDb and the Spam datasets we use
accuracy, for the Spouse dataset we use the macro F1 score of the positive
class. To check statistical significance we use randomized testing Yeh (2000).
Results are considered significant if $\rho$ < 0.05.

###  3.3 Results

The results of the experiments are shown in Table 1. For the TF-IDF setup
KnowMAN TF-IDF outperforms the baselines across all datasets. We find the
optimal $\lambda$ values as follows: Spam/Spouse/IMDb = 2/5/4.9. Using the
additional feature extractor layer (Feature TF-IDF) is beneficial compared to
direct logistic regression for all datasets. Snorkel TF-IDF can outperform the
other two baselines for the Spouse dataset only.

Fine tuning of DistilBERT can not outperform our best KnowMAN. However, for
the Spam dataset Fine-tuned DistilBERT gives better results than KnowMAN
DistilBERT but still is worse than KnowMAN TF-IDF. Using WS DistilBERT gives
the same results for the Spam dataset and slightly better results for IMDb,
when compared to WS TF-IDF, for Spouse the performance decreases. Snorkel
DistilBERT can outperform the other two baselines for the Spam dataset only.
The low performance of Snorkel on IMDb (for both DistilBERT and TF-IDF) might
be explained by the very large amount of LF for this dataset. The KnowMAN
DistilBERT results across datasets are in line with the TF-IDF setup - KnowMAN
can outperform all baselines for the Spouse and IMDb dataset. We observe that
$\lambda=5$ for Spouse and $\lambda=1$ for IMDb is most beneficial when using
DistilBERT. For the Spam dataset we observe that KnowMAN (with $\lambda=2$)
outperforms all the baselines, except for the fine-tuned DistilBERT model.

Discussion The performance drop we observe with DistilBERT for KnowMAN
compared to the tf-idf setup of the IMDb dataset could be explained by
implementation details. Due to memory issues we have to truncate the input
when using DistilBERT. Since the movie reviews from IMDb are rather long this
could harm performance. Since the Spam dataset is very small a single wrongly
classified instance can have great impact on the results. This could explain
why KnowMAN TF-IDF outperforms KnowMAN DistilBERT here as well. In general we
could not perform hyperparameter optimization for the DistilBERT experiments
due to memory issues. Therefore the results for that experiments might not
have reached their optimum. However, the results show the value of using
KnowMAN though. Overall our results confirm the assumption that KnowMAN
enables a focus shift of the shared feature extractor from the signals of the
LFs towards signals of other valuable information. KnowMAN consistently
improves over the other experiments significantly - except for the Spam
dataset. We assume that the dataset size is too small to see significant
changes in the results. Compared to the implementation of Chen and Cardie
(2018a) we could not use the specialized domain feature extractor for our
datasets in the experiments. This is due to the fact that our test sets do not
contain information about LF matches. However, we will address this issue by
integrating a mixture of experts module for the specialized feature extractor
as recommended by Chen et al. (2019).

##  4 Related Work

Adversarial neural networks have been used to reduce the divergence between
distributions, such as Goodfellow et al. (2014), Chen et al. (2018) and Ganin
and Lempitsky (2015). The latter proposed an architecture for gradient
reversal and a shared feature extractor. Unlike us, they focused on a binary
domain discriminator. Similarly, Chen and Cardie (2018a) use an adversarial
approach in a multinomial scenario for domain adaptation.

Some works on adversarial learning in the context of weak supervision focus on
different aspects and only share similarity in name with our approach: Wu et
al. (2017) use _virtual adversarial training_ Miyato et al. (2017) for
perturbing input representations, which can be viewed as a general
regularization technique not specific to weakly supervised learning. Qin et
al. (2018); Zeng et al. (2018) use generative adversarial mechanisms for
selecting _negative_ training instances that are difficult to discriminate
from heuristically annotated ones for a classifier.

Several approaches have focused on denoising the labels for weakly supervised
learning Takamatsu et al. (2012); Manning et al. (2014); Lin et al. (2016).
Snorkel Ratner et al. (2020) is one of the most general approaches in this
line of work. However, Snorkel only models biases and correlations of LFs, and
does not consider problems of weak supervision that may stem from biases in
the features and learned representations.

A recent approach that focuses on denoising weakly supervised data is Sedova
et al. (2021). Knodle is a framework for comparison of different methods that
improve weakly supervised learning. We use some of their datasets for our
approach but denoise the signals of the LFs during training.

##  5 Conclusion

We propose KnowMAN - an adversarial neural network for training models with
noisy weakly supervised data. By integrating a shared feature extractor that
learns labeling function invariant features, KnowMAN can improve results on
weakly supervised data drastically across all experiments and datasets in our
setup. The experiments also show that the adverse effect of labeling function-
specific signals is highly dependent on the datasets and their properties.
Therefore, it is crucial to fine-tune the $\lambda$ parameter on a validation
set to find the optimal degree of blurring the labeling function signals.
Since the modules in the KnowMAN architecture are easily exchangeable, KnowMAN
can be applied to any architecture and dataset labeled with heuristic labeling
functions.

## Acknowledgements

This research was funded by the WWTF through theproject ”Knowledge-infused
Deep Learning for Nat-ural Language Processing” (WWTF Vienna ResearchGroup
VRG19-008), by the Deutsche Forschungs-gemeinschaft (DFG, German Research
Foundation) -RO 5127/2-1.

## References

  * Alberto et al. (2015) Tulio Alberto, Johannes Lochter, and Tiago Almeida. 2015.  Tubespam: Comment spam filtering on youtube.  pages 138–143.
  * Alt et al. (2019) Christoph Alt, Marc Hübner, and Leonhard Hennig. 2019.  Improving relation extraction by pre-trained language representations.  In _Automated Knowledge Base Construction (AKBC)_.
  * Chen et al. (2019) Xilun Chen, Ahmed Hassan Awadallah, Hany Hassan, Wei Wang, and Claire Cardie. 2019\.  Multi-source cross-lingual model transfer: Learning what to share.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 3098–3112, Florence, Italy. Association for Computational Linguistics.
  * Chen and Cardie (2018a) Xilun Chen and Claire Cardie. 2018a.  Multinomial adversarial networks for multi-domain text classification.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 1226–1240, New Orleans, Louisiana. Association for Computational Linguistics.
  * Chen and Cardie (2018b) Xilun Chen and Claire Cardie. 2018b.  Multinomial adversarial networks for multi-domain text classification.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 1226–1240, New Orleans, Louisiana. Association for Computational Linguistics.
  * Chen et al. (2018) Xilun Chen, Yu Sun, Ben Athiwaratkun, Claire Cardie, and Kilian Weinberger. 2018\.  Adversarial deep averaging networks for cross-lingual sentiment classification.  volume 6, pages 557–570.
  * Corney et al. (2016) D. Corney, M. Albakour, Miguel Martinez-Alvarez, and Samir Moussa. 2016.  What do a million news articles look like?  In _NewsIR@ECIR_.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Ganin and Lempitsky (2015) Yaroslav Ganin and Victor Lempitsky. 2015.  Unsupervised domain adaptation by backpropagation.  In _Proceedings of the 32nd International Conference on International Conference on Machine Learning - Volume 37_ , ICML’15, page 1180–1189. JMLR.org.
  * Go et al. (2009) Alec Go, Richa Bhayani, and Lei Huang. 2009.  Twitter sentiment classification using distant supervision.  _CS224N project report, Stanford_ , 1(12):2009.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014.  In _NIPS_.
  * Hu and Liu (2004) Minqing Hu and Bing Liu. 2004.  Mining and summarizing customer reviews.  In _Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’04, page 168–177, New York, NY, USA. Association for Computing Machinery.
  * Kingma and Ba (2014) Diederik Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _International Conference on Learning Representations_.
  * Lin et al. (2016) Yankai Lin, Shiqi Shen, Zhiyuan Liu, Huanbo Luan, and Maosong Sun. 2016.  Neural relation extraction with selective attention over instances.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 2124–2133.
  * Loshchilov and Hutter (2018) Ilya Loshchilov and Frank Hutter. 2018.  Fixing weight decay regularization in adam.
  * Manning et al. (2014) Christopher D Manning, Mihai Surdeanu, John Bauer, Jenny Rose Finkel, Steven Bethard, and David McClosky. 2014.  The stanford corenlp natural language processing toolkit.  In _Proceedings of 52nd annual meeting of the association for computational linguistics: system demonstrations_ , pages 55–60.
  * Mintz et al. (2009) Mike Mintz, Steven Bills, Rion Snow, and Dan Jurafsky. 2009.  Distant supervision for relation extraction without labeled data.  In _Proceedings of the Joint Conference of the 47th Annual Meeting of the ACL and the 4th International Joint Conference on Natural Language Processing of the AFNLP_ , pages 1003–1011.
  * Miyato et al. (2017) Takeru Miyato, Andrew M. Dai, and Ian Goodfellow. 2017.  Adversarial training methods for semi-supervised text classification.
  * Qin et al. (2018) Pengda Qin, Weiran Xu, and William Yang Wang. 2018.  DSGAN: Generative adversarial training for distant supervision relation extraction.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 496–505, Melbourne, Australia. Association for Computational Linguistics.
  * Ratner et al. (2020) Alexander Ratner, Stephen H. Bach, Henry R. Ehrenberg, Jason A. Fries, Sen Wu, and Christopher Ré. 2020.  Snorkel: rapid training data creation with weak supervision.  _VLDB J._ , 29(2-3):709–730.
  * Sanh et al. (2019) Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf. 2019.  Distilbert, a distilled version of BERT: smaller, faster, cheaper and lighter.  _CoRR_ , abs/1910.01108.
  * Sedova et al. (2021) Anastasiya Sedova, Andreas Stephan, Marina Speranskaya, and Benjamin Roth. 2021\.  Knodle: Modular weakly supervised learning with pytorch.  _CoRR_ , abs/2104.11557.
  * Snoek et al. (2012) Jasper Snoek, Hugo Larochelle, and Ryan Adams. 2012.  Practical bayesian optimization of machine learning algorithms.  In _Proc. NIPS_.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. 2014.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958.
  * Takamatsu et al. (2012) Shingo Takamatsu, Issei Sato, and Hiroshi Nakagawa. 2012.  Reducing wrong labels in distant supervision for relation extraction.  In _Proceedings of the 50th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 721–729.
  * Wu et al. (2017) Yi Wu, David Bamman, and Stuart Russell. 2017.  Adversarial training for relation extraction.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , pages 1778–1783.
  * Yeh (2000) Alexander Yeh. 2000.  More accurate tests for the statistical significance of result differences.  In _COLING 2000 Volume 2: The 18th International Conference on Computational Linguistics_.
  * Zeng et al. (2018) Daojian Zeng, Yuan Dai, Feng Li, R Simon Sherratt, and Jin Wang. 2018.  Adversarial learning for distant supervised relation extraction.  _Computers, Materials & Continua_, 55(1):121–136.

##  Appendix A Appendix

###  A.1 Dataset statistics

The datasets used for the KnowMAN experiments have different properties.
Especially the numer of labeling functions and the dataset sizes varies a lot.

dataset | classes | train/test samples | lfs
---|---|---|---
Spam | 2 | 1586/250 | 10
Spouse | 2 | 22254/2701 | 9
IMDb | 2 | 40000/5000 | 6786
Table 2: Dataset statistics for KnowMAN experiments. Lfs are labeling
functions.

###  A.2 Hyperparameter optimization

We perform hyperparameter tuning using Bayesian optimization (Snoek et al.,
2012). Bayesian Optimization is an approach that uses the Bayes Theorem to
direct the search in order to find the minimum or maximum of a black-box
objective function. In comparison with random search and grid search, it tends
to obtain better hyperparameters in fewer steps by making a proper balance
between exploration and exploitation steps. Our hyperparameter space includes
batch size, dropout, number of iterations over $\mathcal{D}$, the shared
hidden size of the models, learning rate for $\mathcal{D}$ and
$\mathcal{F}_{s},\mathcal{C}$ and the number of layers of
$\mathcal{C},\mathcal{D}$ and $\mathcal{F}_{s}$. We implemented two ways of
evaluation and best model saving during training: i) evaluate after each batch
and save the best model, ii) evaluate after a certain number of steps in
between the batches and save the best model. We also optimized the number of
steps if logging in between a batch.

We evaluated the models for IMDb and Spouse on the respective validation set.
For the Spam dataset, there is no development set available and we used the
following hyperparameters for KnowMAN TF-IDF following the parameters used in
Chen and Cardie (2018b): Batch size: 32, dropout: 0.4, n critic: 5, lambda:
2.0, shared hidden size: 700, learning rate C & F: 0.0001, learning rate D:
0.0001 , number of F layers: 1, number of C layers: 1, number of D layers: 1.

Figure 2: Sampled hyperparameters for KnowMAN TF-IDF on IMDb. Optimal
hyperparameters are indicated in red.
Batch size: 895, dropout: 0.275, n critic: 50, lambda: 4.9, shared hidden
size: 585, learning rate C & F: 0.0001, learning rate D: 0.0001, number of F
layers: 1 , number of C layers: 1, number of D layers: 10\.
Histograms on the diagonal show how, for each hyperparameter, how many samples
have been drawn during optimization. Figure 3: Sampled hyperparameters for
KnowMAN DistilBERT on Spouse. Optimal hyperparameters are indicated in red.
Batch size: 16, dropout: 0.379, n critic: 1, lambda: 5.0, shared hidden size:
988, learning rate C & F: 0.0005, learning rate D: 0.001 , number of F layers:
5, number of C layers: 10, number of D layers: 1\.
Histograms on the diagonal show how, for each hyperparameter, how many samples
have been drawn during optimization.

###  A.3 Experimental details

We ran our experiments on a DGX-1 server with one V100 GPU per experiment. The
runtime of one model depends on the dataset: 0.25 hours for the Spam dataset,
0.25 hours for the Spouse dataset, and 8 hours for the IMDb dataset.

Please find our implementation at https://github.com/LuisaMaerz/KnowMAN.