Robust Angular Synchronization via Directed Graph Neural Networks
The angular synchronization problem aims to accurately estimate (up to a
constant additive phase) a set of unknown angles $\theta_1, \dots,
\theta_n\in[0, 2\pi)$ from $m$ noisy measurements of their offsets
$\theta_i-\theta_j \;\mbox{mod} \; 2\pi.$ Applications include, for example,
sensor network localization, phase retrieval, and distributed clock
synchronization. An extension of the problem to the heterogeneous setting
(dubbed $k$-synchronization) is to estimate $k$ groups of angles
simultaneously, given noisy observations (with unknown group assignment) from
each group. Existing methods for angular synchronization usually perform poorly
in high-noise regimes, which are common in applications. In this paper, we
leverage neural networks for the angular synchronization problem, and its
heterogeneous extension, by proposing GNNSync, a theoretically-grounded
end-to-end trainable framework using directed graph neural networks. In
addition, new loss functions are devised to encode synchronization objectives.
Experimental results on extensive data sets demonstrate that GNNSync attains
competitive, and often superior, performance against a comprehensive set of
baselines for the angular synchronization problem and its extension, validating
the robustness of GNNSync even at high noise levels.
Yixuan He, Gesine Reinert, David Wipf, Mihai Cucuringu
# Robust Angular Synchronization via Directed Graph Neural Networks

###### Abstract

The angular synchronization problem aims to accurately estimate (up to a constant additive phase) a set of unknown angles \(\theta_{1},\ldots,\theta_{n}\in[0,2\pi)\) from \(m\) noisy measurements of their offsets \(\theta_{i}-\theta_{j}\) mod \(2\pi\). Applications include, for example, sensor network localization, phase retrieval, and distributed clock synchronization. An extension of the problem to the heterogeneous setting (dubbed \(k\)-synchronization) is to estimate \(k\) groups of angles simultaneously, given noisy observations (with unknown group assignment) from each group. Existing methods for angular synchronization usually perform poorly in high-noise regimes, which are common in applications. In this paper, we leverage neural networks for the angular synchronization problem, and its heterogeneous extension, by proposing GNNSync, a theoretically-grounded end-to-end trainable framework using directed graph neural networks. In addition, new loss functions are devised to encode synchronization objectives. Experimental results on extensive data sets demonstrate that GNNSync attains competitive, and often superior, performance against a comprehensive set of baselines for the angular synchronization problem and its extension, validating the robustness of GNNSync even at high noise levels.

## 1 Introduction

The group synchronization problem has received considerable attention in recent years, as a key building block of many computational problems. Group synchronization aims to estimate a collection of group elements, given a small subset of potentially noisy measurements of their pairwise ratios \(\Upsilon_{i,j}=g_{i}\,g_{j}^{-1}\). Typical applications (for different underlying groups) include rotation-averaging in 3D computer vision (Arrigoni & Fusiello, 2020; Janco & Bendory, 2022) (over the group SO(3) of 3D rotations), the molecule problem in structural biology (Cucuringu et al., 2012b) (over SO(3)), solving jigsaw puzzles (Huroyan et al., 2020) (over the group \(\mathbb{Z}_{4}\) of the integers \(\{0,1,2,3\}\) with addition mod 4 as the group operation), recovering a global ranking from pairwise comparisons (He et al., 2022a; Cucuringu, 2016) (over the group \(\mathbb{Z}_{n}\), resp., SO(2)), and sensor network localization (Cucuringu et al., 2012a) (over the Euclidean group of rigid motions \(\text{Euc}(2)=\mathbb{Z}_{2}\times\text{SO(2)}\times\mathbb{R}^{2}\)).

An important special case is _angular synchronization_, also referred to as _phase synchronization_, which can be viewed as group synchronization over SO(2). The angular synchronization problem aims at obtaining an accurate estimation (up to a constant additive phase) for a set of unknown angles \(\theta_{1},\ldots,\theta_{n}\in[0,2\pi)\) from \(m\) noisy measurements of their pairwise offsets

Figure 1: Sensor network localization map.
\(\theta_{i}-\theta_{j}\mod\ 2\pi\). This problem has a wide range of applications, such as distributed clock synchronization over wireless networks (Giridhar & Kumar, 2006), image reconstruction from pairwise intensity differences (Yu, 2009; 2011), phase retrieval (Forstner et al., 2020; Iwen et al., 2020), and sensor network localization (SNL) (Cucuringu et al., 2012a). In engineering, the SNL problem seeks to reconstruct the 2D coordinates of a cloud of points from a sparse set of pairwise noisy Euclidean distances; in typical divide-and-conquer approaches that aid with scalability, one first computes a local embedding of nearby points (denoted as _patches_) and is left with the task of stitching the patches together in a globally consistent embedding (Cucuringu et al., 2012a). Fig. 1 is an example of SNL on the U.S. map, where our method recovers city locations (in blue) and aims to match ground-truth locations (in red). Most works in the SNL literature that focus on the methodology development consider only purely synthetic data sets in their experiments; here we consider a real-world data set (actual 2D layout with different levels of densities of cities across the U.S. map), and add synthetic noise to perturb the local patch embeddings for testing the robustness to noise of the angular synchronization component.

An extension of angular synchronization to the heterogeneous setting is \(k\)-_synchronization_, introduced in Cucuringu & Tyagi (2022), and motivated by real-world graph realization problems (GRP) and ranking. GRP aims to recover coordinates of a cloud of points in \(\mathbb{R}^{d}\), from a sparse subset (edges of a graph) of noisy pairwise Euclidean distances (the case \(d=2\) is the above SNL problem). The motivation for \(k\)-_synchronization_ arises in structural biology, where the distance measurements between pairs of atoms may correspond to \(k\) different configurations of the molecule, in the case of molecules with multiple conformations. In ranking applications, the \(k=2\) sets of disjoint pairwise measurements may correspond to two different judges, whose latent rankings we aim to recover.

