# High-Dimensional Kuramoto Models on Stiefel Manifolds

Synchronize Complex Networks Almost Globally\thanksreffootnoteinfo

###### Abstract

The Kuramoto model of a system of coupled phase oscillators and its many variations are frequently used to describe synchronization phenomena in nature. This paper proposes a generalization of the Kuramoto model where each oscillator state lives on the embedding of the compact Stiefel manifold in Euclidean space . Previous work on high-dimensional Kuramoto models have largely been influenced by results and techniques that pertain to the original model. This paper uses techniques from the fields of optimization and control theory to prove a fundamental convergence result: the generalized Kuramoto model on converges to a completely synchronized state for any connected interaction topology and from almost all initial conditions provided the pair satisfies and all oscillator frequencies are equal. This result is unexpected since it could not have been predicted based on knowledge of the Kuramoto model in complex networks on the circle with homogeneous oscillator frequencies. In that case, almost global synchronization is graph dependent; it applies if the interaction topology is given by an acyclic or a sufficiently dense graph. In fact, the problem of characterizing all -synchronizing graphs is still open. This paper hence identifies a key distinguishing property of many high-dimensional generalizations of the Kuramoto model which should have important implications for both modeling and analysis of synchronization phenomena in physics and multi-agent consensus in engineering applications.

XXX]Johan Markdahl, XXX]Johan Thunberg, XXX]Jorge Goncalves

Luxembourg Centre for Systems Biomedicine, University of Luxembourg

Key words: Synchronization; Consensus; Kuramoto model; Complex networks; Stiefel manifold; Multi-agent system; Distributed control.

^{1}

^{1}footnotetext: Corresponding author Johan Markdahl. The work of Johan Markdahl is supported by the University of Luxembourg internal research project pppd.

## 1 Introduction

The Kuramoto model and its many variations are canonical models of systems of coupled phase oscillators (Hoppensteadt and Izhikevich, 2012). As such, they are abstract models that capture the essential properties observed in a wide range of synchronization phenomena (Dörfler and Bullo, 2014). However, many properties of a particular system are lost through the use of such models (Hoppensteadt and Izhikevich, 2012). In this paper we explore the properties of a system of coupled oscillators on the Stiefel manifold that includes the Kuramoto model as a special case. For systems of coupled oscillators that are subject to various constraints, a high-dimensional Stiefel manifold may provide a more faithful approximation of reality than a phase oscillator model. The orientation of a bird in a flock or a fish is a school can e.g., be modeled as an element of the circle, the sphere, or the rotation group—all of which are Stiefel manifolds. For this to be of interest, the high-dimensional system must retain some property of the original system which is lost to phase oscillator models. That is indeed the case; we prove that if the complex network of interactions is connected and a condition on the manifold topology is satisfied, then the system converges to the set of completely synchronized states from almost all initial conditions. The same cannot be said about the Kuramoto model in complex networks on the circle in the case of oscillators with homogeneous frequencies (Rodrigues et al., 2016). Under that model, guaranteed almost global synchronization requires that the complex network can be represented by a graph that is acyclic or sufficiently dense (Dörfler and Bullo, 2014). In fact, to characterize all such graphs is an open problem (Mehta et al., 2015; Canale and Monzón, 2015).

The emergence of global phenomena through local interactions in complex networks is studied in the fields of physics and control engineering from somewhat different perspectives (Tang et al., 2014; Motter, 2015). A key difference lies in the models studied and the inclusion/exclusion of agent self-dynamics (Tang et al., 2014). In physics, the focus is on canonical models of natural phenomena where the drive to synchrony competes with self-dynamics. In control applications it is usually feasible to disregard self-dynamics since engineering systems can often be designed to specification. For this reason it is also permissible to use ad-hoc control design, which may result in rather artificial dynamics being required to meet desired performance criteria. A recent trend is to combine ideas from different disciplines to study various aspects of complex networks (Dörfler and Bullo, 2014; Tang et al., 2014; Motter, 2015; Rodrigues et al., 2016). In particular, it is of interest to investigate the convergence properties of generalized Kuramoto models without self-dynamics, or equivalently, systems of coupled oscillators with homogeneous frequencies.

