A global synchronization theorem for oscillators on a random graph

Networks of coupled oscillators have long been studied in biology, physics, engineering, and nonlinear dynamics. winfree1967biological ; winfree1980geometry ; kuramoto1984chemical ; mirollo1990synchronization ; watanabe1994constants ; pikovsky2003synchronization ; strogatz2003sync ; dorfler2014synchronization ; pikovsky2015dynamics ; rodrigues2016kuramoto ; abrams2016introduction ; boccaletti2018synchronization ; matheny2019exotic Recently they have begun to attract the attention of other communities as well. For example, oscillator networks have been recognized as having the potential to yield “beyond-Moore’s law” computational devices CsabaPorod for graph coloring, Wu image segmentation, Novikov and approximate maximum graph cuts. Steinerberger They have also become a model problem for understanding the global convergence of gradient descent in nonlinear optimization. ling2018landscape

In such settings, global issues come to the fore. When performing gradient descent, for instance, one typically wants to avoid getting stuck in local minima. Conditions to enforce convergence to the global minimum then become desirable. Likewise, if an oscillator network has only a single, globally attracting state, we know exactly how the system will behave in the long run.

We say that a network of oscillators globally synchronizes if it converges to a state for which all the oscillators are in phase, starting from all initial conditions except a set of measure zero. Until recently, only a few global synchronization results were known for networks of oscillators. mirollo1990synchronization ; watanabe1994constants These results were restricted to complete graphs, in which each oscillator is coupled to all the others. In the past decade, however, several advances have been made for a wider class of network structures, starting with work by Taylor, taylor2012there who proved that if each oscillator in a network of identical Kuramoto oscillators is connected to at least 94 percent of the others, the network will fall into perfect synchrony for almost all initial conditions, no matter what the topology of the network is like in other respects. Taylor’s result was strengthened by Ling, Xu, and Bandeira ling2018landscape in 2018 and further progress has been made since then. lu2020synchronization ; kassabovtownsendstrogatz

Ling et al. ling2018landscape also made a seminal advance in the study of random networks. They considered identical Kuramoto oscillators on an Erdős–Rényi random graph, ErdosRenyi in which any two oscillators are coupled with unit strength, independently and at random, with probability $0\leq p\leq 1$; otherwise they are uncoupled. They showed that if $p\gg\log(n)/n^{1/3}$, then with high probability the network is globally synchronizing as $n\rightarrow\infty$.

The open question is to find and prove the sharpest result along these lines. Intuitively, as one increases the value of $p$ from $0$ to $1$, one expects to find a critical threshold above which global synchrony is extremely likely and below which it is extremely rare. At the very least, for global synchrony to be ubiquitous, $p$ must be large enough to ensure that the random network is connected, and from the theory of random graphs ErdosRenyi we know that connectedness occurs with high probability once $p>(1+\varepsilon)\log(n)/n$ for any $\varepsilon>0$. So the critical threshold for global synchronization cannot be any smaller than this connectivity threshold, and is apt to be a little above. On the basis of numerical evidence, Ling et al. ling2018landscape conjectured, and we agree with them, that if $p\gg\log(n)/n$ then Erdős–Rényi graphs are globally synchronizing as $n\rightarrow\infty$, but nobody has proven that yet. The challenge is to see how close one can get. In this paper, we come within a factor of $\log(n)$ and prove that $p\gg\log^{2}(n)/n$ is sufficient to give global synchrony with high probability. With the aid of a computer, we have convincing evidence that $p\gg\log(n)/n$ is sufficient as $n\rightarrow\infty$ (see Section VI.1).

We study the homogeneous Kuramoto model jadbabaie2004stability ; wiley2006size ; taylor2012there in which each oscillator has the same frequency $\omega$. By going into a rotating frame at this frequency, we can set $\omega=0$ without loss of generality. Then phase-locked states in the original frame correspond to equilibrium states in the rotating frame. So, to explore the question that concerns us, it suffices to study the following simplified system of identical Kuramoto oscillators:

where $\theta_{j}(t)$ is the phase of oscillator $j$ (in the rotating frame) and the adjacency matrix $A$ is randomly generated. In particular, $A_{jk}=A_{kj}=1$ with probability $p$, with $A_{jk}=A_{kj}=0$ otherwise, independently for $1\leq j,k\leq n$. Thus, all interactions are assumed to be symmetric, equally attractive, and of unit strength. We take the unusual convention that the network can have self-loops so that oscillator $i$ is connected to itself with probability $p$, i.e., $\mathbb{P}\left[A_{jj}=1\right]=p$. Since $\sin(0)=0$, this convention does not alter the dynamics of the oscillators but it does make proof details easier to write down.