A key limitation of existing methods for angular synchronization is their poor performance in the presence of considerable noise. High noise levels are not unusual; measurements in SO(3) can have large outliers in certain biological settings (cryo-EM and NMR spectroscopy), see for example Cucuringu et al. (2012b). Therefore, we need new methods to push the boundary of signal recovery when there is a high level of noise. While neural networks (NNs), in principle, could be trained to address high noise regimes, the angular synchronization problem is not directly amenable to a standard NN architecture due to the directed graph (digraph) structure of the underlying data measurement process and the underlying group structure; hence the need for a customized graph neural network (GNN) architecture and loss function for this task. Here we propose a GNN method called GNNsync for angular synchronization, with a novel cycle loss function, which downweights noisy observations, and explicitly enforces cycle consistency as a quality measure. GNNsync's novelty does not lie in simply applying a data-driven NN to this task, but rather in proposing a framework for handling the pairwise comparisons encoded in a digraph, accounting for the underlying SO(2) group structure, and designing a loss function for increased robustness to noise and outliers, with theoretical support.

Our main contributions are summarized as follows.

\(\bullet\) We demonstrate how the angular synchronization problem can be recast as a theoretically-grounded directed graph learning task by first incorporating the inductive biases of classical estimators within the design of a more robust GNN architecture, called GNNsync, and then pairing with a novel training loss that exploits cycle consistency to help disambiguate unknown angles.

\(\bullet\) We perform extensive experiments comparing GNNsync with existing state-of-the-art algorithms from the angular synchronization and \(k\)-synchronization literature, across a variety of synthetic outlier models at various density and noise levels, and on a real-world application. GNNsync attains leading performance, especially in high noise regimes, validating its robustness to noise.

## 2 Related work

### Angular synchronization

The seminal work of Singer (2011) introduced spectral and semidefinite programming (SDP) relaxations for angular synchronization. For the spectral relaxation, the estimated angles are given by the eigenvector corresponding to the largest eigenvalue of a Hermitian matrix \(\mathbf{H}\), whose entries are given by \(\mathbf{H}_{i,j}=\exp(\iota\mathbf{A}_{i,j})\mathbbm{1}\left(\mathbf{A}_{i,j} \neq 0\right)\), where \(\iota\) is the imaginary unit, and \(\mathbf{A}_{i,j}\) is the observed potentially noisy offset \(\theta_{i}-\theta_{j}\) mod \(2\pi\). Singer (2011) also provided an SDP relaxation involving the same matrix \(\mathbf{H}\), and empirically demonstrated that the spectral and SDP relaxations yield similar experimental results. A row normalization was introduced to \(\mathbf{H}\) prior to the eigenvector computation by Cucuringu et al. (2012a), which showed improved results. Cucuringu et al. (2012b) generalized this approach to
the 3D setting \(\text{Euc}(3)=\mathbb{Z}_{2}\times\text{SO}(3)\times\mathbb{R}^{3}\), and incorporated into the optimization pipeline the ability to operate in a semi-supervised setting, where certain group elements are known a-priori. Cucuringu and Tyagi (2022) extended the angular synchronization problem to a heterogeneous setting, to the so-called \(k\)_-synchronization_ problem, whose goal is to estimate \(k\) sets of angles simultaneously, given only the graph union of noisy pairwise offsets, which we also explore in our experiments. The key idea in their work is to estimate the \(k\) sets of angles from the top \(k\) eigenvectors of the angular embedding matrix \(\mathbf{H}\).

Boumal (2016) modeled the angular (phase) synchronization problem as a least-squares non-convex optimization problem, and proposed a modified version of the power method called the Generalized Power Method (GPM), which is straightforward to implement and free of parameter tuning. GPM often attains leading performance among baselines in our experiments, and the iterative steps in the GPM method motivated the design of the projected gradient steps in our GNNSync architecture. However, GPM is not directly applicable to \(k\)-synchronization with \(k>1\) while GNNSync is. For \(k=1,\) GNNSync tends to perform significantly better than GPM at high noise levels. Bandeira et al. (2017) studied the tightness of the maximum likelihood semidefinite relaxation for angular synchronization, where the maximum likelihood estimate is the solution to a nonbipartite Grothendieck problem over the complex numbers. A truncated least-squares approach was proposed by Huang et al. (2017) that minimizes the discrepancy between the estimated angle differences and the observed differences under some constraints. Gao and Zhao (2019) tackled the angular synchronization problem with a multi-frequency approach. Liu et al. (2023) unified various group synchronization problems over subgroups of the orthogonal group. Filbir et al. (2021) provided recovery guarantees for eigenvector relaxation and semidefinite convex relaxation methods for weighted angular synchronization. Lerman and Shi (2022) applied a message-passing procedure based on cycle consistency information, to estimate the corruption levels of group ratios and consequently solve the synchronization problem, but the method is focused on the restrictive setting of adversarial or uniform corruption and sufficiently small noise. In addition, Lerman and Shi (2022) requires post-processing based on the estimated corruption levels to obtain the group elements, while GNNSync is trained end-to-end. Maunu and Lerman (2023) utilized energy minimization ideas, with a variant converging linearly to the ground truth rotations.

### Directed graph neural networks

Digraph node embeddings can be effectively learned via directed graph neural networks (He et al., 2022c). For learning such an embedding, Tong et al. (2020) constructed a GNN using higher-order proximity. Zhang et al. (2021) built a complex Hermitian Laplacian matrix and proposed a spectral digraph GNN. He et al. (2022b) introduced imbalance objectives for digraph clustering. Our GNNSync framework can readily incorporate any existing digraph neural network.

### Relationship with other group synchronization methods