In this paper we show that a canonical model of coupled oscillators on (without self-dynamics) displays a strong reliability property that usually requires ad-hoc control design (Sarlette and Sepulchre, 2009; Sepulchre, 2011): almost global convergence for all connected graph topologies. From a control theory perspective this is highly desirable for two main reasons: (i) scalability; it implies that the probability of reaching consensus, i.e., converge to complete synchronization, is independently of the number of agents and (ii) reliability; since the convergence is achieved for any connected network topology the failure of a single agent is unlikely to bring about the failure of the whole system. Note that if convergence is only guaranteed from a subset of the manifold, i.e., property (i) is not satisfied as e.g., in Zhu (2013); Lageman and Sun (2016); Thunberg et al. (2018a), then the probability of guaranteed convergence to the consensus manifold decreases exponentially with .

Our convergence result hinges on the manifold topology; it holds provided the pair satisfies . This inequality is e.g., satisfied by the sphere , but not by the circle . The possible implications of this difference is intriguing. In nature, the main function of flocks (which often tend to synchrony) is to protect prey animals from predators (Caraco et al., 1980). On a speculative note, it would be interesting to investigate if aquatic and airborne animals, whose orientation in space is largely described by a point on , may have evolved to exhibit flocking to a greater extent than non-flying terrestrial animals, whose orientation is largely described by a point on , due to them having an inherit topological advantage (relatively speaking) against selective pressure. This question is provoked by the results of our paper; answering it is however outside of the scope of the present work.

### 1.1 The Stiefel manifold

The compact, real Stiefel manifold is the set of -frames in -dimensional Euclidean space (Absil et al., 2009). It can be embedded in as an analytic matrix manifold given by . The Stiefel manifold encompasses a number of important manifolds including the -sphere , the special orthogonal group , and the orthogonal group . In particular, since for all , the Stiefel manifold is of interest in modelling multi-agent systems where the norms the agents’ states are constant. Intuitively speaking, this also implies that can be conceptualized as a nowhere dense subset of a high-dimensional sphere with radius . Useful resources for working with are provided in James (1976); Edelman et al. (1998); Absil et al. (2009); Chikuse (2012). We note that in the spirit of optimization by means of dynamical systems (Helmke and Moore, 2012), our model provides a novel (distributed) averaging algorithm on with guaranteed convergence, which is of interest in signal processing applications (Kaneko et al., 2013).

### 1.2 Literature review

The main challenge for control design on nonlinear manifolds is to achieve the desired performance on a global level (Sarlette, 2009). To that aim, and a theory of almost global consensus, i.e., complete synchronization from almost all initial conditions, in complex networks on homogeneous manifolds was established (Sarlette, 2009; Sarlette and Sepulchre, 2009). Particular attention was given to the circle (Scardovi et al., 2007) and the special orthogonal group (Tron et al., 2012). In the survey Sepulchre (2011), a distinction is made between three different approaches. The first, potential shaping, uses gradient descent flows of special potential functions (Sarlette, 2009; Tron et al., 2012). Its application is restricted to connected undirected graphs since the Hessian of the potential function, a symmetric matrix, encodes the graph topology. The second, gossip algorithms, is probabilistic in nature and designed for use in discrete-time systems. It is mainly intended to model opinion dynamics, although application to networks of quantum systems is possible (Mazzarella et al., 2015). The third, dynamic consensus, introduces auxilary variables that are communicated between agents to overcome the non-convexity of the manifold (Scardovi et al., 2007; Sarlette and Sepulchre, 2009; Thunberg et al., 2018b). Gossip algorithms and dynamic consensus can accommodate quasi-strongly connected digraphs but their application is largely limited to opinion dynamics and engineering systems respectively. To these three approaches we add a forth category: canonical consensus protocols. These are obtained by projecting linear consensus protocols on the tangent space of the manifold. For certain manifolds, including the -sphere for (Markdahl et al., 2018b), almost global convergence is obtained for all connected graphs.

This paper proves that the Stiefel manifold also belongs to the category of manifolds on which canonical consensus protocols over complex network converge almost globally, provided the pair satisfies . While the protocols as are already known, their strong convergence properties are not. The reason for this was the perhaps undue attention given to the pathological cases of and . Because the Stiefel manifold encompass both nice cases, i.e., -sphere for , and pathological cases, i.e., for all , it was unclear how should be categorized. This paper establishes that the forth category encompasses the majority of instances of .