Finally, since the adjacency matrix $A$ is symmetric, we know that (1) is a gradient system. wiley2006size ; jadbabaie2004stability Thus, all the attractors of (1) are equilibrium points, which means we do not need to concern ourselves with the possibility of more complicated attracting invariant sets like limit cycles, tori, or strange attractors.

We find it helpful to visualize the oscillators on the unit circle, where oscillator $j$ is positioned at the coordinate $(\cos(\theta_{j}),\sin(\theta_{j}))$. From this perspective, a network is globally synchronizing if starting from any initial positions (except a set of measure zero), all the oscillators eventually end up at the same point on the circle. Due to periodicity, we assume that the phases, $\theta_{j}$, take values in the interval $[-\pi,\pi)$.

One cannot determine whether a network is globally synchronizing by numerical simulation of (1), as it is impossible to try all initial conditions. Of course, one can try millions of random initial conditions of the oscillators’ phases and then watch the dynamics of (1). But even if all observed initial states eventually fall into the all-in-phase state, one cannot conclude the system is globally synchronizing because other stable equilibria could still exist; their basins might be minuscule but could nevertheless have positive measure.

With that caveat in mind, we note that such numerical experiments have been conducted, and they suggest that $p\gg\log(n)/n$ is sufficient for global synchronization. ling2018landscape In this paper, we investigate global synchrony via theoretical study. We show that $p\gg\log^{2}(n)/n$ is good enough to ensure global synchronization with high probability, improving on $p\gg\log(n)/n^{1/3}$ proved by Ling et al., ling2018landscape and bringing us closer to the connectivity threshold of Erdős–Rényi graphs.

Although we focus on Erdős–Rényi graphs, many of the inequalities we derive hold for any random or deterministic network. To highlight this, we state many of our findings for a general graph $G$ and a general parameter $p\in\mathbb{R}$. In the end, we restrict ourselves back to Erdős–Rényi graphs and take $p$ to be the probability of a connection between any two oscillators.

Our results depend on both the adjacency matrix $A$ and the graph Laplacian matrix $L=D-A$, where $D$ is a diagonal matrix and $D_{ii}$ is the degree of vertex $i$ (counting self-loops). For any $p\in\mathbb{R}$, denote the shifted adjacency and graph Laplacian matrix by

where $J_{n}$ is the $n\times n$ matrix of all ones and $I_{n}$ is the $n\times n$ identity matrix. It is worth noting that for Erdős–Rényi graphs, the shifts are precisely the expectation of the matrices as $\mathbb{E}\left[A\right]=pJ_{n}$ and $\mathbb{E}\left[L\right]=pJ_{n}-npI_{n}$. Remarkably, we show that the global synchrony of a network can be guaranteed by ensuring that the spectral norms $\|\Delta_{A}\|$ and $\|\Delta_{L}\|$ satisfy particular inequalities. The spectral norm of a symmetric matrix is the maximum eigenvalue in absolute magnitude, and $\|\Delta_{A}\|$ and $\|\Delta_{L}\|$ are extensively studied in the random matrix literature. FK ; Vu What is remarkable is that no other information about the network’s structure is needed; the norms of these two matrices alone encapsulate the structural aspects that matter for global synchronization. We also find it appealing that the spectral norm of the graph Laplacian matrix appears naturally in our analysis, as it has been used previously to study the dynamics of networks of oscillators Pecora as well as diffusion on graphs.

While we focus in this paper on global synchronization of identical oscillators, Medvedev and his collaborators have considered other aspects of synchronization for non-identical oscillators on Erdős–Rényi graphs, as well as on other graphs such as Cayley and Watts–Strogatz graphs, often in the continuum limit. Medvedev1 ; Medvedev2 ; Medvedev3 ; ChibaMedvedevMizuhara Their findings have a similar overarching message that synchronization occurs spontaneously above a critical threshold.

To prove that a random graph is globally synchronizing with high probability, we bridge the gap between spectral graph theory and coupled oscillator theory. The literature contains many good probabilistic estimates for the spectral norm of a random graph’s adjacency and graph Laplacian matrices, which we use to control the long-time dynamics of (1). The key to our proof is to establish conditions on these two spectral norms that force any stable equilibrium to have phases that lie within a half-circle. Confining the phases in this way then guarantees global synchronization, because it is known that the only stable equilibrium of (1) with phases confined to a half-circle is the all-in-phase synchronized state. jadbabaie2004stability ; dorfler2014synchronization

Our first attempt at controlling the phases is the inequality we state below as (9), which can be used to guarantee that very dense Erdős–Rényi graphs are globally synchronizing with high probability; as an aside, when $p=1$ this inequality provides a new proof that a complete graph of identical Kuramoto oscillators is globally synchronizing. watanabe1994constants A similar inequality to (9) was derived by Ling, Xu, and Bandeira ling2018landscape to show that $p\gg\log(n)/n^{1/3}$ is sufficient for global synchrony with high probability.