Angular synchronization outputs can be used to obtain global rankings by using a one-dimensional ordering. To this end, recovering rankings of \(n\) objects from pairwise comparisons can be viewed as group synchronization over \(\mathbb{Z}_{n}\). To recover global rankings from pairwise comparisons, GNNRank (He et al., 2022a) adopted an unfolding idea to add an inductive bias from Fogel et al. (2014) to the NN architecture. Inspired by He et al. (2022a), we adapt their framework to borrow strength from solving a related problem. We adapt their "innerproduct" variant to \(k\)-synchronization, remove the 1D ordering at the end of the GNNRank framework, and rescale the estimated quantities to the range \([0,2\pi)\). We also borrow strength from the projected gradient steps in GPM (Boumal, 2016) and add projected gradient steps to our GNNSync architecture. Another key novelty is that we devise novel objectives, which reflect the angular structure of the data, to serve as our training loss functions. The architectures are also very different: While in GNNRank the proximal gradient steps play a vital role from an unrolling perspective, and the whole architecture could be viewed as an unrolling of the SerialRank algorithm, here, although we borrow strength from the GPM method, the whole architecture is different from merely unrolling GPM. Furthermore, the baselines serve as initial guesses for the "proximal baseline" variant in GNNRank, but serve as input node features in our approach.

Other methods have been introduced for group synchronization, but mostly in the context of SO(3). Shi and Lerman (2020) proposed an efficient algorithm for synchronization over SO(3) under high levels of corruption and noise. Shi et al. (2022) provided a novel quadratic programming formulation for estimating the corruption levels, but again its focus is on SO(3). Unrolled algorithms (which are NNs) were introduced for SO(3) in Janco and Bendory (2022). While an adaptation to SO(2) may be possible in principle, as its objective functions are based on the level of agreement between the
estimated angles and ground-truth, its experiments require ground-truth during training, usually not available in practice. In contrast, our GNNsync framework can be trained without any known angles.

## 3 Problem definition

The _angular synchronization_ problem aims at obtaining an accurate estimation (up to a constant additive phase) for a set of \(n\) unknown angles \(\theta_{1},\ldots,\theta_{n}\in[0,2\pi)\) from \(m\) noisy measurements of their offsets \(\theta_{i}-\theta_{j}\) mod \(2\pi\), for \(i,j\in\{1,\ldots,n\}\). We encode the noisy measurements in a digraph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\), where each of the \(n\) elements of the node set \(\mathcal{V}\) has as attribute an angle \(\theta_{i}\in[0,2\pi)\). The edge set \(\mathcal{E}\) represents pairwise measurements of the angular offsets \((\theta_{i}-\theta_{j})\) mod \(2\pi\). The weighted directed graph has a corresponding adjacency matrix \(\mathbf{A}\) with \(\mathbf{A}_{i,j}=(\theta_{i}-\theta_{j})\) mod \(2\pi\geq 0\) where at most one of \(\mathbf{A}_{i,j}\) and \(\mathbf{A}_{j,i}\) can be nonzero by construction; this representation is used to reduce computational complexity. In general, the matrix \(\mathbf{A}\) is not symmetric. Estimating the unknown angles from noisy offsets amounts to assigning an estimate \(r_{i}\in[0,2\pi)\) to each node \(i\in\mathcal{V}\).

An extension of the above problem to the heterogeneous setting is the \(k\)_-synchronization_ problem, which is defined as follows. We are given only the graph union of \(k\) digraphs \(\mathcal{G}_{1},\ldots,\mathcal{G}_{k}\), with the same node set and disjoint edge sets, which encode noisy measurements of \(k\) sets \((\theta_{i,l}-\theta_{j,l})\) mod \(2\pi\), for \(l\in\{1,\ldots,k\},i,j\in\{1,\ldots,n\}\), of angle differences modulo \(2\pi\). Its adjacency matrix is denoted by \(\mathbf{A}\). The problem is to estimate these \(k\) sets of \(n\) unknown angles \(\theta_{i,l}\in[0,2\pi)\), \(\forall l\in\{1,\ldots,k\},i\in\{1,\ldots,n\},\) simultaneously. Note that we are given only \(\mathcal{G}=\mathcal{G}_{1}\cup\cdots\cup\mathcal{G}_{k}\) and the value of \(k\), and each edge in \(\mathcal{G}\) belongs to exactly one of \(\mathcal{G}_{1},\ldots,\mathcal{G}_{k}\). To unify notations, we view the normal angular synchronization problem as a special case of the more general \(k\)-synchronization problem where \(k=1\).

## 4 Loss and evaluation

### Loss and evaluation for angular synchronization

For a vector \(\mathbf{r}=[r_{1},\ldots,r_{n}]^{\top}\) with estimated angles as entries, we define \(\mathbf{T}=[(\mathbf{r}\mathbf{1}^{\top}-\mathbf{1}\mathbf{r}^{\top})\text{ mod }2\pi]\in\mathbb{R}^{n\times n}\). Then \(\mathbf{T}_{i,j}=(r_{i}-r_{j})\) mod \(2\pi\) estimates \(\mathbf{A}_{i,j}\). We only compare \(\mathbf{T}\) with \(\mathbf{A}\) at locations where \(\mathbf{A}\) has nonzero entries. We introduce the residual matrix \(\mathbf{M}\) with entries

\[\mathbf{M}_{i,j}=\min\left((\mathbf{T}_{i,j}-\mathbf{A}_{i,j})\text{ mod }2\pi,(\mathbf{A}_{i,j}-\mathbf{T}_{i,j})\text{ mod }2\pi\right)\]

if \(\mathbf{A}_{i,j}\neq 0\), and \(\mathbf{M}_{i,j}=0\) if \(\mathbf{A}_{i,j}=0\). Then our upset loss is defined as

\[\mathcal{L}_{\text{upset}}=\left\|\mathbf{M}\right\|_{F}/t, \tag{1}\]

where the subscript \(F\) means Frobenius norm, and \(t\) is the number of nonzero elements in \(\mathbf{A}\). Despite the non-differentiability of the loss function, using the concept of a limiting subdifferential from Li et al. (2020) we can give the following theoretical guarantee on the minimization of eq. (1); its proof is in Appendix (App.) A.1, where also the case of general \(k\) is discussed.

**Proposition 1**.: _Every local minimum of eq. (1) is a directional stationary point of eq. (1)._