Since the Stiefel manifold includes the -sphere, special orthogonal group, and orthogonal group as special cases, there is a considerable literature of synchronization on the Stiefel manifold. Previous works that explicitly address synchronization on the Stiefel manifold are however, to the best of our knowledge, limited to Sarlette and Sepulchre (2009); Thunberg et al. (2018b), which both rely on the dynamic consensus approach. Some preliminary results of this paper appears in Markdahl et al. (2018a). As such, the model we study has not appeared in the literature previously. Several high-dimensional Kuramoto models have been proposed in the literature (Olfati-Saber, 2006; Lohe, 2009, 2010; Zhu, 2013; Ha et al., 2018). Most previous work focus on agents that evolve on the -sphere. That system is sometimes called the Lohe model (Chi et al., 2014; Zhang et al., 2018; Ha et al., 2018) in reference to Lohe (2009, 2010). Numerous related high-dimensional oscillator models exist. The framework of synchronization on a nonlinear space has even been formulated in terms of the intrinsic geometry on (rather) general Riemannian manifolds (Tron et al., 2013).

There are a number of works that discuss collective motion on the -sphere, including Ritort (1998); Olfati-Saber (2006); D. A. Paley (2009); Lohe (2009, 2010); Zhu (2013); Li and Spong (2014); Chi et al. (2014); Tanaka (2014); Caponigro et al. (2015); Li (2015); Markdahl (2015); Lageman and Sun (2016); Aydogdu et al. (2017); Song et al. (2017); Al-Abri and Zhang (2017); Markdahl et al. (2018b); Lohe (2018); Zhang et al. (2018); Thunberg et al. (2018a); Chandra et al. (2018); Ha et al. (2018). The combination of graph topology and stability results (if present) are either digraphs and convergence from initial conditions restricted to a hemisphere (Zhu, 2013; Caponigro et al., 2015; Lageman and Sun, 2016; Thunberg et al., 2018a; Zhang et al., 2018) or a specific class of graphs (acyclic (Markdahl, 2015; Zhang et al., 2018), essentially a single cycle, (Markdahl, 2015; Song et al., 2017), or complete (Olfati-Saber, 2006; Zhu, 2013; Li and Spong, 2014; Lohe, 2018; Zhang et al., 2018; Chandra et al., 2018; Ha et al., 2018)) and almost global consensus. Actually, almost global consensus holds for any connected graph provided (Markdahl et al., 2018b). Some papers address the case of self-dynamics (Chi et al., 2014; Chandra et al., 2018; Ha et al., 2018). Applications for consensus on include autonomous reduced attitude synchronization and balancing (Song et al., 2017), synchronization of interacting tops (Ritort, 1998), modeling of collective motion in flocks (Couzin et al., 2002; Al-Abri and Zhang, 2017), synchronization in planetary scale sensor networks (D. A. Paley, 2009), and the consensus problem in opinion dynamics (Caponigro et al., 2015; Aydogdu et al., 2017). Applications on include synchronization of quantum bits (Lohe, 2010) and full attitude synchronization by means of the quaternion representation of a rotation matrix, although the resulting closed-loop system is unwinding on (S. P. Bhat and D. S. Bernstein, 2000).

There is a vast literature about consensus on , in particular for synchronization of the Kuramoto model. We refer the interested reader to the survey papers Arenas et al. (2008); Dörfler and Bullo (2014); Rodrigues et al. (2016). For synchronization on , there is also a large literature. The fact that some of the systems considered are direct generalizations of the Kuramoto model is rarely made explicit. This is partly due to the widespread use of various representations of which obscure such relations and the preference for rigid-body attitude dynamics based on Euler’s equations of motion over models that presume angular velocity based actuation. Examples of systems that are generalization of the Kuramoto model include Sarlette and Sepulchre (2009) for the case of a kinematic model and Sarlette et al. (2009) for the case of a dynamic model. A model of quantum synchronization over quantum networks on the unitary group also exists (Lohe, 2009, 2010; Ha et al., 2018).

## 2 Problem Formulation

### 2.1 Preliminaries