To improve on (9), our argument becomes more intricate. We carefully examine the possible distribution of edges between oscillators whose phases lie on different arcs of the circle, and show that any equilibrium is destabilized if there are too many edges between oscillators that have disparate phases (see Lemma 8 and Figure 2).

An important quantity in the study of Kuramoto oscillators is the so-called complex order parameter, $\rho_{1}$. The magnitude of $\rho_{1}$ is between $0$ and $1$ and measures the synchrony of the oscillators in the network. We find it useful to also look at second-order moments of the oscillator distribution for analyzing the synchrony of random networks. Higher-order moments are also called Daido order parameters and can be used to analyze oscillators with all-to-all coupling, corresponding to a complete graph. daido1992order ; Perez_97 ; pikovsky2015dynamics ; terada2020linear

For an equilibrium $\theta=(\theta_{1},\ldots,\theta_{n})$, we define the first- and second-order moments as

(For convenience, we use the notation $\sum_{j}$ to mean $\sum_{j=1}^{n}$.) Without loss of generality, we may assume that the complex order parameter $\rho_{1}$ is real-valued and non-negative. To see this, write $\rho_{1}=|\rho_{1}|e^{i\psi}$ for some $\psi$. Then, $\hat{\theta}=(\theta_{1}-\psi,\ldots,\theta_{n}-\psi)$ is also an equilibrium of (1) with the same stability properties as $\theta$ since (1) is invariant under a global shift of all phases by $\psi$. Therefore, for the rest of this paper, we assume that $\psi=0$ for any equilibrium of interest with $0\leq\rho_{1}\leq 1$. When $\rho_{1}=1$, the oscillators are in pure synchrony and when $\rho_{1}\approx 0$, the phases are scattered around the unit circle without a dominant phase.

To avoid working with complex numbers, it is convenient to consider the quantity $|\rho_{2}|^{2}$. For $m=1,2$, we have

By analyzing $\rho_{1}^{2}$ and $|\rho_{2}|^{2}$, one hopes to witness the rough statistics of an equilibrium to understand its potential for synchrony without concern for its precise pattern of phases.

Let $q_{\theta}=(e^{i\theta_{1}},\dots,e^{i\theta_{n}})^{\top}$ and note that

where $A$ is the adjacency matrix, $\overline{q_{\theta}}$ is the complex conjugate of $q_{\theta}$, and ^⊤ denotes the transpose of a vector. Since $\cos^{2}(\theta_{k}-\theta_{j})=\frac{1}{2}(\cos(2(\theta_{k}-\theta_{j}))+1)$, we have

where $\mathbf{1}$ is the vector of all ones and $q_{2\theta}=(e^{i2\theta_{1}},\dots,e^{i2\theta_{n}})^{\top}$.

We would like to know when a stable equilibrium is close to the all-in-phase state, and we know this when $\rho_{1}$ is close to $1$. Therefore, we begin by deriving a lower bound on $\rho_{1}$ for any stable equilibrium of (1). Similar inequalities involving $\rho_{1}^{2}$ and $|\rho_{2}|^{2}$ have been used to demonstrate that sufficiently dense Kuramoto networks are globally synchronizing. ling2018landscape ; kassabovtownsendstrogatz

To maximize the lower bound on $\rho_{1}^{2}$ from Lemma 1, one can optimize over $p$. For a random graph where each edge has a fixed probability of being present, independently of the other edges, we usually just select $p$ to be that probability. Regardless, to make the lower bound on $\rho_{1}^{2}$ in Lemma 1 useful, we need to find a nontrivial lower bound on $|\rho_{2}|^{2}$, since this quantity appears in the right hand side of (4). To obtain such a bound we use the following technical lemma.

To see how Lemma 2 can be used to derive a lower bound on $|\rho_{2}|^{2}$ for any stable equilibrium state, we start by dropping the second sum in Lemma 2. By using the upper bound $\|\Delta_{A}q_{\theta}\|^{2}\leq n\|\Delta_{A}\|^{2}$, we find that $\sum_{j}\sin^{2}(\theta_{j})\leq\|\Delta_{A}\|^{2}/(np^{2}\rho_{1}^{2})$. Since $n{\rm Re}\!\left(\rho_{2}\right)=\sum_{j}\cos(2\theta_{j})=\sum_{j}\left(1-2\sin^{2}(\theta_{j})\right)$, we have the following lower bound on $|\rho_{2}|$:

From (4) and (7), we observe that when $\|\Delta_{A}\|\ll np$, then $\rho_{1}$ and $|\rho_{2}|$ must both be close to $1$. Intuitively, this should mean that the corresponding stable equilibrium, $\theta$, must be close to the all-in-phase state with the possible exception of a small number of stray oscillators. However, our goal is to prove global synchrony, which is a more stringent condition, and we must completely rule out the existence of stray oscillators.