For _evaluation_, we employ a Mean Square Error (MSE) function with angle corrections, considered in Singer & Shkolnisky (2011). As the offset measurements are unchanged if we shift all angles by a constant, denoting the ground-truth angle vector as \(\mathbf{R}\), this evaluation function can be written as

\[\mathcal{D}_{\text{MSE}}(\mathbf{r},\mathbf{R})=\min_{\theta_{0}\in[0,2\pi)} \sum_{i=1}^{n}[\min(\delta_{i}\text{ mod }2\pi,(-\delta_{i})\text{ mod }2\pi)]^{2}, \tag{2}\]

where \(\delta_{i}=r_{i}+\theta_{0}-\theta_{i},\ \forall i=1,\ldots,n\). Additional implementation details are provided in App. C.4.

### Cycle consistency relation

For noiseless observations, every cycle in the angular synchronization problem (\(k=1\)) or every cycle whose edges correspond to the same offset graph \(\mathcal{G}_{l}\) (\(k>1\)) satisfy the _cycle consistency_ relation that the angle sum mod \(2\pi\) is \(0\). For 3-cycles \((i,j,q)\), such that \(\mathbf{A}_{i,j}\cdot\mathbf{A}_{j,q}\cdot\mathbf{A}_{q,i}>0\), this leads to

\[(\mathbf{A}_{i,j}+\mathbf{A}_{j,q}+\mathbf{A}_{q,i})\text{ mod }2\pi=(\theta_{i}- \theta_{j}+\theta_{j}-\theta_{q}+\theta_{q}-\theta_{i})\text{ mod }2\pi=0,\]

as \((a+b\text{ mod }m)=\{(a\text{ mod }m)+(b\text{ mod }m)\text{ mod }m\}\). Hence we obtain the 3-cycle condition

\[(\mathbf{A}_{i,j}+\mathbf{A}_{j,q}+\mathbf{A}_{q,i})\text{ mod }2\pi=0,\forall(i,j,q)\text{ \ such that }\mathbf{A}_{i,j}\cdot\mathbf{A}_{j,q}\cdot\mathbf{A}_{q,i}>0. \tag{3}\]

With \(\mathbb{T}=\{(i,j,q):\mathbf{A}_{i,j}\cdot\mathbf{A}_{j,q}\cdot\mathbf{A}_{q,i} >0\}\), we define the _cycle inconsistency level_\(\frac{1}{|\mathbb{T}|}\sum_{(i,j,q)\in\mathbb{T}}[(\mathbf{A}_{i,j}+ \mathbf{A}_{j,q}+\mathbf{A}_{q,i})\text{ mod }2\pi]\). We devise a loss function to minimize the cycle inconsistency level with reweighted edges.
### Loss and evaluation for general k-synchronization

The upset loss for general \(k\) is defined similarly as in Sec. 4.1. Recall that the observed graph \(\mathcal{G}\) has adjacency matrix \(\mathbf{A}\). Given \(k\) groups of estimated angles \(\{r_{1,l},\ldots,r_{n,l}\}\), \(l=1,\ldots,k\), we define the matrix \(\mathbf{T}^{(l)}\) with entries \(\mathbf{T}^{(l)}_{i,j}=(r_{i,l}-r_{j,l})\) mod \(2\pi\), for \(i,j\in\{1,\ldots,n\},l\in\{1,\ldots,k\}\). We define \(\mathbf{M}^{(l)}\) by \(\mathbf{M}^{(l)}_{i,j}=\min((\mathbf{T}^{(l)}_{i,j}-\mathbf{A}_{i,j})\) mod \(2\pi,(\mathbf{A}_{i,j}-\mathbf{T}^{(l)}_{i,j})\) mod \(2\pi)\) if \(\mathbf{A}_{i,j}\neq 0\), and \(\mathbf{M}^{(l)}_{i,j}=0\) if \(\mathbf{A}_{i,j}=0\). Define \(\mathbf{M}\) by \(\mathbf{M}_{i,j}=\min_{le\{1,\ldots,k\}}\mathbf{M}^{(l)}_{i,j}\). The upset loss is as in eq.1, \(\mathcal{L}_{\text{upset}}=\left\|\mathbf{M}\right\|_{F}/t\).

In addition to \(\mathcal{L}_{\text{upset}}\), we introduce another option as a loss function based on the cycle consistency relation from Sec. 4.2, which adds a regularization that helps in guiding the learning process for certain challenging scenarios (e.g., with sparser \(\mathcal{G}\) or larger \(k\)). Since measurements are typically noisy, we first estimate the corruption level by entries in \(\mathbf{M}\), and use them to construct a confidence matrix \(\mathbf{C}\) by \(\mathbf{C}_{i,j}=\frac{1}{1+\mathbf{M}_{i,j}}\mathbbm{1}\left(\mathbf{A}_{i,j }\neq 0\right)\), then normalize the entries by \(\tilde{\mathbf{C}}_{i,j}=\mathbf{C}_{i,j}\frac{\sum_{u,v}\mathbf{A}_{u,v}}{ \sum_{u,v}\mathbf{A}_{u,v}\cdot\mathbf{C}_{u,v}}\). The normalization is chosen such that \(\sum_{i,j}\mathbf{A}_{i,j}\tilde{\mathbf{C}}_{i,j}=\sum_{u,v}\mathbf{A}_{u,v}\). Keeping the sum of edge weights constant is carried out in order to avoid reducing the cycle inconsistency level by only rescaling edge weights but not their relative magnitudes. Based on the confidence matrix \(\tilde{\mathbf{C}}\), we reweigh edges in \(\mathcal{G}\) to obtain an updated input graph, whose adjacency matrix is the Hadamard product \(\mathbf{A}\odot\tilde{\mathbf{C}}\). This graph attaches larger weights to edges \(\mathbf{A}_{i,j}\) for which \(\mathbf{T}^{(l)}_{i,j}\) is a good estimate when the edge \((i,j)\) belongs to graph \(\mathcal{G}_{l}\).

As the graph assignment of an edge \((i,j)\) is not known during training, we estimate it by