The Euclidean (or Frobenius) inner product of is . The norm of is given by . Our approach is based on the canonical embedding of the compact, real Stiefel manifold as an analytic matrix manifold in (Edelman et al., 1998; Absil et al., 2009), . The most important instances of Stiefel manifolds are: the -sphere , the special orthogonal group , and the orthogonal group . The pair forms a smooth Riemannian manifold. Note that is not the metric tensor obtained from the canonical embedding of as a quotient manifold in (Edelman et al., 1998). The canonical embedding in and its corresponding metric, the so-called chordal distance , is often preferred from a computational point of view (Conway et al., 1996).

Define the projections and . The tangent space of at is given by

The projection onto the tangent space, , is given by

(1) |

From a computational perspective, it is sometimes preferable to use the equivalent expression

(2) |

since it consists of three terms compared to the four in (1).

The intrinsic gradient on (in terms of the Euclidean inner product ) of a function is given by , where is any smooth extension of on , and denotes the extrinsic gradient in Euclidean space. The state of the system we consider is an -tuple of Stiefel matrices. The gradient of a scalar with respect to the state is thus also an -tuple of matrices.

A graph is a pair where and is a set of subsets of , where each subset is of cardinality two. Throughout this paper, if an expression depends on an edge and two nodes , then it is implicitly understood that . Moreover, any graph under consideration is assumed to be connected, i.e., it contains a tree subgraph with . Each element also corresponds to an agent. Items associated with a specific agent carry the subindex ; we let denote the state of an agent, the projection onto the tagent space , the neighbor set of , the gradient of with respect to , etc.

### 2.2 Gradient descent flow on the Stiefel manifold

The synchronization set, or consensus manifold, of the -fold product of a Stiefel manifold is defined as

The synchronization set is a (sub)manifold; it is diffeomorphic to by the map . Let be the chordal distance between agent and . Given a graph , define the potential function by

(3) |

where satisfies for all . Note that is a real-analytic function and that . Let . Let be a (any) smooth extension of obtained by relaxing the requirement to for all . The system we study is the gradient descent flow on given by

(4) | ||||

(5) | ||||

where for all .

### 2.3 Problem statement

The aim of this paper is classify each instance of as either satisfying or not satisfying the following property: the gradient descent flow (4), (5) with interaction topology given by any connected graph converges to the consensus manifold from almost all initial conditions. We know that this property is satisfied by the -sphere for (Markdahl et al., 2018b). We also know that it is not satisfied on the circle (Sarlette, 2009) nor on (Tron et al., 2012; Markdahl et al., 2018b). Note that and . The main result of our paper states the high-dimensional Kuramoto model synchronizes complex network almost globally provided the pair satisfies .

### 2.4 High-dimensional Kuramoto model

The high-dimensional Kuramoto model in complex networks over the Stiefel manifold is given by

(6) |

where , . The variables and are generalizations of the frequency term that appear in the self-dynamics of the Kuramoto model.

###### Proposition 1

PROOF. To get (7), let and set , for all . Note that is given by . The restriction of implies that for all .

To get (8), let and set for all . Note that is given by .

### 2.5 Distributed control design on the Stiefel manifold

A distributed system is a network of subsystems, each of which is conceptualized as an autonomous agent (Mesbahi and Egerstedt, 2010). Distributed control protocols serve to automate decision making at each node of the network so that some overall goal can be achieved on a global level. For reasons such as limited sensing and communication capabilities of the agents, scalability in terms of the number of agents , information security, or energy efficiency, it is desirable to design control protocols that require a minimum of both information access on the behalf of each agent and communication channels between agents. Essentially, if the network is modeled by a sparse, possible acyclic, connected graph, then the system should still be able to achieve its goal.

On the Stiefel manifold, we assume the motion of each agent to be generated on a kinematic level. The dynamics of each agent hence take the form of a single integrator, where is the input signal. To reach consensus, an agent acts based on its local view of the current degree of synchronization in the system. Agent does not know , but we assume it can calculate

at the current time. Symmetry of gives whereby it follows that . Since agent can evaluate at its current position, it is reasonable to assume that it can also calculate . The agent can hence realize its own role, given by (5), in the overall dynamics given by (4).

### 2.6 Local stability

