Directed acyclic decomposition of Kuramoto equations

TL;DR

This paper introduces a framework for decomposing Kuramoto oscillator networks into directed acyclic subnetworks, enabling a divide-and-conquer approach to analyze complex synchronization behaviors in large, heterogeneous networks.

Contribution

It presents a novel method for decomposing Kuramoto networks into directed acyclic subnetworks, facilitating analysis of synchronization configurations.

Findings

Decomposition method for Kuramoto networks into acyclic subnetworks

Framework enables scalable analysis of large networks

Supports divide-and-conquer approach for synchronization study

Abstract

The Kuramoto model is one of the most widely studied model for describing synchronization behaviors in a network of coupled oscillators, and it has found a wide range of applications. Finding all possible frequency synchronization configurations in a general non-uniform heterogeneous sparse network is an important yet difficult problem due to the complicated nonlinear interactions. In this paper, we develop a general framework for decomposing a Kuramoto network into smaller directed acyclic subnetworks that will form the foundation of a divide-and-conquer approach for studying the frequency synchronization configurations of large Kuramoto networks.

Tables

Table 1. Table 1: Number of facet subnetworks and the root counts contributed by each classified by the size (n.o. edges) of subnetworks

Facet subnetwork size	9	10	11	12	13
N.o. facet subnetworks	116	414	540	384	144
Root count per subnetwork	1	2	4	8	16

Table 2. Table 2: Time consumed by a preliminary implementation based on a modified version of Hom4PS-3 for solving the synchronization system derived from the Kuramoto network displayed in Figure   8 using different numbers of processor cores.

Threads/cores used	1	2	3	4	6	8
Wall clock time	9.1s	4.6s	3.3s	2.6s	1.8s	1.4s

Equations

$${\omega }_i-\sum _{j\in {\mathcal{N}}_G(i)}{k}_{ij}\sin ({\theta }_i-{\theta }_j)=0\text{for }i=1,\dots ,n.$$

$${\dot{\theta }}_i(t)={\omega }_i-\sum _{j\in {\mathcal{N}}_G(i)}{k}_{ij}\sin ({\theta }_i(t)-{\theta }_j(t)),\text{for }i=0,\dots ,n,$$

$${\omega }_i-\sum _{j\in {\mathcal{N}}_G(i)}{k}_{ij}\sin ({\theta }_i(t)-{\theta }_j(t))=0\text{for }i=0,\dots ,n.$$

$$x_i=e^{\mathbf{i}{z}_i}=e^{r_i+\mathbf{i}{\theta }_i}\text{for }i=1,\dots ,n\text{ and }{x}_0=e^{0+0\mathbf{i}}=1.$$

$$\sin (z_i-z_j)=\frac{e^{z_i\mathbf{i}}{e}^{-z_j\mathbf{i}}-e^{z_j\mathbf{i}}{e}^{-z_i\mathbf{i}}}{2\mathbf{i}}=\frac{1}{2\mathbf{i}}\left(\frac{x_i}{x_j}-\frac{x_j}{x_i}\right),$$

$${\omega }_i-\sum _{j\in {\mathcal{N}}_G(i)}{a}_{ij}^{\prime }\left(\frac{x_i}{x_j}-\frac{x_j}{x_i}\right)=0\text{for }i=1,\dots ,n$$

$$c_k-\sum _{(i,j)\in \mathcal{E}(G)}{a}_{ijk}\left(\frac{x_i}{x_j}-\frac{x_j}{x_i}\right)=0\text{for }k=1,\dots ,n$$

$$c_1$$

$$c_2$$

$$c_1$$

$$c_2$$