\[g(i,j)=\arg\min_{le\{1,\ldots,k\}}\mathbf{M}^{(l)}_{i,j},\text{ and set }g(j,i)=g(i,j), \tag{4}\]

thus obtaining our estimated graphs \(\tilde{\mathcal{G}}_{1},\ldots,\tilde{\mathcal{G}}_{k}\), which are also edge disjoint. Next, aiming to minimize 3-cycle inconsistency of the updated input graph given our graph assignment estimates, we introduce a loss function denoted as the _cycle inconsistency loss_\(\mathcal{L}_{\text{cycle}}\); for simplicity, we only focus on 3-cycles (triangles). We interpret the matrix \(\tilde{\mathbf{A}}=(\mathbf{A}\odot\tilde{\mathbf{C}}-(\mathbf{A}\odot \tilde{\mathbf{C}})^{\top})\) mod \(2\pi\) as the adjacency matrix of another weighted directed graph \(\tilde{\mathcal{G}}\). The entry \(\tilde{A}_{i,j}\) of the new adjacency matrix approximates angular differences of a reweighted graph, with noisy observations downweighted. Note that we only reweigh the adjacency matrix in the cycle loss definition, but do not update the input graph. The underlying idea is that this updated denoised graph may display higher cycle consistency than the original graph. From our graph assignment estimates, we obtain estimated adjacency matrices \(\tilde{\mathbf{A}}^{(l)}\) for \(l\in\{1,\ldots,k\}\), where \(\tilde{\mathbf{A}}^{(l)}_{i,j}=1(g(i,j)=l)\tilde{\mathbf{A}}_{i,j}\). Let \(\mathbb{T}^{(l)}=\{(i,j,q):\tilde{A}^{(l)}_{i,j}\tilde{A}^{(l)}_{j,q}\tilde{A} ^{(l)}_{q,i}>0\}\) denote the set of all triangles in \(\tilde{\mathcal{G}}_{l}\), and set \(S^{(l)}_{i,j,q}=\tilde{\mathbf{A}}^{(l)}_{i,j}+\tilde{\mathbf{A}}^{(l)}_{j,q} +\tilde{\mathbf{A}}^{(l)}_{q,i}\), for \((i,j,q)\in\mathbb{T}^{(l)}\). We define

\[\mathcal{L}^{(l)}_{\text{cycle}}=\frac{1}{|\mathbb{T}^{(l)}|}\sum_{(i,j,q)\in \mathbb{T}^{(l)}}\min(S^{(l)}_{i,j,q}\text{ mod }2\pi,(-S^{(l)}_{i,j,q})\text{ mod }2\pi) \tag{5}\]

and set \(\mathcal{L}_{\text{cycle}}=\frac{1}{k}\sum_{l=1}^{k}\mathcal{L}^{(l)}_{\text{ cycle}}\). The default training loss for \(k\geqslant 2\) is \(\mathcal{L}_{\text{cycle}}\) or \(\mathcal{L}_{\text{upset}}\) alone; in the experiment section, we also report the performance of a variant based on \(\mathcal{L}_{\text{upset}}+\mathcal{L}_{\text{cycle}}\).

For evaluation, we compute \(\mathcal{D}_{\text{MSE}}\) with eq.2, for each of the \(k\) sets of angles, and consider the average. As the ordering of the \(k\) sets can be arbitrary, we consider all permutations of \(\{1,\ldots,k\}\), denoted by \(perm(k)\). Denoting the ground-truth angle matrix as \(\mathbf{R}\), whose \((i,l)\) entry is the ground-truth angle \(\theta_{i,l}\), and the \(l\)-th entry of the permutation \(pe\) by \(pe(l)\), the final MSE value is

\[\mathcal{D}_{\text{MSE}}(\mathbf{r},\mathbf{R})=\frac{1}{k}\min_{peperm(k)}\sum _{l=1}^{k}\mathcal{D}_{\text{MSE}}(\mathbf{r}_{\cdot,pe(l)},\mathbf{R}_{ \cdot,l}). \tag{6}\]

Note that the MSE loss is not used during training as we do not have any ground-truth supervision; the MSE formulation in eq.2 is only used for evaluation. The lack of ground-truth information in the presence of noise is precisely what renders this problem very difficult. If any partial ground-truth information is available, then this can be incorporated into the loss function.

## 5 GNNSync architecture

### Obtaining directed graph embeddings

For obtaining digraph embeddings, any digraph GNN that outputs node embeddings can be applied, e.g. DIMPA by He et al. (2022b), the inception block model by Tong et al. (2020), and MagNet by
Zhang et al. (2021). Here we employ DIMPA; details are in A.2. Denoting the final node embedding matrix by \(\mathbf{Z}\in\mathbb{R}^{n\times kd}\), the embedding vector \(\mathbf{z}_{i}\) for a node \(i\) is \(\mathbf{z}_{i}=(\mathbf{Z})_{(i,:)}\in\mathbb{R}^{kd},\) the \(i^{\text{th}}\) row of \(\mathbf{Z}\).

### Obtaining initial estimated angles

To obtain the initial estimated angles for the angular synchronization problem, we introduce a trainable vector \(\mathbf{a}\) with dimension equal to the embedding dimension, then calculate the unnormalized estimated angles by the inner product of \(\mathbf{z}_{i}\) with \(\mathbf{a}\), plus a trainable bias \(b\), followed by a sigmoid layer to force positive values, and finally rescale the angles to \([0,2\pi)\); in short: \(r_{i}^{(0)}=2\pi\text{ sigmoid}(\mathbf{z}_{i}\cdot\mathbf{a}+b)\).