As we stated in the introduction, the main challenge for feedback control on compact manifolds is to ensure good performance on a global level (Sarlette, 2009). In this paper we focus on almost global asymptotical stability. The local stability properties of the system given by (4) are summarized in Proposition 2. This result states some rather generic properties of gradient descent flows. We do not give a proof, but refer the interested reader to Absil and Kurdyka (2006); Lageman (2007); Helmke and Moore (2012).

###### Proposition 2

The gradient descent flow (4) converges to a critical point of . The manifold is separated by a nonzero constant distance from all other critical points. All global minimizers of belong to . The sublevel sets

are forward invariant.

Let be the set of critical points of . Let belong to the same connected component of as does. From Proposition 2, it follows that the region of attraction of contains all sets such that . Moreover, the synchronization manifold is stable.

### 2.7 -Synchronizing graphs

###### Definition 3

A graph is -synchronizing if all minimizers of belong to .

This definition is a direct extension of the concept of -synchronizing graphs with regard to the Kuramoto model in complex networks (Canale and Monzón, 2007).

###### Definition 4

An equilibrium set of system (4) is referred to as almost globally asymptotically stable (AGAS) if it is stable and attractive from all initial conditions , where has Haar measure zero on .

It is not possible to globally stabilize an equilibrium set on a compact manifold by means of continuous, time-invariant feedback (S. P. Bhat and D. S. Bernstein, 2000). This obstruction, which is due to topological reasons, does not exclude the possibility of AGAS. The property of a graph being -synchronizing implies that is AGAS:

###### Proposition 5

The consensus manifold is an AGAS equilibrium set of the dynamics (4) when satisfies if the graph is -synchronizing.

Because of this implication, we may establish that the consensus manifold is AGAS by solving an optimization problem, i.e., by showing that

Here, it is understood that the operator yields all local minimum points, even if they are not global.

## 3 Main Result

###### Theorem 6

Let the pair satisfy and suppose that is connected. The consensus manifold is almost globally asymptotically stable. All minimizers of the potential function given by (3) belong to .

PROOF. The main ideas of the proof of the first statement are rather basic, although the calculations involved are extensive. We begin with a brief sketch. The details of every step is provided in Section 4.1 to 4.6.

Let denote the quadratic form obtained from the intrinsic Hessian of evaluated at a critical point on . The second-order necessary conditions for constrained optimization on imply that is weakly positive at any minimum of . Our goal is to exclude minimality of all equilibria by finding such that . All elements of are global minimizers since for all arguments and .

We consider perturbations towards , i.e., for all and some . We do not need to determine an exact expression for the desired perturbation to show that , it suffices to prove that it exists. Showing that can assume negative values is done by solving a nonconvex constrained optimization problem to minimize an upper bound of over . This is done using the necessary second-order optimality conditions. For any critical point and pair such that , we can show that assumes negative values. Such a cannot be a minimum point of . Throughout these steps, we never utilize any particular property of the graph topology. Connectedness is however required for complete synchronization. It follows that all connected graphs are -synchronizing. By Proposition 5, this implies that is AGAS.

## 4 Detailed Proof of Main Result

This section contains the proof of our main result, Theorem 6. A sketch of the proof is provided in Section 3.

### 4.1 Equilibria are critical points

We start by characterizing the equilibria of system (4), (5). At an equilibrium,

Since the two terms in this expression are orthogonal, we get

(10) |

Assume (10) holds. Define . Since , it follows that . Hence for some . Moreover, since , we find that is symmetric. The matrix is closely related to the Lagrange multipliers for the constraints for all . Note that equilibria and critical points coincide for gradient descent flows (Helmke and Moore, 2012), so we do not need to formulate the Lagrangian.

### 4.2 The intrinsic Hessian

The first step in the proof sketch of Theorem 6 is to determine the intrinsic Hessian . Let be a (any) smooth extension of obtained by relaxing the constraint to for all . Take a and calculate

Using the rules governing derivatives of inner products with respect to matrices, introducing , after a few calculations, we obtain

Evaluate at an equilibrium, where and is symmetric by Section 4.1, to find

The intrinsic Hessian is a -tensor consisting of blocks , which are obtained by projecting the extrinsic Hesssian on the tangent space of

### 4.3 The quadratic form

Consider the quadratic form obtained from the intrinsic Hessian evaluated at an equilibrium for a perturbation