To precisely control the number of stray oscillators, we define a set of indices for oscillators whose phases lie outside of a sector of half-angle $\phi$ (centered about the all-in-phase state):

for any angle $0\leq\phi\leq\pi$ (see Figure 1). If we can prove that $C_{\pi/2}$ is empty, then we know that all the phases in the equilibrium state lie strictly inside a half-circle. That would give us what we want, because by a basic theorem for the homogeneous Kuramoto model, the only stable equilibrium $\theta$ with all its phases confined to a half-circle is the all-in-phase state. jadbabaie2004stability ; dorfler2014synchronization In terms of bounds, since $|C_{\pi/2}|$ is integer-valued, if we can show that $|C_{\pi/2}|<1$ then we know that $C_{\pi/2}$ is the empty set.

From Lemma 2 and $\|\Delta_{A}q_{\theta}\|^{2}\leq n\|\Delta_{A}\|^{2}$, we find that

Therefore, by plugging in $\phi=\pi/2$, we see that $|C_{\pi/2}|\leq\|\Delta_{A}\|^{2}/(np^{2}\rho_{1}^{2})$. Thus, if $\|\Delta_{A}\|^{2}<np^{2}\rho_{1}^{2}$ then $C_{\pi/2}$ must be the empty set and the network is globally synchronizing.

Unfortunately, this kind of reasoning is not sufficient to prove global synchrony for graphs of interest to us here. For example, for an Erdős–Rényi random graph, we know that $\|\Delta_{A}\|^{2}\approx 4p(1-p)n$ with high probability for large $n$ (see Section VI). So the upper bound on $|C_{\pi/2}|$ is approximately $4(1-p)/(p\rho_{1}^{2})$, which for $p<1/2$ is certainly not good enough to conclude that $C_{\pi/2}$ is empty. Instead, we must further improve our bounds on $|C_{\pi/2}|$ by using a recursive refinement strategy that we refer to as an “amplification" argument (see Section V).

The precise amplification argument that we use requires bounds on the number of edges and sizes of vertex sets of a graph expressed in terms of the spectral norms of $\Delta_{A}$ and $\Delta_{L}$. It is worth noting that these bounds hold for both deterministic and random graphs. For a vertex set $C$ of a graph $G$, we denote the characteristic vector by $v_{C}$, i.e., $(v_{C})_{j}=1$ if $j\in C$ and $(v_{C})_{j}=0$ if $j\not\in C$. We use $v_{C}^{\top}$ to denote the vector transpose of $v_{C}$. We write $|C|$ to be the number of vertices in $C$ and denote the number of edges in $G$ between two vertex sets $C$ and $C^{\prime}$ as $E_{C,C^{\prime}}$. For Erdős–Rényi graphs, one expects to have $E_{C,C^{\prime}}\approx p|C||C^{\prime}|$. However, for our argument to work, expectations are not adequate; instead we need to have bounds on the difference between $E_{C,C^{\prime}}$ and $p|C||C^{\prime}|$. The results in this section are proved by classical techniques and closely related bounds are well known. We give the proofs of the precise inequalities that we need so that the paper is self-contained.

We first show that for any vertex set $C$, the number of edges connecting a vertex in $C$ to another vertex in $C$ deviates from $p|C|^{2}$ by at most $\|\Delta_{A}\||C|$, for any $0\leq p\leq 1$.

Lemma 3 controls the number of edges connecting a set of vertices. The following result controls the number of edges leaving a vertex set. We denote the vertices of $G$ that are not in $C_{1}$ by $V_{G}\setminus C_{1}$.

When a bound on $\|\Delta_{L}\|$ is available, Lemma 4 can be used to ensure that $G$ is connected. In particular, Lemma 4 tells us that a graph of size $n$ is connected if $\|\Delta_{L}\|<np$ for some $0\leq p\leq 1$.

Since $0\leq E_{C_{1},C_{2}}\leq|C_{1}||C_{2}|$, we know that $\left|E_{C_{1},C_{2}}-p|C_{1}||C_{2}|\right|\leq\max\{p,1-p\}|C_{1}||C_{2}|$. Therefore, Lemma 4 is only a useful bound when $\|\Delta_{L}\|\ll np$. We can take Lemma 4 a step further and bound the number of edges between any two sets of vertices of $G$ using $\|\Delta_{L}\|$. We denote the union of the vertex sets $C_{1}$ and $C_{2}$ by $C_{1}\sqcup C_{2}$.

Lemma 3 and Lemma 5 can be combined to obtain the following result.

We are finally ready for our amplification argument, which is a way to improve the bounds on $|C_{\pi/2}|$ from (9). We first write down a new inequality that holds for any stable equilibrium. We write it down using a kernel function $K$ that later allows us to improve our argument with the aid of a computer (see Section VI.1).