For general \(k\)-synchronization, we apply independent \(\mathbf{a},b\) values to obtain \(k\) different groups of initial angle estimates based on different columns of the node embedding matrix \(\mathbf{Z}\). In general, denote \(\mathbf{Z}_{i,w:\text{x}}\) as the \((v-u+1)\)-vector whose entries are from the \(i\)-th row and the \(u\)-th to \(v\)-th columns of the matrix \(\mathbf{Z}\). With a trainable vector \(\mathbf{a}^{(l)}\) for each \(l\in\{1,\dots,k\}\) with dimension equal to \(d\), we obtain the unnormalized estimated angles by the inner product of \(\mathbf{Z}_{i,(l-1)d+1:ld}\) with \(\mathbf{a}^{(l)}\), plus a trainable bias \(b_{l}\), followed by a sigmoid layer to force positive angle values, then rescale the angles to \([0,2\pi)\); in short: \(r_{i,l}^{(0)}=2\pi\text{ sigmoid}(\mathbf{z}_{i,(l-1)d+1:ld}\cdot\mathbf{a}^{(l)}+b_{l})\).

### Projected gradient steps for final angle estimates

Our final angle estimates are obtained after applying several (default: \(\Gamma=5\)) projected gradient steps to the initial angle estimates. In brief, projected gradient descent for constrained optimization problems first takes a gradient step while ignoring the constraints, and then projects the result back onto the feasible set to incorporate the constraints. Here the projected gradient steps are inspired by Boumal (2016). We construct \(\mathbf{H}\) by \(\mathbf{H}_{i,j}=\exp(\iota\mathbf{A}_{i,j})\mathbb{1}(\mathbf{A}_{i,j}\neq 0),\) and update the estimated angles using Algo. 1, where \(\mathbf{r}_{:,l}\) denotes the \(l\)-th column of \(\mathbf{r}\). In Algo. 1 the gradient step is on line 6, while the projection step on line 7 projects the updated matrix to elementwise angles. Fig. 2 shows the GNNsync framework.

If graph assignments can be estimated effectively right after the GNN, one can replace \(\mathbf{H}\) with \(\mathbf{H}^{(l)}\) for each \(l=1,\dots,k\) separately, where \(\mathbf{H}_{i,j}^{(l)}=\exp(\iota\mathbf{A}_{i,j}^{(l)})\mathbb{1}(\mathbf{A}_ {i,j}^{(l)}\neq 0),\) and \(\mathbf{A}_{i,j}^{(l)}=\mathbb{1}(g(i,j)=1)\mathbf{A}_{i,j}\) is the estimated adjacency matrix for graph \(\mathcal{G}_{l}\) using network assignments from \(g(i,j)\) from eq. (4) applied to the initial angle estimates \(\mathbf{r}^{(0)}\). Yet, separate \(\mathbf{H}^{(l)}\)'s may make the architecture sensitive to the accuracy of graph assignments after the GNN, and hence for robustness we simply choose a single \(\mathbf{H}\). We also find in Sec. 6.4 that the use of \(\mathbf{H}\) instead of separate \(\mathbf{H}^{(l)}\)'s is essential for satisfactory performance. Besides, it is possible to make Algo. 1 parameter-free by further fixing the \(\{\alpha_{\gamma}\}\) values (default: \(\alpha_{\gamma}=1,\forall\gamma\)); we find that using fixed \(\{\alpha_{\gamma}\}\) does not strongly affect performance in our experiments. GNNsync executes the projected gradient descent steps at every training iteration as part of a unified end-to-end training process, but one could also use Algo. 1 to post-process predicted angles without putting the steps in the end-to-end framework. We find that putting Algo. 1 in our end-to-end training framework is usually helpful.

Figure 2: GNNsync overview: starting from an adjacency matrix \(\mathbf{A}\) encoding (noisy) pairwise offsets and an input feature matrix \(\mathbf{X}\), GNNsync first applies a directed GNN to learn node embeddings \(\mathbf{Z}\). It then calculates the inner product with a learnable vector (or \(k\) learnable vectors for \(k>1\)) to produce the initial estimated angles \(r_{i,l}^{(0)}\in[0,2\pi)\) for \(l\in\{1,\dots,k\}\), after rescaling. It then applies several projected gradient steps to the initial angle estimates to obtain the final angle estimates, \(r_{i,l}\in[0,2\pi)\). Let the ground-truth angle matrix be \(\mathbf{R}\in\mathbb{R}^{n\times k}\). The loss function is applied to the output angle matrix \(\mathbf{r}\), given \(\mathbf{A}\), while the final evaluation is based on \(\mathbf{R}\) and \(\mathbf{r}\). Orange frames indicate trainable vectors/matrices, green squares fixed inputs, the red square the final estimated angles (outputs), and the yellow circles the loss function and evaluation.
### Robustness of GNNSync

Measurement noise that perturbs the edge offsets can significantly impact the performance of group synchronization algorithms. To this end, we demonstrate the robustness of GNNSync to such noise perturbations, with the following theoretical guarantee, proved and further discussed in App. A.2

**Proposition 2**.: _For adjacency matrices \(\mathbf{A},\hat{\mathbf{A}}\), assume their row-normalized variants \(\mathbf{A}_{s},\hat{\mathbf{A}}_{s},\mathbf{A}_{t},\hat{\mathbf{A}}_{t}\) satisfy \(\left\|\mathbf{A}_{s}-\hat{\mathbf{A}}_{s}\right\|_{F}<\epsilon_{s}\) and \(\left\|\mathbf{A}_{t}-\hat{\mathbf{A}}_{t}\right\|_{F}<\epsilon_{t},\) where subscripts \(s,t\) denote source and target, resp. Assume further their input feature matrices \(\mathbf{X},\hat{\mathbf{X}}\) satisfy \(\left\|\mathbf{X}-\hat{\mathbf{X}}\right\|_{F}<\epsilon_{f}\). Then their initial angles \(\mathbf{r}^{(0)},\hat{\mathbf{r}}^{(0)}\) from a trained GNNSync using DIMPA satisfy \(\left\|\mathbf{r}^{(0)}-\hat{\mathbf{r}}^{(0)}\right\|_{F}<B_{s}\epsilon_{s}+ B_{t}\epsilon_{s}+B_{f}\epsilon_{f}\), for values \(B_{s},B_{t},B_{f}\) that can be bounded by imposing constraints on model parameters and input._

## 6 Experiments

Implementation details are in App. C and extended results in App. D.

### Data sets and protocol