$$\{\begin{array}{r}\frac{c_1}{t}-\left[a_{101}\left(\frac{x_1}{x_0}-\frac{x_0}{x_1}\right)+a_{121}\left(\frac{x_1}{x_2}-\frac{x_2}{x_1}\right)+a_{201}\left(\frac{x_2}{x_0}-\frac{x_0}{x_2}\right)\right]\\ \frac{c_2}{t}-\left[a_{102}\left(\frac{x_1}{x_0}-\frac{x_0}{x_1}\right)+a_{122}\left(\frac{x_1}{x_2}-\frac{x_2}{x_1}\right)+a_{202}\left(\frac{x_2}{x_0}-\frac{x_0}{x_2}\right)\right]\end{array}.$$

$$x_1$$

$$x_{2}t$$

$${\nabla }_G=\mathrm{conv}\left\{{\mathbf{e}}_i-{\mathbf{e}}_j\mid (i,j)\in \mathcal{E}(G)\right\},$$

$${\mathcal{V}}_F$$

$${\mathcal{E}}_F$$

$$c_k-\sum _{(i,j)\in {\mathcal{E}}_F}{a}_{kij}\left(\frac{x_i}{x_j}\right)=0\text{for each }k\in {\mathcal{V}}_F.$$

$$a_{ij}^{\prime }\frac{x_i}{x_j}=\frac{k_{ij}}{2\mathbf{i}}{e}^{\mathbf{i}({\theta }_i-{\theta }_j)}=\frac{k_{ij}}{2}[\sin ({\theta }_i-{\theta }_j)-\mathbf{i}\cos ({\theta }_i-{\theta }_j)]$$

$$⟨\mathbf{p},\alpha ⟩$$

$$⟨{\mathbf{e}}_{i_k}-{\mathbf{e}}_{i_{k+1}},\mathbf{\alpha }⟩={\alpha }_{i_k}-{\alpha }_{i_{k+1}}=h<0\text{for each }k=1,\dots ,m,$$

$${\alpha }_{i_1}<{\alpha }_{i_2}<\cdots <{\alpha }_{i_m}<{\alpha }_{i_1},$$

$$V_1$$

$$V_2$$

$$⟨{\mathbf{e}}_{i_k}-{\mathbf{e}}_{i_{k+1}},\mathbf{\alpha }⟩=h\text{for each }k=0,\dots ,m-1$$

$$⟨{\mathbf{e}}_{i_0}-{\mathbf{e}}_{i_m},\mathbf{\alpha }⟩=mh.$$

$$⟨{\mathbf{e}}_{j_0}-{\mathbf{e}}_{j_{\ell }},\mathbf{\alpha }⟩=\ell h.$$

$$⟨{\mathbf{e}}_{i_k}-{\mathbf{e}}_{i_{k+1}},\mathbf{\alpha }⟩=h\text{for each }k=0,\dots ,m-1.$$

$$⟨{\mathbf{e}}_{i_1}-{\mathbf{e}}_{i_m},\mathbf{\alpha }⟩=mh.$$

$$⟨{\mathbf{e}}_{i_{k+1}}-{\mathbf{e}}_{i_k},\mathbf{\alpha }⟩\ge h\text{for each }k=m,\dots ,\ell -1.$$

$$⟨{\mathbf{e}}_{i_1}-{\mathbf{e}}_{i_m},\mathbf{\alpha }⟩\ge (\ell -m)h$$

$$⟨{\mathbf{e}}_j-{\mathbf{e}}_i,-\mathbf{\alpha }⟩=⟨{\mathbf{e}}_i-{\mathbf{e}}_j,\mathbf{\alpha }⟩=h$$

Full text

The spontaneous synchronization of a network of oscillators is an emergent phenomenon that naturally appears in many seemingly independent complex systems including mechanical, chemical, biological, and even social systems. The Kuramoto model is one of the most widely studied and successful mathematical models for analyzing synchronization behaviors. While much is known about the macro-scale question of whether or not a Kuramoto network can be synchronized, detailed analysis of the possible configurations of the oscillator once it has reached synchronization remains difficult for large networks partly due to the nonlinear interactions involved. In this work, we demonstrate that by dividing the link between two oscillators into two one-way interactions, complex networks can indeed be decomposed into much simpler subnetwork. This is a crucial step toward fully understanding synchronization configurations in large networks.

I Introduction

Mathematical modeling and analysis of spontaneous synchronization have found many important applications in physics, chemistry, engineering, biology, and medical science. Dörfler and Bullo (2014) Originally introduced to describe chemical oscillators, the Kuramoto model Kuramoto (1975, 2012) has become one of the most widely studied models for describing synchronization behaviors in a wide range of situations. This paper focuses on the study of frequency synchronization, which describes a particular type of synchronization behavior where oscillators are tuned into the same frequency. The central objective is to understand the set of all such configurations on potentially non-uniform and non-homogeneous Kuramoto network. With the appropriate frame of reference, this problem is equivalent to the study of the full set of solutions of the system of nonlinear equations

[TABLE]

The main contribution of this paper is the development of a general framework for decomposing a Kuramoto network into subnetworks that can be studied more easily. This decomposition is enabled by the key insights gained through an abstract polytope that encodes the network topology. The way in which a network is broken down into subnetworks mirrors the process through which how the boundary of that polytope is broken down into facets. From the view point of dynamics, this decomposition comes from the limit behavior of Kuramoto networks as the differences of the oscillators’ natural frequencies are amplified to infinity. The framework presented in this paper marks a crucial step toward a divide-and-conquer approach to studying large Kuramoto networks.

Using a complex algebraic formulation, we deform the synchronization equations under which a Kuramoto network degenerates into a union of those of simpler subnetworks supported by directed acyclic subgraphs. By extending our domain to complex phase angles, we can form a new rational system that also captures the synchronization configurations as solutions. This procedure, detailed in section III, allows for the introduction of powerful tools stemming from complex algebraic geometry. The main ideas are briefly illustrated through a simple example in section IV. Then, in section V, through a construction known as the adjacency polytope, we decompose the Kuramoto network into simpler subnetworks induced by facets of this polytope. Each “facet subnetwork” corresponds to a directed acyclic subnetwork of the original network. We also explore the topological properties of these subnetworks. These “facet subnetworks” preserves and reveals many important properties of the original network. We demonstrate in theorem 2 that the number of synchronization configurations that the original network has is bounded by the total root count of the facet subnetworks. In theorem 3, we demonstrate that the decomposition of a network into facet subnetworks can be understood as a a smooth deformation of the synchronization equations that can degenerate into facet subsystems. Among the facet subnetworks, simplest type is known as a “primitive subnetwork” and is analyzed in section VI. Each primitive subnetwork has a unique synchronization configuration, and hence form the basic building block of the decomposition. Finally, we illustrate the decomposition scheme via concrete examples in section VII and conclude with remarks on future directions in section VIII.

II Kuramoto model

The Kuramoto model Kuramoto (1975) is a mathematical model used to study behavior of a network of coupled oscillators. An oscillator is simply an object that can continuously vary between two states. In isolation, each oscillator has its own natural frequency. When we consider networks of coupled oscillators, however, rich and complex dynamic behaviors emerge. The oscillators are coupled with one another by idealized springs with stiffness characterized by their coupling strength. For a pair of oscillators $(i,j)$, the real number $k_{ij}$ quantifies the strength of coupling between them. The topology of a network of $N=n+1$ oscillators is modeled by a graph $G=(V,E)$, in which nodes $V=\{0,\dots,n\}$ represent the oscillators and edges $E$ represent their connections. The coupling strengths $K=\{k_{ij}\}$ and the natural frequencies $\omega=(\omega_{0},\dots,\omega_{n})$ capture the quantitative information of the network. In its simplest form, the Kuramoto model is a differential equation that describes the nonlinear interactions in such a network given by

[TABLE]

where $\theta_{i}(t)$ is the phase angle of the oscillator $i$ as a function of time $t$, and $\mathcal{N}_{G}(i)$ denotes the set of its neighboring nodes.

The network described by $(G,K,\omega)$ together with the above differential equation is referred to as a Kuramoto network. Central to this paper are the special configurations in which the angular velocity of all oscillators become perfectly aligned, known as frequency synchronization configurations. That is, $\frac{d\theta_{i}}{dt}=c$ for $i=0,\dots,n$ and a constant $c$. By adopting a proper frame of reference, we can assume $c=0$, and the (frequency) synchronization configurations are simply equilibria of the Kuramoto model (1) given by

[TABLE]

For the removal of the inherent degree of freedom induced by uniform translation, the standard practice is to choose node 0 to be the reference node and set $\theta_{0}=0$. Assuming the couplings are symmetric, i.e., $k_{ij}=k_{ji}$, the $n+1$ equations are then linearly dependent. This allows for the elimination of one equation and produces a system of $n$ equations in $n$ unknowns $\theta_{1},\dots,\theta_{n}$. This synchronization system is the main focus of this paper. The structure of the solution set to this system has been the subject of intense research since the 1970s. Baillieul and Byrnes (1982); Kuramoto (1975) While earlier studies have focused on macro-scale and statistical analyses of large (often infinite) networks,Kuramoto (1984) recent research has gradually shifted toward precise analysis, of small finite networks.Strogatz (2000); Aeyels and Rogge (2004); Mirollo and Strogatz (2005); Mehta et al. (2015); Chen et al. (2016); Rogge and Aeyels (2004); Ochab and Góra (2010); Delabays, Coletta, and Jacquod (2016, 2017); Manik, Timme, and Witthaut (2017) The complex algebraic approachBaillieul and Byrnes (1982) started by Baillieul and Byrnes has been particularly successful in this regard and forms the foundation of this work.

The key question we set out to answer is whether or not the full set of solutions to this system can be understood through a study of the simpler subnetworks of the original Kuramoto network.

III Complex algebraic formulation

Even though the original formulation of the synchronization system (2) considers only real phase angles, it is useful to expand the domain to more general complex phase angles, as this would allow us to apply powerful tools from complex algebraic geometry. Consider complex phase angle $z_{i}=\theta_{i}-r_{i}\mathbf{i}$ for $i=0,\dots,n$ where $\mathbf{i}=\sqrt{-1}$. Using the identity $\sin(z)=\frac{1}{2\mathbf{i}}(e^{\mathbf{i}z}-e^{-\mathbf{i}z})$, we define the new complex variables

[TABLE]

These variables represent the phases of the oscillator as points on the complex plane. If $r_{i}=0$, then $x_{i}$ lies on the unit circle as the original formulation requires. If $r_{i}\neq 0$, $x_{i}$ deviates from the unit circle and no longer represents real solutions of the original equations. However, such extraneous solutions (i.e. non-real solutions) can be identified easily. With these new variables, the transcendental terms in (2) can be converted to rational functions

[TABLE]

and the system can be transformed into the algebraic synchronization system

[TABLE]

in the complex variables $x_{1},\dots,x_{n}$ where $a_{ij}^{\prime}=\frac{k_{ij}}{2\mathbf{i}}$. In this formulation, we also allow for complex coupling strength and natural frequencies. This algebraic formulation has been used in Ref. Chen, 2019; Chen, Davis, and Mehta, 2018, and it is similar to the formulation in Ref. Baillieul and Byrnes, 1982. However, the above formulation directly connects to the construction of “adjacency polytopes.”

Denote the rational functions on the left hand sides of the above system by $\mathbf{f}=(f_{1},\dots,f_{n})^{\top}$. It is clear that the two systems $\mathbf{f}=\mathbf{0}$ and $M\mathbf{f}=\mathbf{0}$ have the exact same solution set for any nonsingular $n\times n$ matrix $M=[M_{ij}]$. Therefore, without loss of generality, we can consider the equivalent system $M\mathbf{f}=\mathbf{0}$, which is of the form

[TABLE]

with $c_{k}=\sum_{i=1}^{n}M_{ki}\omega_{i}$. This system is referred to as the unmixed form of the algebraic synchronization system (simply unmixed synchronization system hereafter). Replacing a system of equations $\mathbf{f}=\mathbf{0}$ by $M\,\mathbf{f}=\mathbf{0}$ for a nonsingular square matrix $M$ is a common practice in numerical methods for solving nonlinear systems (e.g., the randomization procedure in numerical algebraic geometry Sommese and Wampler (2005)).

In this form, every monomial appears in every equation, and each equation no longer represents the balancing condition on a single oscillator. Instead, each equation becomes a linear combination of all balancing conditions in the network. The rest of the paper focuses on this system.

IV A example

We briefly outline the basic idea behind the proposed decomposition scheme with a simple example of a network of three coupled oscillators (Figure 1a). Though this case is simple enough to be solved directly by algebraic manipulations, it can still illuminate main idea of the proposed framework. In this case, the unmixed synchronization system (5) is

[TABLE]

Existing root count resultsBaillieul and Byrnes (1982); Chen, Davis, and Mehta (2018) dictate that there are no more than six complex synchronization configurations for this Kuramoto network of three oscillators. Our goal is to understand these synchronization configurations by examining simpler subnetworks of this network.

Figure 2 displays the proposed decomposition scheme of this Kuramoto network into six subnetworks, supported by directed acyclic graphs (i.e., asymmetric couplings) with each subnetwork corresponding to one of the possible complex synchronization configurations.

Each of the subnetworks has its own synchronization system. For instance, the subnetwork displayed in Figure 1b corresponds to the unmixed synchronization system given by

[TABLE]

which is a system of equations involving a subset of terms in (6). Here, each edge is interpreted as a directed edge, and the directed edge $(i,j)$ corresponds to the term $x_{i}/x_{j}$. It is easy to verify that for generic choices of coefficients, the above system has a unique complex solution. The same analysis can be applied to the rest of the subnetworks presented in Figure 2, and we can verify that the synchronization system induced by each subnetwork has a unique solution under the assumption of generic coefficients. Moreover, these solutions are in one-to-one correspondence with the six complex frequency synchronization configurations of the original network in Figure 1a via a homotopy function $H(x_{1},x_{2},t):\mathbb{C}^{2}\times\mathbb{R}\to\mathbb{C}^{2}$ given by

[TABLE]

At $t=1$, the equation $H(x_{1},x_{2},1)=0$ is equivalent to the original unmixed synchronization system (6). Under a mild “genericity” assumption, it can be proved that as $t$ varies continuously from 1 to 0 (without reaching 0), the complex solutions of $H(x_{1},x_{2},t)=0$ also move smoothly, forming smooth paths emanating from the solutions at $t=1$.

As $t$ approaches 0, $H$ becomes undefined, and the six solution paths degenerate into the complex frequency synchronizations configurations of the six subnetworks in Figure 1a in the sense that the six complex synchronization configurations of the six subnetworks are the limit points of the following reparametrization of the six paths:

[TABLE]

That is, as $t\to 0$, along each path, $(y_{1},y_{2})$ converge to the unique solution of the synchronization systems of one of the subnetworks shown in the decomposition Figure 2.

As we will demonstrate in section VI, the unique solution of each of the six subsystems (e.g. (7)) can be computed easily and explicitly. This gives rise to a practical algorithm that can be used to find all (complex) solutions to system  (6). After finding the solutions to the six subsystems, one can trace the smooth paths, starting from these solutions and reach the solutions to the original synchronization system (6).

The deformation induced by $H(x_{1},x_{2},t)$ works by moving the parameter $t$ between 1 and 0 with $t=1$ corresponding to the original system. This has the effect of amplifying the difference in natural frequencies of the oscillators. As $t\to 0$, the difference in natural frequencies goes to infinity, at which point synchronization of the entire network is no longer possible (even considering complex configurations) and the network breaks down into six simpler subnetworks.

One important feature of this decomposition scheme is that it is guaranteed to be able to find all complex solutions to (6), which include all real solutions to the original Kuramoto equations (2). In the following sections, we develop the general construction.

V Adjacency polytope and facet networks

The decomposition illustrated above is constructed from the geometric information encoded in a polytope — the “adjacency polytope.” Chen (2019); Chen, Davis, and Mehta (2018) In this section, we briefly review the definition. A polytope is the convex hull of a finite lists of points in $\mathbb{R}^{n}$ (a bounded geometric object with finitely many flat faces). An adjacency polytope is a polytope that can capture the topology of a network, count the number of synchronization configurations, and guide the decomposition of the network.

Definition 1** (ChenChen (2019)).**

Given a Kuramoto network $(G,K,\omega)$, the adjacency polytope of this network is the polytope

[TABLE]

and the adjacency polytope bound of this network is the normalized volume $\operatorname{NVol}(\nabla_{G})=n\,!\operatorname{vol}_{n}(\nabla_{G})$.

Here, $\mathbf{e}_{i}$ is the $i$-th standard basis of $\mathbb{R}^{n}$ and $\mathbf{e}_{0}=\mathbf{0}$. The adjacency polytope is the convex hull of entries in the (directed) adjacency list of the graph as points in $\mathbb{R}^{n}$, and it is a geometric encoding of the topology of the network. This construction has been studied in the context of Kuramoto model and power flow equations.Chen (2019); Chen and Mehta (2018) Similar constructions have also appeared in other contexts.Matsui et al. (2011); Higashitani, Kummer, and Michałek (2016); Delucchi and Hoessly (2016) The adjacency polytope bound is a simplification and relaxation of the Bernshtein-Kushnirenko-Khovanskii bound.Bernshtein (1975) Recently, the author proved that the this bound remains a sharp upper bound for the total number of complex synchronization configurations for certain graphs.Chen (2019) This bound is sharp in the sense that it is attainable when generic coupling strength and natural frequencies are used. The explicit formulae for this bound in the case of cycle and tree graphs have also been established.Chen, Davis, and Mehta (2018)

V.1 Facet subnetworks and subsystems

The main goal of this paper is to demonstrate that much more information of the original Kuramoto network can be extracted from this polytope. Of particular importance are the facets of the adjacency polytope, which are proper faces of maximum dimension. Each facet gives rise to a directed acyclic subnetwork of the original network. To define such generalized Kuramoto networks, we first need to reinterpret synchronization equations (4) and (5) from the view point of directed graphs. In (4), an (undirected) edge $\{i,j\}$ corresponds to the term $a_{ij}^{\prime}(x_{i}/x_{j}-x_{j}/x_{i})$. We can split this undirected edge into directed edges $(i,j)$ and $(j,i)$ with the corresponding terms $a_{ij}^{\prime}x_{i}/x_{j}$ and $a_{ij}^{\prime}x_{j}/x_{i}$ respectively, as illustrated in Figure 3. This allows us to define more general directed Kuramoto networks induced by special directed graphs:

Definition 2**.**

For a facet $F$ of $\nabla_{G}$, we define

[TABLE]

With these, we define the facet subnetwork associated with $F$ to be the (directed) Kuramoto network $(\mathcal{V}_{F},\mathcal{E}_{F})$ and the corresponding facet subsystem to be the system of equations

[TABLE]

A facet subsystem is a system of equations involving a subset of terms in the unmixed synchronization system (5) of the original network. The selection of the terms depends on the nodes and edges that appear in the corresponding facet subnetwork. In this context, we consider the edge $(i,j)$ to be a directed edge. That is, $(i,j)\in\mathcal{E}_{F}$ does not imply $(j,i)\in\mathcal{E}_{F}$. As we shall prove, the directed edges $(i,j)$ and $(j,i)$ can never be in the same facet subnetwork.

Remark 1**.**

The splitting of an undirected edge into two directed edges displayed in Figure 3 can also be interpreted as an explicit use of an imaginary interaction term. In considering the real frequency synchronization configurations defined by (2), we require $|x_{i}|=1$ for all $i$ since $x_{i}=e^{r_{i}+\mathbf{i}\theta_{i}}$. In this case,

[TABLE]

which can be interpreted as the interaction term along the directed edge $(i,j)$. If both $(i,j)$ and $(j,i)$ are present, then the imaginary parts cancel each other out, leaving only the real part $k_{ij}\sin(\theta_{i}-\theta_{j})$ as the combined interaction term that matches the original synchronization equations (2). This interpretation of the unidirectional interaction is similar to yet different from the interpretations used in the recent studies about Kuramoto models on directed graphs. Delabays, Jacquod, and Dörfler (2018); Rogge and Aeyels (2004)

Remark 2** (Facet computation).**

In this section, we have established an one-to-one correspondence between the facets of an adjacency polytope and the facet subnetworks of a Kuramoto network. As we demonstrate in the following subsections, a great deal of information about all possible synchronization configurations can be extracted from the facet networks. The proposed analysis therefore hinges on finding all facets of adjacency polytopes. Fortunately, enumerating all facets of a convex polytope in high dimensions is a well-studied problem in convex geometry, and there exist several software packages that can carrying out such calculations efficiently.Fukuda, Liebling, and Margot (1997); Büeler, Enge, and Fukuda (2000); Avis and Jordan (2018); Gawrilow and Joswig (2000) Although this problem is known to be of class #P in general, given the symmetries built into the definition of adjacency polytopes, it is plausible that facets of this type of polytope can be listed much more efficiently. For instance, the theory for exploiting the symmetry with respect to the origin ($\mathbf{x}\in\nabla_{G}$ if and only if $-\mathbf{x}\in\nabla_{G}$) is well understood.

V.2 Topological properties of facet subnetworks

A natural question to ask here is as follows: what special topological properties do facet subnetworks posses? We now establish property restrictions and symmetries.

Theorem 1**.**

Given a Kuramoto network, the facet subnetworks (definition 2) have the following properties:

• •

Each facet subnetwork is acyclic;

• •

Each facet subnetwork contains all the nodes;

• •

Each facet subnetwork is weakly connected;

• •

Any two paths between a pair of nodes in facet subnetwork are of the same length;

• •

Each path in a subnetwork contains no more than half of the edges of any cycle in $G$ containing it; and

• •

For each facet subnetwork, its transpose is also a facet subnetwork.

Proof.

Let $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{R}^{n}$ be the inner normal vector of a facet $F$ of $\nabla_{G}$, then there exists an $h\in\mathbb{R}$ such that

[TABLE]

However, $\mathbf{0}$ is an interior point of $\nabla_{G}$. Therefore $h<\langle\mathbf{0},\alpha\rangle=0$.

(Acyclic) Suppose a facet subnetwork $G_{F}$ contains a directed cycle $i_{1}\to i_{2}\to\cdots\to i_{m}\to i_{1}$ for some $\{i_{1},\dots,i_{m}\}\subset\{0,\dots,n\}$. Let $i_{m+1}=i_{1}$, then $\mathbf{e}_{i_{k}}-\mathbf{e}_{i_{k+1}}\in F$ for $k=1,\dots,m$. Since $\bm{\alpha}$ is the inner normal of $F$, we have

[TABLE]

and

[TABLE]

which is a contradiction.

(Containing all nodes) For $n=1$, there are exactly two directed edges corresponding to the two facets of $\nabla_{G}$ which is simply a line segment. The statement clearly holds.

For $n>1$, since $\dim\nabla_{G}=n$, the facet $F$ is $(n-1)$-dimensional. Suppose there is a node $i\not\in G_{F}$. Without loss of generality, we can assume $i\neq 0$ by reassigning the reference node. Let $G^{\prime}$ be the graph obtained from $G$ by removing the node $i$ and all its edges. Then, $\nabla_{G^{\prime}}\subset\nabla_{G}$, and the set of vertices of $\nabla_{G^{\prime}}$ is a subset of the vertices of $\nabla_{G}$. It is straightforward to verify that $F$ remains a face of $\nabla_{G^{\prime}}$. Therefore, $\dim F<\dim\nabla_{G^{\prime}}\leq n-1$, which contradicts with the assumption that $\dim F=n-1$.

(Weakly connected) Suppose $G_{F}$ is not weakly connected. Let $C_{1}\ni 0$ be the nodes in one weakly connected component and let $C_{2}=\mathcal{V}_{F}\setminus C_{1}$. By relabeling the nodes, we can assume $C_{1}=\{0,1,\dots,m\}$ and $C_{2}=\{m+1,\dots,n\}$. Define

[TABLE]

Since $C_{1}$ and $C_{2}$ belong to different weakly connected components of $G_{F}$, i.e., there is no $\pm(\mathbf{e}_{i}-\mathbf{e}_{j})\in F$ with $i\in C_{1}$ and $j\in C_{2}$, we can see that $V_{1}\cap V_{2}=\varnothing$.

By construction, $V_{1}\subset\mathbb{R}^{m}\times\{\mathbf{0}\}$ and $V_{2}\subset\{\mathbf{0}\}\times\mathbb{R}^{n-m}$. Since points in $V_{2}$ are points of the form $\mathbf{e}_{i}-\mathbf{e}_{j}$, we can see that they actually belong to the smaller subspace defined by $\langle\;\cdot\;,\mathbf{1}\rangle=0$. Therefore, treating points in $V_{1}$ and $V_{2}$ as vectors, we get $\dim(\operatorname{span}(V_{1})\leq m)$ and $\dim(\operatorname{span}(V_{2})<n-m)$. Consequently, the pyramid formed by $F=\operatorname{conv}(V_{1}\cup V_{2})$ and $\mathbf{0}$ is of a dimension strictly less than $m+n-m=n$, and therefore, the dimension of the facet $F$ itself must be strictly less than $n-1$. This contradicts with the fact that $\nabla_{G}$ is full-dimensional and $F$ is one of its facets.

(Equal length) Suppose there are two paths $i_{0}\to i_{2}\to\dots\to i_{m}$ and $j_{0}\to j_{2}\to\dots\to j_{\ell}$ with the same starting and end points, that is, $i_{0}=j_{0}$ and $i_{m}=j_{\ell}$. Then

[TABLE]

Summing these $m$ equations, we obtain

[TABLE]

Similarly, summing the $\ell$ equations along the path $j_{0}\to j_{2}\to\dots\to j_{\ell}$ produces

[TABLE]

However, $j_{0}=i_{0}$ and $i_{m}=j_{\ell}$. Thus we must have $mh=\ell h$, i.e., $m=\ell$, and these two paths must have the same length.

(Half cycle) Suppose $G$ contains an undirected cycle $i_{0}\leftrightarrow\dots\leftrightarrow i_{m}\leftrightarrow\dots\leftrightarrow i_{\ell}=i_{0}$. Among these edges, the path $i_{0}\to i_{2}\to\dots\to i_{m}$ is in $G_{F}$ we can show that $m$ is no more than $\ell/2$. By construction, $\mathbf{e}_{i_{k}}-\mathbf{e}_{i_{k+1}}\in F$ for each $k=0,\dots,m-1$, and hence,

[TABLE]

Summing these $m$ equations, we obtain

[TABLE]

Similarly, consider the path $i_{\ell}\to i_{\ell-1}\to\cdots\to i_{m}$, which is in $G$. We must have

[TABLE]

Summing these equations produces

[TABLE]

Recall that $\left\langle\mathbf{e}_{i_{1}}-\mathbf{e}_{i_{m}}\,,\,\bm{\alpha}\right\rangle=mh$. Thus we must have $mh\geq(\ell-m)h$. Since $h<0$, this implies that $m\leq\ell-m$, i.e., $2m\leq\ell$.

Since $\nabla_{G}$ is centrally symmetric, $-F=\{-\mathbf{x}\mid\mathbf{x}\in F\}\subset\nabla_{G}$, and $\dim(-F)=\dim F=n-1$. Moreover, for $\mathbf{e}_{j}-\mathbf{e}_{i}\in-F$, $\mathbf{e}_{i}-\mathbf{e}_{j}\in F$ by assumption, thus

[TABLE]

with the inner normal $\bm{\alpha}$ of the facet $F$. Similarly,

[TABLE]

for any $-\mathbf{p}\in\nabla_{G}\setminus(-F)$. Therefore we can conclude that $-F$ is also a facet of $\nabla_{G}$. We can verify that the facet network $G_{-F}$ is the transpose of $G_{F}$ by definition. ∎

V.3 Facet subnetworks as maximal flow networks

The facet subnetworks can also be understood from the point of view of flow networks. Indeed, they form subnetworks consisting of edges on which certain pseudo-flow is maximized.

Given a vector $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{R}^{n}$, we can define the pseudo-flow $f:V\times V\to\mathbb{R}$ given by

[TABLE]

with $\alpha_{0}=0$. Without loss of generality, we further require that the maximum flow be 1, which can be used as the capacity for all edges. By combining this setup with the topological properties established in Theorem 1, we can see that according to this interpretation, a facet subnetwork is the subnetwork consisting of all vertices $\{0,1,\dots,n\}$, but only the edges with no residual capacity, i.e., edges $(i,j)$ with $f_{\bm{\alpha}}(i,j)=1$. Simply put, a facet subnetwork is the largest subnetwork on which the pseudo-flow $f_{\bm{\alpha}}$ is maximized. The problem of listing all facet subnetworks is therefore equivalent to the problem of finding all possible assignments of $\bm{\alpha}$ that induce such a maximized flow-subnetworks. Classic flow network algorithms based on linear programming problems can thus be used here.

V.4 Algebraic properties of facet subsystems

By adopting the complex formulation (4) of the synchronization system, we paid the price of potentially introducing extraneous solutions (synchronization with imaginary phase angles). However, this effort allows us to employ more powerful tools from complex algebraic geometry. A fundamental fact in complex algebraic geometry is that, as far as complex root count is concerned, the generic behavior of a family of algebraic systems coincides with the maximal behavior. That is, if we consider the facet system as a family of algebraic systems parameterized by the coefficients, then the generic complex root count coincides with the maximum root count. This is an important consequence of the Bertini’s Theorem, and it forms the foundation of much of numerical algebraic geometry (see Ref. Sommese and Wampler, 2005, Theorem 7.1.1). Applying this to our context, we obtain the following proposition.

Proposition 1**.**

*For a generic choice of coefficients $\{c_{k}\}$ and $\{a_{ijk}\}$, the total number of isolated111 Isolated solutions here refer to geometrically isolated solutions. A solution of a system of equations is said to be geometrically isolated if there is a nonempty open set in which it is the only solution.

complex solutions to a given facet subsystem (9) induced by the facet $F$ of an adjacency polytope is a constant. Denote this constant by $\mathcal{N}(F)$, and then $\mathcal{N}(F)$ is also the maximum number of isolated complex solutions this facet subsystem could have.*

More specifically, this proposition can be understood as a direct application of the cheater’s homotopy Li, Sauer, and Yorke (1989) or the parameter homotopy Morgan and Sommese (1989) theory to the facet subsystems. Moreover, Kushnirenko’s Theorem Kouchnirenko (1976) provides us with the explicit formula for the maximum root count in the form of normalized volume.

Proposition 2**.**

For a facet $F$ of an adjacency polytope $\nabla_{G}$,

[TABLE]

Proof.

Since the directed graph $(\mathcal{V}_{F},\mathcal{E}_{F})$ associated with a facet subnetwork is necessarily acyclic by theorem 1, edges $(i,j)$ and $(j,i)$ cannot both be in $\mathcal{E}_{F}$. Consequently, the coefficients in the facet subsystem (9) are independent. That is, generic choices of the parameters $\{c_{k}\}$ and $\{a_{ijk}\}$ for the network correspond to independent generic choices of the coefficients for the facet subsystem.

In addition, it is important to note that the support of the facet subsystem — the collection of all the exponent vectors appeared in the system — is the set of vertices in $F$ together with $\mathbf{0}$. According to Kushnirenko’s Theorem Kouchnirenko (1976), $\mathcal{N}(F)=\operatorname{NVol}(\operatorname{conv}(F\cup\{\mathbf{0}\}))$. ∎

$\mathcal{N}(F)$ can be considered a generalization of the adjacency polytope bound to facet subnetworks. Moreover, the sum of the root count for all facet subsystems gives us a root count for the whole system. It is critical to note here that the factor of $n!$ in $\mathcal{N}(F)$ above is a part of the definition of normalized volume and it is used to compensate for the shrinking of the Euclidean volume itself (the volume of a minimum simplex with integral vertices is $1/n!$ and goes to 0 as $n\to\infty$). It does not indicate factorial growth in the root count for facet subsystems. As section VI will demonstrate, the smallest facet subsystem has exactly one solution regardless of the dimension.

Theorem 2**.**

Given a Kuramoto network, let $\nabla_{G}$ be the adjacency polytope defined above. Then the total number of isolated complex synchronization configurations this network has is bounded above by the sum of generic complex root counts of the facet subsystems, that is, it is bounded by

[TABLE]

where $\mathcal{F}(\nabla_{G})$ is the set of facets of the polytope $\nabla_{G}$.

Proof.

$\mathbf{0}$ is an interior point of $\nabla_{G}$. Thus the collection of pyramids

[TABLE]

form a subdivision of $\nabla_{G}$. The normalized volume of $\nabla_{G}$ is the sum of the normalized volume of these pyramid. That is,

[TABLE]

When this result is combined with theorem 1, it is plausible that the adjacency polytope bound $\operatorname{NVol}(\nabla_{G})$ can be computed simply by listing all possible facet subnetworks and computing the adjacency polytope bound for each of such simpler subnetworks. This process is potentially easier than the $\#P$ problem of volume computation in general. As we demonstrate in the examples in section VII, this can be done for certain classes of networks.

V.5 Degeneration into facet subsystems

The true value of the facet subsystems lies in their roles as destinations for a deformation of the synchronization equations. That is, we can form a continuous deformation (in the sense of homotopy) of the unmixed synchronization equations that can degenerate into simpler facet subsystems and reduce the problem of finding synchronization configurations to the problem of solving each individual facet subsystem. This can be done through a specialized homotopy.

Consider the homotopy $H(\mathbf{x},t)=(H_{1},\dots,H_{k})$ given by

[TABLE]

The equation $H(\mathbf{x},1)=\mathbf{0}$ is exactly the unmixed synchronization system (5). Since $H(\mathbf{x},t)$ is smooth in $t$ for $t\in(0,1]$, as $t$ varies between 0 and 1, the equation $H(\mathbf{x},t)=\mathbf{0}$ represents a continuous deformation of (5). With proper choices of coefficients, the solutions of $H(\mathbf{x},t)=\mathbf{0}$ also move smoothly, as $t$ varies within $(0,1]$, forming smooth solution paths. As $t\to 0$, however, this smooth deformation breaks down, and the equation $H(\mathbf{x},t)=\mathbf{0}$ degenerates into the facet subsystems.

Theorem 3**.**

Given a Kuramoto network $(G,K,\omega)$ with generic coupling strengths and natural frequencies, the solution set of the system of equations $H(\mathbf{x},t)=\mathbf{0}$ contains a finite number of smooth curves parameterized by $t$ such that:

The set of limit points of these curves as $t\to 1$ contains all isolated complex synchronization configurations of this network; and 2. 2.

The set of limit points of these curves as $t\to 0$ reaches all complex synchronization configurations of the facet subsystems at “toric infinity” in the sense that for each of these complex solutions $\mathbf{y}=(y_{1},\dots,y_{n})$, there exists a curve defined by $H(\mathbf{x},t)=\mathbf{0}$, along which $x_{i}(t)=y_{i}t^{\alpha_{i}}+o(t)$ for some $\alpha_{i}\in\mathbb{Q}$ and $t$ sufficiently close to 0.

In this case, “generic” coupling strength and natural frequencies can be understood as coefficients with arbitrarily small perturbation. Given any set of $K=[k_{ij}]$, $\omega=(\omega_{1},\dots,\omega_{n})$, and a threshold $\epsilon>0$, there exists a new set of coefficients no more than $\epsilon$-distance away from $(K,\omega)$ for which the above statement holds true. This can also be interpreted from a probabilistic point of view: if $K$ and $\omega$ are selected at random, then the above statement holds true with probability one.

“Toric infinity” is a term from complex algebraic geometry that is used to describe space outside $(\mathbb{C}\setminus\{0\})^{n}$ but inside a certain closure. In this context, as $t\to 0$, some coordinates in $\mathbf{x}=(x_{1},\dots,x_{n})$ either goes to 0 or infinity. In either case, the point $\mathbf{x}$ escapes $(\mathbb{C}\setminus\{0\})^{n}$, known as the “algebraic torus”.

Proof.

(1) This statement describes a direct consequence of the parameter homotopy method Morgan and Sommese (1989). If we consider the unmixed synchronization system (5) as a family $\mathbf{f}(\mathbf{x};\mathbf{c})$ parameterized by constant terms $\mathbf{c}=(c_{1},\dots,c_{n})$, then the parameter $\mathbf{c}/t$ remains generic for generic choices of $\mathbf{c}$ and every $t\in(0,1]$. According to the theory of parameter homotopy Morgan and Sommese (1989), the solution set of $H(\mathbf{x},t)=\mathbf{0}$ in $\mathbb{C}^{n}\times(0,1]$ consists of smooth paths parameterized by $t$ whose limit points as $t\to 1$ include all solutions of $H(\mathbf{x},1)\equiv\mathbf{f}(\mathbf{x};\mathbf{c})=\mathbf{0}$.

(2) Let $C\subset\mathbb{C}^{n}\times(0,1]$ be a curve defined by $H(\mathbf{x},t)=\mathbf{0}$. By the previous part, $C$ can be expressed as a path $\mathbf{x}(t)$, smoothly and analytically parameterized by $t$ for $t\neq 0$. Although the smoothness breaks down at $t=0$, according to the theory of polyhedral homotopy method Huber and Sturmfels (1995) as well as the Newton-Puiseux Theorem, there exists a smooth function $\mathbf{y}(t)=(y_{1},\dots,y_{n})$ and a vector $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{Q}^{n}$ such that $\mathbf{y}(t)$ has a limit point in $\mathbb{C}^{n}$ as $t\to 0^{+}$ and

[TABLE]

for sufficiently small positive $t$ values. Moreover, $\hat{\bm{\alpha}}=(\bm{\alpha},1)\in\mathbb{Q}^{n+1}$ must be an inner normal vector of a facet of the Newton polytope $\operatorname{Newt}(H)$ of $H(\mathbf{x},t)$ which is the pyramid formed by $\nabla_{G}\subset\mathbb{R}^{n}\subset\mathbb{R}^{n+1}$ together with $(\mathbf{0},-1)\in\mathbb{R}^{n+1}$, i.e., $\hat{\bm{\alpha}}$ is an inner normal vector of

[TABLE]

In this case, $\bm{\alpha}$ is an inner normal vector of a facet $F$ of $\nabla_{G}$. Consequently, the initial form $\operatorname{init}_{\alpha}(H)$ is exactly the facet subsystem induced by $F$ as displayed in (9), and $\lim_{t\to 0^{+}}\mathbf{y}$ exists and is a solution to this system. ∎

The above theorem reduces the problem of solving a complicated synchronization system (5) to the problem of solving simpler facet subsystems: Once the solutions to each facet subsystem are found, they become the starting points of the smooth paths defined by $H(\mathbf{x},t)=\mathbf{0}$. Then, standard numerical “path tracking” algorithms Allgower (1981); Chen and Li (2015); Davidenko (1953); Sommese and Wampler (2005) can be applied to trace these smooth paths and reach the desired synchronization configurations defined by $H(\mathbf{x},1)=\mathbf{f}(\mathbf{x})=\mathbf{0}$. This homotopy method has the advantage of being pleasantly parallel in the sense that every configuration can be found independently.

Remark 3** (Connection to polyhedral homotopy).**

*The homotopy constructed from the above theorem can be viewed as a highly specialized polyhedral homotopy method Huber and Sturmfels (1995) with a special lifting function that takes the value of $-1$ on constant terms and 0 everywhere else. Moreover, instead of degenerating into binomial systems that depend on lifting, our homotopy degenerates into facet subsystems that may be found through an examination of the network topology. *

Remark 4** (Pushing Kuramoto networks to the limit).**

As noted in the end of section VII, the degeneration of a given network into facet subnetworks can be thought of as pushing the network to its breaking point. Since the construction of the homotopy $H$ replaces constant terms in (5) by $c_{k}/t$, decreasing $t$ from 1 to 0 has the effect of amplifying the differences in the natural frequencies. We expect this to make synchronization more difficult resulting in some real synchronization configurations degenerating into non-real complex configurations. As $t\to 0$, the differences in the natural frequencies are amplified to infinity, at which point the network is not even able to support complex synchronization configurations. The network breaks apart, and the facet subnetworks describe the pieces in the sense of limits.

VI Primitive facet subnetworks

The general facet decomposition scheme outlined above decomposes a Kuramoto network into a collection of smaller directed acyclic networks involving the original oscillators — the facet subnetworks. Among these subnetworks, the basic building blocks are “primitive” subnetworks which are, in a sense, the smallest Kuramoto networks that could have a frequency synchronization configuration.

Definition 3**.**

In a given Kuramoto network, a facet subnetwork is called primitive if the underlying graph $G$ is a weakly connected acyclic graph that contains all the vertices ${0,\dots,n}$ and exactly $n=N-1$ directed edges.

Since such a primitive subnetwork contains exactly $N$ node and $N-1$ directed edges, removing any edge creates disconnected components. It is thus a minimum facet subnetwork. Yet, as we demonstrate in this section, it is sophisticated enough to have a complex synchronization configuration, and this configuration must be unique. Equally as important, this unique synchronization configuration can be computed quickly and easily using only $O(n)$ complex multiplication and division, and no additional memory is needed. Once solved, the solutions of these primitive facet subsystems can be used to bootstrap the homotopy continuation procedure described in the previous section.

Theorem 4**.**

For $N>1$, the facet subsystem (9) supported by a primitive facet subnetwork has a unique complex solution for which $x_{i}\neq 0$ for each $i=1,\dots,n$.

Proof.

Since $G$ has exactly $n$ directed edges, the corresponding synchronization system has exactly $n$ nonconstant terms in each equation. Under the assumption of generic coefficients, we can reduce the system to

[TABLE]

via Gaussian elimination with no cancellation of the terms with $\{c_{k}^{\prime}\}$ and $\{a_{k}^{\prime}\}$ representing the resulting coefficients. The invertibility of the Gaussian elimination process ensures that all $a_{k}^{\prime}\neq 0$. Having exactly $n=N-1$ edges and $N$ nodes also ensures that the underlying undirected graph is an undirected tree. Therefore, for each $i=1,\dots,n$, there is a unique path between node [math] and node $i$ through the nodes $i_{0},i_{1},\dots,i_{m}$ with $i_{0}=0$ and $i_{m}=i$, where $m$ is the length of this path.

[TABLE]

We should recall that node [math] is the reference node, i.e., $x_{i_{0}}=x_{0}=1$. Solving this system via a forward substitution process, we can determine the value of $x_{i_{m}}=x_{i}\neq 0$. Since the path between node $i$ and node [math] is unique, there is no possibility of inconsistency in this process, and hence the value of each $x_{i}$ is uniquely determined. We can therefore conclude that a solution to the facet system exists, and it is unique. ∎

Aside from offering starting points for the homotopy continuation method that will identify all complex synchronization configurations, the existence of primitive facet subnetworks also directly contributes to an estimation of the total number of complex synchronization configurations. Since each primitive facet subsystem has a unique solution, by  theorem 2, the total number of complex synchronization configurations is greater than or equal to the number of primitive facet subnetworks.

Remark 5** (Computational complexity).**

From a computational view point, the true value of primitive subnetworks lies in the ease with which it can be solved. The above proof is constructive and suggests a practical algorithm for solving the reduced facet subsystem (12) that only requires $O(n)$ complex multiplication and division, and no additional memory once the initial elimination step is done. This is of great importance in the homotopy method for solving the synchronization system described in the previous section, since the solutions to the facet subsystems are the starting points of solution paths that can lead to the desired synchronization configurations of the original Kuramoto network.

Remark 6** (Toric interpretation).**

The above theorem can also be interpreted in the language of toric algebraic geometry as follows. The facet subsystem supported by a primitive facet subnetwork can be reduced to an equivalent binomial system whose associated exponent matrix (the matrix whose columns are the exponent vectors for the nonconstant terms) is a unimodular matrix. Consequently, it defines a (normal) irreducible toric variety of zero dimension which must be a single point. From this toric geometry view point, we can also see that primitive subnetworks are the smallest building blocks in the proposed decomposition scheme.

Remark 7** (Real synchronization configurations).**

Recall that algebraic synchronization equations are obtained through a change of variables $x_{i}=e^{\mathbf{i}\theta_{i}}$ for $i=1,\dots,n$, so real solutions to the original Kuramoto equations corresponds to complex solutions to the algebraic system with $|x_{i}|=1$ for $i=1,\dots,n$ (i.e., the real torus solutions in $(S^{1})^{n}=(\{z\in\mathbb{C}\mid|z|=1\})^{n}$). Since a primitive facet subsystem (12) can be solved using only complex multiplication and division, we can ensure $|x_{i}|=1$ for $i=1,\dots,n$ by choosing coefficients that have unit absolute values. In this case, the unique solution corresponds to a real solution. That is, there is always a way to choose the coefficients so that the unique solution of a given primitive facet subsystem gives rise to a real solution. Under the interpretation of facet subnetworks as acyclic directed Kuramoto networks with one-way interaction (fig. 3), this real solution can be viewed as a generalized real synchronization configuration that represents the limit behavior of the Kuramoto network as the differences in the natural frequencies are amplified to infinity.

As the next section demonstrates, primitive subnetworks appear naturally in the facet decomposition of many types of Kuramoto networks, e.g., trees, cycles, and chordal graphs. In some cases, all facet subnetworks are primitive.

VII Examples

In this section, we illustrate the facet decomposition scheme using concrete examples. In all examples, the facets of the adjacency polytopes are computed using the well test software package Polymake,Gawrilow and Joswig (2000) which can compute the list of all facets efficiently even for complicated convex polytopes in very high dimension.

VII.1 A tree network

We start with the simplest types of networks — trees. For a tree graph containing five nodes, displayed in Figure 4a, the corresponding adjacency polytope (definition 1) has 16 facets which produces 16 facet subnetworks, displayed in Figure 5. Every facet subnetwork is primitive in this case. Since each primitive facet subnetwork has a unique complex synchronization configuration, we can conclude that the original tree network has at most 16 complex synchronization configurations. This aligns with the well known result that on a tree network of $N$ oscillators, there are at most $2^{N-1}$ complex synchronization configurations. Adjacency polytope bounds for trees networks in general were analyzed in a closely related workChen, Davis, and Mehta (2018) using convex geometry method. With the framework proposed in this paper, we can see that this bound actually has a topological origin — an acyclic cover of a tree network. In this case, the topological constraints given in theorem 1 actually determine the list of facet subnetworks. That is, without actually computing the facets of the adjacency polytope, one can easily enumerate all of the facet subnetworks by listing the acyclic subnetworks.

Interestingly, this upper bound on the number of complex synchronization configurations coincide with the tightest possible upper bound on the number of real configurations as it is well known that there could also be as many as $2^{N-1}$ real synchronization configurations. Baillieul and Byrnes (1982)

VII.2 A cycle network

For a cycle network, there is a similar decomposition, but with more subnetworks. This is consistent with the observation that cycle networks could have more possible synchronization configurations than tree networks. Unlike the previous case, not all acyclic subgraphs give rise to facet subnetworks. Figure 6 displays all 30 of the facet subnetworks from a cycle graph of five nodes (fig. 4b). Each subnetwork is primitive and hence has a unique complex synchronization configuration. Consequently, a cycle network of five oscillators has at most 30 complex synchronization configurations. This number also coincide with the maximum number of real synchronization configurations obtained by a recent study,Zachariah et al. (2018) i.e., all complex solutions could be real.

The obvious symmetry shown in Figure 6 is not accidental. In addition to the transpose symmetries stated in Theorem 1 (for a facet subnetwork $G_{F}$, its transpose $G_{F}^{\top}$ is also a facet subnetwork), we can also see that the set of facet subnetworks inherited the full range of cyclic symmetry from the original network.

Remark 8**.**

This computationally obtained root count for the Kuramoto equations induced by cycle networks has also been rigorously analyzed using techniques from convex geometry in a related work,Chen, Davis, and Mehta (2018) where an explicit formula for the root count is established. Moreover, a follow-up paper Chen and Davis (2018) by Robert Davis and the author demonstrated that for cycle networks with an odd number of oscillators, the facet decomposition scheme always produces the most refined decomposition into primitive subnetworks, and there is an one-to-one correspondence between these primitive subnetworks and the complex synchronization configurations of the original network under the genericity conditions on the coefficients. For cycle networks with an even number of oscillators, however, the facet decomposition scheme alone is not sufficient to produce only primitive subnetworks. A refined decomposition scheme Chen and Davis (2018) is proposed that can further decompose cycle networks into only primitive subnetworks.

VII.3 A chordal network

We now consider a chordal network of four oscillators, which consists of a cycle of four nodes together with a “chord” edge, as displayed in Figure 4c. Alternatively, such a graph can also be considered as two cycle graphs sharing common edges. Figure 7 shows all 12 of the facet subnetworks. Eight of them are primitive. The remaining four are non-primitive but less complicated than the original network. In total, this chordal network of four oscillators has 16 complex synchronization configurations: The eight primitive facet subnetworks each contribute one configuration, and the four non-primitive subnetworks each contribute two.

As in the previous case, the set of facet subnetworks exhibits symmetric structures. In additional to the transpose symmetry stated in Theorem 1, the set inherits from the original network the reflection symmetries given by $\{0\mapsto 2\}$, $\{1\mapsto 3\}$, etc.

VII.4 A wheel network

Figure 8 displays a “wheel network” of 10 oscillators which is simply a cycle network consisting of 9 oscillators together with an additional central node that is directly connected to all other nodes. This network has a total of 1598 facet subnetworks. 116 of them are primitive facet subnetworks (see for example Figure 9a and Figure 9b). The number of edges in the subnetworks range from 9 (primitive) to 13. Figure 9 includes some examples from each of these classes, and Table 1 summarizes the distribution of these subnetworks by size.

Transpose symmetry can be seen from fig. 9a and fig. 9b — one could be obtained from another by reversing the direction of all the edges. Similar to the case of cycle networks, the facet subnetworks exhibit cyclic symmetries inherited from the original network.

By computing the volume of each facet of $\nabla_{G}$, the contribution of each facet subnetwork to the total root count can be obtained, and they are summarized in table 1. For instance, each of the 116 primitive subnetworks contribute one complex root while each of the 414 non-primitive subnetworks having 10 edges contribute two complex roots. In total, the maximum number of complex synchronization configurations on this network is 8480, and this bound is exact for generic choices of coefficients.

In spite of having a large number of potential complex synchronization configurations, a modern homotopy continuation solver can still find all configurations quickly. For instance, a modified version222Hom4PS-3 implements the general polyhedral homotopy method.Huber and Sturmfels (1995) As noted in Remark 3, the proposed homotopy method can be considered as a highly specialized version of the unmixed form of the polyhedral homotopy method. The preliminary implementation used here combines Hom4PS-3 with a pre-processor that analyzes the adjacency polytope, generates suitable facet systems as starting systems, and bootstraps the homotopy method. This pre-processor relies on facet information computed by Polymake.Gawrilow and Joswig (2000) of Hom4PS-3 can solve the corresponding synchronization system in seconds on a moderate workstation. More importantly, the homotopy method scales almost linearly as more processor cores are used. Table 2 displays the computation time consumed on a workstation equipped with 16 GB of memory and Intel Xeon E5-2690 processor running at 2.9 GHz and up to 8 cores.

VIII Conclusion and future directions

This paper focuses on the study of frequency synchronization configurations in Kuramoto models for networks of coupled oscillators. It lays the foundation for a divide-and-conquer approach to analyzing large and complex Kuramoto networks from the view point of algebraic geometry. The main contribution of this paper is a general decomposition scheme for Kuramoto networks that can reduce a complicated Kuramoto network into a collection of simpler directed acyclic subnetworks which are generalized Kuramoto networks that allow one-way interactions among oscillators (c.f. Ref.Delabays, Jacquod, and Dörfler, 2018). The challenge of finding all possible frequency synchronization configurations is thus reduced to the problem of fully understanding these much simpler subnetworks.

Starting from a complex algebraic re-formulation of the underlying transcendental equations, the proposed framework uses the geometric information extracted from the “adjacency polytope” of the network, which is a polytope that encodes the network topology. This polytope is also the convex hull of the union of the Newton polytopes derived from the algebraic synchronization equations. The subnetworks are then in one-to-one correspondence with the facets (maximal faces) of the adjacency polytope. Since the adjacency polytopes are highly symmetric, the facets can be enumerated efficiently. Associated with each of these facet subnetworks is a much simpler, generalized synchronization equation, a “facet subsystem,” which only involves a fraction of the terms from the original algebraic synchronization equations. These subsystems are expected to be easier to solve than the original system. The facet subsystems associated with primitive subnetworks can be solved in linear time with no additional memory required.

From a computational point of view, the proposed decomposition framework gives rise to a homotopy continuation method. It induces a continuous deformation of the algebraic synchronization system that can degenerate into a collection of simpler facet subsystems. Homotopy continuation methods have been proved to be highly scalable as each solution can be computed independently. The can therefore be parallelized easily on modern high performance computing hardware.

The resulting homotopy method can be viewed as a highly specialized polyhedral homotopy method with fixed lifting functions (valuation functions). In the case that all subnetworks are primitive, this homotopy construction circumvents two of the main computational bottlenecks of the polyhedral homotopy method:

Since the facet subnetworks corresponds to facets of a well-known family of polytopes, they can be listed more easily, thereby allowing one to skip the costly step of “mixed cells computation”. 2. 2.

The starting system — facet subsystems — can be solved in linear time and does not require the costly step of computing Smith Normal Forms or Hermite Normal Forms of exponent matrices.

The examples in the section VII demonstrate that for certain types of networks, some facet subnetworks may not be primitive. The next step is to develop algorithms in order to refine our facet-based decomposition efficiently so that all resulting subnetworks are primitive. Recent fellow-up works have presented promising developments in this direction. For instance, in a recent paper by Robert Davis and the author, the directed acyclic decomposition scheme proposed here is significantly refined for cycle graphs so that all resulting subnetworks are primitive. This refinement is equivalent to a triangulation process for the facets of the adjacency polytope, but with the derivation of the explicit formula, the complexity of the refinement step for each facet is only linear in the number of oscillators. The refinement scheme for other types of networks remains an open problem.

In this work, we have expanded our domain to include complex synchronization configurations in order to take advantage of the powerful tools of complex algebraic geometry. These complex configurations can be of great value even when only real configurations are needed. First, the complex configurations necessarily include all real configurations. Through an examination of the imaginary parts of the phase angles, the extraneous (non-real) solutions can be filtered out easily. Similarly, by examining the eigenvalues of the Jacobian matrices at real solutions unstable solutions can also be filtered out if only stable solutions are desired. More importantly, the collection of complex configurations forms a reservoir of real configurations as one vary the parameters (natural frequencies and coupling strength). As the parameters move in the space of all possible parameters, for almost all parameters, the total number of complex configurations (counting multiplicity) remains the same, while the real configurations can increase or decrease. Non-real configurations may collide and form real configurations, and real configurations may degenerate into non-real ones. This phenomenon is a central question in numerical algebraic geometry and algebraic geometry in general. Numerical algorithms using all complex solutions as seeds for the exploration of the landscape of real solutions have been developed in the context of mechanical engineering.Griffin and Hauenstein (2015) Similarly numerical algorithms targeting real solutions of certain stability types have been studied in the context of chemical clusters.Chen and Mehta (2015) Combining these techniques and the framework developed in this work can shed new light on the possible real or even stable synchronization configurations of Kuramoto models.

Acknowledgments

This work was funded in part by the AMS-Simons Travel Grant and Auburn University at Montgomery Grant-In-Aid program and NSF under Grant Nos. 1115587 and 1923099.

The author would like to thank Dhagash Mehta for introducing the interesting topic of the Kuramoto model; Daniel Molzahn for the helpful discussions on the closely related power-flow equations; and Wuwei Lin for asking the question of whether or not a Kuramoto network can be decomposed into smaller subnetwork in a reasonable way — a question that started this research. The author also greatly benefited from discussions with Robert Davis. Indeed, discussions with Robert Davis produced a significant refinement of the framework developed in this paper for cycle graphs, and this extension is discussed in a separate follow-up paper Chen and Davis (2018).

*This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in Chaos: An Interdisciplinary Journal of Nonlinear Science (Vol.29, Issue 9) and may be found at https://doi.org/10.1063/1.5097826. *