Previous works in angular synchronization typically only consider synthetic data sets in their experiments, and those applying synchronization to real-world data do not typically publish the data sets. To bridge the gap between synthetic experiments and the real world, we construct synthetic data sets with both correlated and uncorrelated ground-truth rotation angles, using various measurement graphs and noise levels. In addition, we conduct sensor network localization on two data sets.

For synthetic data, we perform experiments on graphs with \(n=360\) nodes for different measurement graphs, with edge density parameter \(p\in\{0.05,0.1,0.15\}\), noise level \(\eta\in\{0,0.1,\ldots,0.9\}\) for \(k=1\), and \(\eta\in\{0,0.1,\ldots,0.7\}\) for \(k\in\{2,3,4\}\). The graph generation procedure is as follows (with further details in App. B.1): *1) Generate \(k\) group(s) of ground-truth angles. One option is to generate each angle from the same Gamma distribution with shape 0.5 and scale \(2\pi\). We denote this option with subscript "1". As angles could be highly correlated in practical scenarios, we introduce a more realistic but challenging option "2", with multivariate normal ground-truth angles. The mean of the ground-truth angles is \(\pi\), with covariance matrix for each \(l\in\{1,\ldots,k\}\) defined by \(\mathbf{ww}^{\top}\), where entries in \(\mathbf{w}\) are generated independently from a standard normal distribution. We explore two more options in the SI. We then apply mod \(2\pi\) to all angles. * 2) Generate a noisy background adjacency matrix \(\mathbf{A}_{\text{noise}}\in\mathbb{R}^{n\times n}\). * 3) Construct a complete adjacency matrix where \(\eta\) portion of the entries are noisy and the rest represent true angular differences. * 4) Generate a measurement graph and sparsify the complete adjacency matrix by only keeping the edges in the measurement graph.

We construct 3 types of measurement graphs from NetworkX (Hagberg et al., 2008) and use the following notations, where the subscript \(o\in\{1,2,3,4\}\) is the option mentioned in step 1) above:

\(\bullet\) Erdos-Renyi (ER) Outlier model: denoted by \(\text{ERO}_{o}(p,k,\eta)\), using as the measurement graph the ER model from NetworkX, where \(p\) is the edge density parameter for the ER measurement graph;

\(\bullet\) Barabasi Albert (BA) Outlier model: denoted by \(\text{BAO}_{o}(p,k,\eta)\), where the measurement graph is a BA model with the number of edges to attach from a new node to existing nodes equal to \([np/2]\), using the standard implementation from NetworkX Hagberg et al. (2008); and

\(\bullet\) Random Geometric Graph (RGG) Outlier model: denoted by \(\text{RGGO}_{o}(p,k,\eta)\), with NetworkX parameter "distance threshold value (radius)" \(2p\) for the RGG measurement graph. For \(k=1\), we omit the value \(k\) and subscript \(o\) in the notation, as the two options coincide in this special case.

For real-world data, we conduct sensor network localization on the U.S. map and the PACM point cloud data set (Cucuringu et al., 2012a) with a focus on the SO(2) component, as follows, with data processing details provided in App. B.2. *1) Starting with the ground-truth locations of \(n=1097\) U.S. cities (resp., \(n=426\) points), we construct patches using each city (resp., point) as a central node and add its \(50\) nearest neighbors to the corresponding patch. *2) For each patch, we add noise to each node's coordinates independently. *3) We then rotate the patches using random rotation angles (ground-truth angles generated as in 1) for synthetic models). For each pair of patches that have at least \(6\) overlapping nodes, we apply Procrustes alignment (Gower, 1975) to estimate the rotation angle based on these overlapping nodes and add an edge to the observed measurement adjacency matrix. *4) We perform angular synchronization to obtain the initial estimated angles and update the estimated angles by shifting by the average pairwise differences between the estimated and ground-truth angles, to eliminate the degree of freedom of a global rotation. *5) Finally, we apply the estimated rotations to the noisy patches and estimate node coordinates by averaging the estimated locations for each node from all patches that contain this node.
### Baselines

In our numerical experiments for angular synchronization, we compare against **7 baselines**, where results are averaged over 10 runs: \(\bullet\) Spectral Baseline (Spectral) by Singer (2011), \(\bullet\) Row-Normalized Spectral Baseline (Spectral_RN) by Cucuringu et al. (2012a), \(\bullet\) Generalized Power Method (GPM) by Boumal (2016), \(\bullet\) TranSync by Huang et al. (2017), \(\bullet\) CEMP_GCW, \(\bullet\) CEMP_MST by Lerman & Shi (2022), and \(\bullet\) Trimmed Averaging Synchronization (TAS) by Maunu & Lerman (2023).

For more general \(k\)-synchronization, we compare against two baselines from Cucuringu & Tyagi (2022), which are based on the top \(k\) eigenvectors of the matrix \(\mathbf{H}\) or its row-normalized version. We use names \(\bullet\) Spectral and \(\bullet\) Spectral_RN to denote them as before. To show that GNNSync (as well as the baselines) deviate from trivial or random solutions, we include an additional baseline denoted "Trivial" for each \(k\), where all angles are predicted equal (with value 1, for simplicity).

### Main experimental results

By default, we use the output angles of the baseline "Spectral_RN" as input features for GNNSync, and thus \(d_{\text{in}}=k\). The main experimental results are shown in Fig. 3 for \(k=1\), and Fig. 4 for general \(k\in\{2,3,4\}\), with additional results reported in App. D. For \(k>1\), we use "GNNSync-cycle", "GNNSync-upset" and "GNNSync-sum" to denote GNNSync variants when considering the training loss function \(\mathcal{L}_{\text{cycle}},\mathcal{L}_{\text{upset}}\), and \(\mathcal{L}_{\text{upset}}+\mathcal{L}_{\text{cycle}}\), respectively.

Figure 4: MSE performance on \(k\)-synchronization for \(k\in\{2,3,4\}\). \(p\) is the network density and \(\eta\) is the noise level. Error bars indicate one standard deviation. Dashed lines highlight GNNSync variants.

Figure 3: MSE performance on angular synchronization (\(k=1\)). Error bars indicate one standard deviation. Dashed lines highlight GNNSync variants.
From Fig. 3 (with additional figures in App. D Fig. 6-8), we conclude that GNNSync produces generally the best performance compared to baselines, in angular synchronization (\(k=1\)). From Fig. 4 (see also App. D Fig. 9-17), we again conclude that GNNSync variants attain leading performance for \(k>1\). The first two columns of Fig. 4 compare the performance of the two options of ground-truth angles on RGGO models. In columns 3 and 4, we show the effect of varying density parameter \(p\), and different synthetic models under various measurement graphs.

For \(k>1\), GNNSync-upset performs better than both baselines in most cases, with \(\mathcal{L}_{\text{upset}}\) simple yet effective to train. GNNSync-cycle generally attains the best performance. As the problems become harder (with increasing \(\eta\), decreasing \(p\), increasing \(k\), more complex measurement graph RGG), GNNSync-cycle outperforms both baselines and other GNNSync variants by a larger margin. The performance of GNNSync-sum lies between that of GNNSync-upset and GNNSync-cycle, but is closer to that of GNNSync-upset, see App. D for more discussions on linear combinations of the two losses. We conclude that while GNNSync-upset generally attains satisfactory performance, GNNSync-cycle is more robust to harder problems than other GNNSync variants and the baselines. Accounting for the performance of trivial guesses, we observe that GNNSync variants are more robust to noise, and attain satisfactory performance even when the competitive baselines are outperformed by trivial guesses. We highlight that there is a clear advantage of using cycle consistency in the pipeline, especially when the problem is harder, thus reflecting the angular nature of the problem. For 3-cycle consistency and the cycle loss \(\mathcal{L}_{\text{cycle}}\), gradient descent in principle drives down the (non-negative) values \(S\) of the sum of three predicted angular differences. To minimize the \(S\) values, we encourage a reweighing process of the initial edge weights so that cycle consistency roughly holds. Unlike \(\mathcal{L}_{\text{upset}}\) which explicitly encourages small \(\mathbf{M}_{i,j}\) values for all edges, \(\mathcal{L}_{\text{cycle}}\) only implicitly encourages small \(\mathbf{M}_{i,j}\) values via the confidence matrix reweighing process for edges with relatively small noise. In an ideal case, we only have large \(\mathbf{M}_{i,j}\) values on noisy edges. In this case, the reweighing process would downweight these noisy edges, which results in a smaller value of the cycle loss function. This is also the underlying reason why \(\mathcal{L}_{\text{cycle}}\) is more robust to noise than \(\mathcal{L}_{\text{upset}}\). For \(k>1\), we hence recommend using the more intricate \(\mathcal{L}_{\text{cycle}}\) function as the training loss function, and we will focus on GNNSync-cycle in the ablation study.

From Fig. 1 (see also App. D Tab. 1-4, Fig. 18-27), we observe that GNNSync is able to align patches and recover coordinates effectively, and is more robust to noise than baselines. GNNSync attains competitive MSE values and Average Normalized Error (ANE) results, where ANE (defined explicitly in App. D.1) measures the discrepancy between the predicted locations and the actual locations.

### Ablation study and discussion

In this subsection, we justify several model choices for all \(k\): \(\bullet\) the use of the projected gradient steps; \(\bullet\) an end-to-end framework instead of training first without the projected gradient steps and then applying Algo. 1 as a post-processing procedure; \(\bullet\) fixed instead of trainable \(\{\alpha_{\gamma}\}\) values. For \(k>1\), we also justify the use of the \(\mathbf{H}\) matrix in Algo. 1 instead of separate \(\mathbf{H}^{(l)}\)'s based on estimated graph assignments of the edges. To validate the ability of GNNSync to borrow strength from baselines, we set the input feature matrix \(\mathbf{X}\) as a set of angles that is estimated by one of the baselines (or \(k\) sets of angles estimated by one of the baselines for \(k>1\)) and report the performance.

Due to space considerations, results for the ablation study are reported in App. D. For \(k=1\), Fig. 22-24 report the MSE performance for different GNNSync variants. Improvements over all possible baselines when taking their output as input features for \(k=1\) are reported in Fig. 34-36. For \(k>1\), we report the results when using \(\mathcal{L}_{\text{cycle}}\) as the training loss function in Fig. 25-33. We conclude that Algo. 1 is indeed helpful in guiding GNNSync to attain lower loss values (we omit loss results for space considerations) and better MSE performance, and that end-to-end training usually attains comparable or better performance than using Algo. 1 for post-processing, even when there is no trainable parameter in Algo. 1. Furthermore, we observe across all data sets, that GNNSync usually improves on existing baselines when employing their outputs as input features, and never performs significantly worse than the corresponding baseline; hence, GNNSync can be used to enhance existing methods. Further, setting \(\{\alpha_{\gamma}\}\) values to be trainable does not seem to boost performance much, and hence we stick to fixed \(\{\alpha_{\gamma}\}\) values. For \(k>1\), using separate \(\mathbf{H}^{(l)}\)'s instead of the whole \(\mathbf{H}\) in Algo. 1 harms performance, which can be explained by the fact that learning graph assignments effectively via GNN outputs is challenging.
Conclusion and outlook

This paper proposed a general NN framework for angular synchronization and a heterogeneous extension. As the current framework is limited to SO(2), we believe that extending our GNN-based framework to the setting of other more general groups is an exciting research direction to pursue, and constitutes ongoing work (for instance, for doing synchronization over the full Euclidean group \(\text{Euc}(2)=\mathbb{Z}_{2}\times\text{SO}(2)\times\mathbb{R}^{2}\)). We also plan to optimize the loss functions under constraints, train our framework with supervision of ground-truth angles (anchor information), and explore the interplay with low-rank matrix completion. Another interesting direction is to extend our SNL example to explore the graph realization problem, of recovering point clouds from a sparse noisy set of pairwise Euclidean distances.