{'fraction_non_alphanumeric': 0.04638423645320197, 'fraction_numerical': 0.015146141215106733, 'mean_word_length': 4.451199427139277, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral Sparsification for Communication-Efficient Collaborative Rotation and Translation Estimation

We propose fast and communication-efficient distributed algorithms for rotation averaging and translation recovery problems that arise from multi-robot simultaneous localization and mapping (SLAM) and distributed camera network localization applications. Our methods are based on theoretical relations between the Hessians of the underlying Riemannian optimization problems and the Laplacians of suitably weighted graphs. We leverage these results to design a distributed solver that performs approximate second-order optimization by solving a Laplacian system at each iteration. Crucially, our algorithms permit robots to employ spectral sparsification to sparsify intermediate dense matrices before communication, and hence provide a mechanism to trade off accuracy with communication efficiency with provable guarantees. We perform rigorous theoretical analysis of our methods and prove that they enjoy (local) linear rate of convergence on the problems of interest. Numerical experiments show that the proposed methods converge to high-precision solutions in a few iterations and that they are significantly more communication-efficient compared to baseline second-order solvers.


Introduction
Multi-robot simultaneous localization and mapping (SLAM) is a fundamental capability that enables consistent spatial perception across a distributed robot team. State-of-the-art systems (e.g., [1][2][3]) typically solve a pose graph optimization (PGO) problem in the estimation back-end, where the robots collaboratively estimate their trajectories in a global frame from noisy relative rotation and translation measurements.
To solve the underlying non-convex optimization problem, recent efforts focus on developing fully distributed algorithms in which robots carry out iterative optimization by passing messages over the communication network [4][5][6][7][8][9]. While these methods are flexible in terms of the required communication architecture, they often suffer from slow convergence due to their first-order nature and the inherent poor conditioning of typical SLAM problems. To resolve the slow convergence issue, an alternative is to pursue a secondorder optimization framework. Cunningham et al. [10][11][12] develop the pioneering work of DDF-SAM where robots exchange marginals over commonly observed landmarks. From an optimization perspective, DDF-SAM achieves second-order updates by performing partial elimination on robots' local Hessian matrices. A shortcoming of this approach is that the resulting matrices are usually dense (even if the original problem is sparse), and hence could be expensive to transmit. Further work has investigated heuristic schemes that replace the dense matrices with sparse approximations [13][14][15][16][17]. However, these methods lack formal performance guarantees, and, in practice, could cause slow convergence or have unpredictable behaviors.
The above discussion reveals a conflict between achieving fast convergence and efficient communication in collaborative SLAM. In this work, we ask the natural question: can we achieve both at the same time? Or to take a step back: can we trade off one performance metric for the other with provable guarantees? We present results towards answering these questions in the affirmative. Specifically, we develop fast and communication-efficient methods for solving the rotation averaging and translation recovery problems that arise from distributed SLAM, camera network localization, and structure-from-motion (SfM) applications, which in practice converge to high-precision solutions in a few iterations. Our approach is based on theoretical relations between the Hessians of the optimization problems and the Laplacians of the underlying graphs. We leverage these theoretical insights to develop a fast collaborative optimization method in which each iteration computes an approximate second-order update by solving a Laplacian system of the form Lx = b. Importantly, during communication, robots use spectral sparsification [18] to sparsify the dense matrices to be transmitted with provable guarantees. By varying the degree of sparsification, our method thus provides a principled way for trading off accuracy with communication efficiency.
Furthermore, the theoretical properties of spectral sparsification allow us to perform rigorous convergence analysis, and establish a linear rate of convergence for the proposed methods. We demonstrate the superior speed and communication efficiency of our methods on large-scale benchmark datasets. In addition, we show that our collaborative rotation averaging and translation recovery solvers produce solutions that closely approximate optimal solutions of the full PGO problem. Together, the empirical results and theoretical guarantees make our approach an appealing choice to implement in real-world multi-robot SLAM systems. Although this work is primarily motivated by collaborative SLAM applications, we believe our results on rotation averaging could be of independent interest to the computer vision community, where rotation averaging plays an important role in applications such as distributed camera networks [4] and large-scale SfM [19,20].
The rest of this paper is organized as follows. In the remainder of this section, we introduce necessary notations and mathematical preliminaries. In Sec. 2, we review related works. In Sec. 3, we formally introduce the rotation averaging and translation recovery problems, and establish theoretical relations between the Hessians and the underlying graph Laplacians. Then, in Sec. 4, we leverage these theoretical results to design fast and communication-efficient solvers for the problems of interest and establish convergence guarantees. Finally, Sec. 5 presents numerical evaluations of the proposed algorithms.

Notations and Preliminaries
Unless stated otherwise, lowercase and uppercase letters are generally used for vectors and matrices, respectively. We define [n] {1, 2, . . . , n} as the set of positive integers from 1 to n.

Linear Algebra and Spectral Approximation
S n and S n + denote the set of n × n symmetric and symmetric positive semidefinite matrices, respectively. For any matrix A, ker(A) and image(A) denote the kernel and image of A. A † denotes the Moore-Penrose inverse of A, which coincides with the inverse A −1 when A is invertible. When A ∈ S n , λ 1 (A), . . . , λ n (A) denote the real eigenvalues of A sorted in increasing order. When A 0, we also define X A tr(X AX) where X is any matrix with compatible dimensions.
Following [21,22], for A, B ∈ S n and > 0, we say that B is an -approximation of A, denoted as A ≈ B, if the following holds, Note that the above relation is symmetric and holds under composition: if A ≈ B and B ≈ δ C, then A ≈ +δ C. Furthermore, if A is singular, the relation (1) implies that B is necessarily singular and ker(A) = ker(B).

Graph Theory
A weighted undirected graph is denoted as G = (V, E, w), where V and E denote the vertex and edge sets, and w : E → R >0 denotes the edge weight function. For a graph G with n vertices, its graph Laplacian L(G; w) ∈ S n + is defined as, if i = j and e = (i, j) ∈ E, 0, otherwise. (2) In (2), E(i) denotes the subset of edges incident to vertex i ∈ V. Our notation L(G; w) serves to emphasize that the Laplacian of G depends on the edge weight w. When the edge weight w is irrelevant or clear from context, we will write the graph as G = (V, E) and its Laplacian as L(G) or simply L. The graph Laplacian L is always positive semidefinite and has a zero eigenvalue, i.e., λ 1 (L) = 0. The second smallest eigenvalue λ 2 (L) is known as the algebraic connectivity, which is always positive for connected graphs.

Riemannian Manifolds
The reader is referred to [23,24] for a comprehensive review of optimization on matrix manifolds. In general, we use M to denote a smooth matrix manifold. For integer n > 1, M n denotes the product manifold formed by n copies of M. T x M denotes the tangent space at x ∈ M. For tangent vectors η, ξ ∈ T x M, their inner product is denoted as η, ξ x , and the corresponding norm is η x = η, η x . In the rest of the paper, we drop the subscript x as it will be clear from context. At x ∈ M, the injectivity radius inj(x) is a positive constant such that the exponential map Exp x : T x M → M is a diffeomorphism when restricted to the domain U = {η ∈ T x M : η < inj(x)}. In this case, we define the logarithm map to be Log x Exp −1 x . Unless otherwise mentioned, we use d(x, y) to denote the geodesic distance between two points x, y ∈ M induced by the Riemannian metric. In addition, it holds that d(x, y) = v where v = Log x (y); see [24,Proposition 10.22].

The Rotation Group SO(d)
The rotation group is denoted as SO(d) = {R ∈ R d×d : R R = I, det(R) = 1}. The tangent space at R is given by T R SO(d) = {RV : V ∈ so(d)}, where so(d) is the space of d × d skew symmetric matrices. In this work, we exclusively work with 2D and 3D rotations. For d = 3, note that T R SO(3) is 3-dimensional. We define a basis for T R SO(3) such that each tangent vector η ∈ T R SO(3) is identified with a vector v ∈ R 3 , Note that (3) defines a bijection between η ∈ T R SO(3) and v ∈ R 3 . For d = 2, we can define a similar basis for the 1-dimensional tangent space T R SO(2). Specifically, each tangent vector η ∈ T R SO(2) is identified by a scalar v ∈ R as, We have overloaded the notation [·] × to map the input scalar or vector to the corresponding skew symmetric matrix in so(2) or so (3). Under the basis given in (3) and (4), the inner product on the tangent space is defined by the corresponding vector dot product, i.e., η 1 , η 2 = v 1 v 2 where v 1 , v 2 ∈ R p are vector representations of η 1 and η 2 , and p = dim SO(d) = d(d − 1)/2. We define the function Exp : R p → SO(d) as, where exp(·) denotes the conventional matrix exponential. Note that Exp :  [10,12] where agents exchange Gaussian marginals over commonly observed landmarks. Paull et al. [13] propose to sparsify the resulting dense information matrix using convex optimization. From a linear algebra perspective, marginalization corresponds to a domain decomposition approach (e.g., see [25,Chapter 14]) where one eliminates a subset of variables using the Schur complement. This idea lies at the heart of state-of-the-art centralized solvers [26][27][28]. Related works use sparse approximations of the resulting matrix (e.g., with tree-based sparsity patterns) to precondition the underlying numerical optimization [14][15][16][17]. Recent work [29] combines domain decomposition with event-triggered transmission to improve communication efficiency during collaborative estimation. Zhang et al. [30] develop a centralized incremental solver for multi-robot SLAM. Fully decentralized solvers for SLAM have also gained increasing attention; see [4][5][6][7][8][9] and also [31] for a recent survey. While fully decentralized methods are the most flexible in terms of the required communication architecture (since they only rely on peer-to-peer message passing), they often have slower convergence due to the first-order nature of the underlying optimization algorithms.

Graph Structure in Rotation Averaging and Pose Graph SLAM
Prior works have investigated the impact of graph structure on rotation averaging and PGO problems from different perspectives. One line of research [32][33][34] adopts an estimation-theoretic approach and shows that the Fisher information matrix is closely related to the underlying graph Laplacian matrix. Eriksson et al. [35] establish sufficient conditions for strong duality to hold in rotation averaging, where the derived analytical error bound depends on the algebraic connectivity of the graph. In a recent work, Bernreiter et al. [36] use tools from graph signal processing to correct onboard estimation errors in multi-robot mapping. Doherty et al. [37] propose a measurement selection approach for pose graph SLAM that seeks to maximize the algebraic connectivity of the underlying graph. This paper differs from the aforementioned works by analyzing the impact of graph structure on the underlying optimization problems, and exploiting the theoretical analysis to design novel optimization algorithms in the multi-robot setting. Among related works in this area, the ones most related to this paper are [38][39][40][41]. Carlone [38] analyzes the influences of graph connectivity and noise level on the convergence of Gauss-Newton methods when solving PGO. In a pair of papers [39,40], Wilson et al. study the local convexity of rotation averaging under the intrinsic (geodesic) distance, by bounding the Riemannian Hessian using the Laplacian of a suitably weighted graph. Recently, Nasiri et al. [41] develop a Gauss-Newton method for rotation averaging under the chordal distance, and show that its convergence basin is influenced by the norm of the inverse reduced Laplacian matrix. Our work differs from [38][39][40][41] by focusing on the development of fast and communicationefficient solvers in multi-robot teams with provable performance guarantees. During this process, we also prove new results on the connections between the Riemannian Hessian and graph Laplacian, and show that they hold under both geodesic (intrinsic) and chordal (extrinsic) distance; see Remark 1.

Spectral Sparsification and Laplacian Solvers
A remarkable property of graph Laplacians is that they admit sparse spectral approximations; see [18] for a survey. Spielman and Srivastava [42] show that every graph with n vertices can be approximated using a sparse graph with O(n log n) edges. This is achieved using a random sampling procedure that selects each edge with probability proportional to its effective resistance, which intuitively measures the importance of each edge to the whole graph. Batson et al. [43] develop a procedure based on the so-called barrier functions for constructing linear-sized sparsifiers. Another line of work [44,45] employs sparsification during approximate Gaussian elimination. Spectral sparsification is one of the main tools that enables recent progress in fast Laplacian solvers (i.e., for solving linear systems of the form Lx = b, where L is graph Laplacian); see [46] for a survey. Peng and Spielman [47] develop a parallel solver that invokes sparsification as a subroutine, which is improved and extended in following works [21,22]. Recently, Tutunov [48] extends the approach in [47] to solve decentralized consensus optimization problems. In this work, we leverage spectral sparsification to design communication-efficient collaborative optimization methods for rotation averaging with provable convergence guarantees; see Remark 4.

Laplacian Systems Arising from Rotation Averaging and Translation Recovery
In state-of-the-art PGO solvers, the cost function often has a separable structure that decouples the rotation and translation terms. For example, SE-Sync [49] considers the formulation, In (6), R i ∈ SO(d) and t i ∈ R d are rotation matrices and translation vectors to be estimated, R ij ∈ SO(d) and t ij ∈ R d are noisy relative rotation and translation measurements, and κ ij , τ ij > 0 are constant measurement weights. Prior works have leveraged the separable structure in (6) to design effective initialization or optimization methods [50,51]. In this section, we also leverage the separability and show that when considering the estimation of rotations and translations as separate problems, the Hessian matrices of the resulting rotation averaging and translation recovery problems are closely related to the Laplacians of suitably weighted graphs. Our results on rotation averaging could be of independent interest to the computer vision community, since the problem plays a fundamental role in applications such as distributed camera networks [4] and large-scale SfM [19,20]. The theoretical relations we establish in this section pave the way for designing fast and communication-efficient solvers in Sec. 4.

Rotation Averaging
We model rotation averaging using an undirected measurement graph G = (V, E). Each vertex i ∈ V = [n] corresponds to a rotation variable R i ∈ SO(d) to be estimated. Each edge (i, j) ∈ E corresponds to a noisy relative measurement of the form, where R i , R j ∈ SO(d) are the latent (ground truth) rotations and R ij ∈ SO(d) is the measurement noise.
In rotation averaging, we aim to estimate the rotations by minimizing the sum of squared measurement residuals. We consider the following generic formulation that permits the use of both geodesic and chordal distance metrics.
Problem 1 (Rotation Averaging). minimize For each edge (i, j) ∈ E, κ ij > 0 is the corresponding measurement weight. The function ϕ is defined as, Note that (9b) corresponds to the rotation part of (6).
In the multi-robot setting, each robot owns a subset of all rotation variables and needs to collaborate with others to solve the full rotation averaging problem; see Fig. 1a for an illustration. Since Problem 1 is nonlinear and non-convex, an iterative optimization framework is required to search for its minimum. Before proceeding, however, one needs to be careful of the inherent gauge-symmetry of rotation averaging: in (8), note that left multiplying each rotation R i ∈ SO(d), i ∈ [n] by a common rotation S ∈ SO(d) does not change the cost function. As a result, each solution R = (R 1 , . . . , R n ) ∈ SO(d) n actually corresponds to an equivalence class of solutions in the form of, From the optimization perspective, the equivalence relation (10) shows that rotation averaging is actually an optimization problem defined over a quotient manifold M = M/ ∼, where M = SO(d) n is called the total space and ∼ denotes the equivalence relation underlying (10); see [24,Chapter 9] for more details. Accounting for the quotient structure is critical for establishing the relation between the Hessian and the graph Laplacian.
In this work, we are interested in applying Newton's method on the quotient manifold M due to its superior local convergence rate. The Newton update can be derived by considering a local perturbation of the cost function. Specifically, let R = {R 1 , . . . , R n } ∈ SO(d) n be our current rotation estimates. For each rotation matrix R i , we seek a local correction to it in the form of Exp(v i )R i , where v i ∈ R p is some vector to be determined and Exp(·) is defined in (5). In (8), replacing each R i with its correction Exp(v i )R i leads to the following local approximation 1 of the optimization problem: In (11), the overall vector v ∈ R pn is formed by concatenating all v i 's. Compared to (8), the optimization variable in (11) becomes the vector v and the rotations R are treated as fixed. Furthermore, we note that the quotient structure of Problem 1 gives rise to the following vertical space [24,Chapter 9.4] that summarizes all directions of change along which (11) is invariant, Intuitively, N captures the action of any global left rotation. Indeed, for any v ∈ N , we have Exp(v i ) = Exp(v j ) for all i, j ∈ [n], and thus the cost function (11) remains constant. Let us denote the gradient and Hessian of (11) as follows, Our notations g(R) and H(R) serve to emphasize that the gradient and Hessian are defined in the total space M and depend on the current rotation estimates R. In [24,Chapter 9.12], it is shown that executing the Newton update on the quotient manifold amounts to finding the minimum norm solution of the following linear system, where P H is the orthogonal projection onto the horizontal space H, defined as the orthogonal complement of the vertical space, i.e., H N ⊥ . Intuitively, including P H in (14) accounts for the gauge symmetry by eliminating the effect of any vertical component from v. The following theorem reveals an interesting connection between H(R) defined in (14) and the Laplacian of the underlying graph.
Theorem 1 (Local Hessian Approximation for Rotation Averaging). Let R ∈ SO(d) n denote the set of ground truth rotations from which the noisy measurements R ij are generated according to (7). For any δ ∈ (0, 1/2), there exist constantsθ, r > 0 such that if, In (16), G = (V, E) is the measurement graph. The edge weight function w is given by w = κ for the squared geodesic distance cost (9a), and w = 2κ under the squared chordal distance cost (9b).
We prove Theorem 1 in Appendix A, and perform numerical validation in Sec. 5.1. While prior works [39,40] have proved related results for analyzing the local convexity of rotation averaging, we differ from [39,40] by using Theorem 1 to design fast optimization algorithms for both single and multi-robot settings. Before proceeding, we note that Theorem 1 directly implies the following bound on the Hessian H(R).

Corollary 1 (Local Hessian Bound and Condition Number for Rotation Averaging)
. Under the assumptions of Theorem 1, define constants µ H = e −δ λ 2 (L(G; w)) and L H = e δ λ n (L(G; w)). Then for all R ∈ B r (R ), In the following, κ H = L H /µ H is referred to as the condition number.
Theorem 1 shows that under small measurement noise, the Hessian near a global minimizer is well approximated by the Laplacian of an appropriately weighted graph. 2 This result directly motivates an approximate Newton method that replaces the Hessian with its Laplacian approximation (see Remark 1 for connections with prior work). Specifically, instead of solving (14), one solves the following approximate Newton system, In the following, it would be more convenient to consider the matrix form of the above linear system. For this purpose, let us define matrices V, B(R) ∈ R n×p , Using properties of the Kronecker product, we can show that (18) is equivalent to, Algorithm 1 shows the complete pseudocode of the approximate Newton algorithm. Compared to the original Newton's method, Algorithm 1 uses a constant matrix across all iterations, and hence could be significantly more computationally efficient. For this reason, we believe that Algorithm 1 could be of independent interest for standard (centralized) rotation averaging. Furthermore, in Sec. 4, we show that Algorithm 1 admits communication-efficient extensions in multi-robot settings. Compute approximate Newton update by solving L(G; w)V k = B(R k ).

3:
Update iterate by R k+1 Remark 1 (Connections with prior work [41]). In [41], Nasiri et al. first developed Algorithm 1 for the special case of chordal rotation averaging using a Gauss-Newton formulation. In contrast, we motivate Algorithm 1 by proving the theoretical approximation relation between the Hessian and the graph Laplacian (Theorem 1) that holds under both geodesic and chordal distances. Furthermore, in Sec. 4 we extend Algorithm 1 to a fast and communication-efficient solver in the multi-robot context. Lastly, the theoretical approximation relation we establish also allows us to prove local linear convergence for our methods.
Remark 2 (Feasibility of the approximate Newton system). Using the properties of the graph Laplacian and the Kronecker product, we see that ker(L(G; w) ⊗ I p ) = N where N is the vertical space defined in (12). Furthermore, in [24,Chapter 9.8], it is shown that g(R) ⊥ N . Thus, we conclude that g(R) ∈ image(L(G; w) ⊗ I p ), i.e., the linear system (18) and its equivalent matrix form (20) are always feasible. In fact, the system is singular and hence admits infinitely many solutions. Similar to the original Newton's method on quotient manifold, we will select the minimum norm solution v which guarantees that v ⊥ N .

Translation Recovery
It is a known result that the translation variables in PGO can be recovered from a given rotation estimate R ∈ SO(d) n by solving a linear least squares problem [49][50][51]. From (6), we define the translation recovery problem to be, Similar to rotation averaging, Problem 2 can be modeled using the undirected measurement graph G = (V, E), where each vertex represents a translation vector to be estimated and each edge represents a relative measurement between two translations. It can be shown that (21) is equivalent to a linear system involving the graph Laplacian L(G; τ ), where edges are weighted by the translation measurement weights τ . Specifically, denote M t = t 1 . . . t n ∈ R n×d as the matrix where each row corresponds to a translation vector to be estimated. One can show that the optimal translations are solutions of, where B t ∈ R n×d is a constant matrix that only depends on R and t ij . Furthermore, each column of B t belongs to the image of the Laplacian L(G; τ ), so (22)  To conclude this section, we note that similar to rotation averaging, translation recovery (Problem 2) is subject to a gauge symmetry. Specifically, two translation solutions M t and M t are equivalent if they only differ by a global translation. Mathematically, this means that M t = M t + 1 n c where 1 n ∈ R n is the vector of all ones and c ∈ R d is some constant vector.

Algorithms and Performance Guarantees
In Sec. 3, we have shown that Laplacian systems naturally arise when solving the rotation averaging and translation recovery problems; see (20) and (22), respectively. Recall that we seek to find the solution X ∈ R n×p to a linear system of the form,  where L ∈ S n + is the Laplacian of a graph that involves multiple robots' variables, and each column of B ∈ R n×p is in the image of L so that (23) is always feasible. In Sec. 4.1, we develop a communicationefficient distributed solver for (23). Then, in Sec. 4.2 and Sec. 4.3, we use the developed solver to design communication-efficient algorithms for collaborative rotation averaging and translation recovery, and establish rigorous theoretical guarantees for both cases. Proofs of all results in this section are provided in Appendices B-C.

A Collaborative Laplacian Solver with Spectral Sparsification
We propose to solve (23) using the domain decomposition framework [25,Chapter 14], which has been utilized in earlier works such as DDF-SAM [10][11][12] to solve collaborative SLAM problems. This is motivated by the fact that in a collaborative SLAM scenario with m robots, there is a natural disjoint partitioning of the vertex set V: where V α contains all vertices (variables) of robot α ∈ [m]. Furthermore, V α can be partitioned as V α = F α C α where C α denotes all separator (interface) vertices and F α denotes all interior vertices of robot α. In multi-robot SLAM, the separators are given by the set of variables that have inter-robot measurements; see Fig. 1a. Note that given the set of all separators C = C 1 . . . C m , robots' interior vertices F α become disconnected from each other. The natural vertex partitioning in (24) further gives rise to a disjoint partitioning of the edge set, For each robot α ∈ [m], its local edge set E α consists of all edges that connect two vertices from V α . In Fig. 1a, the local edges are shown using colors corresponding to the robots. The remaining inter-robot edges form E c , which are highlighted as bold black edges in Fig. 1a.
In domain decomposition, we adopt a variable ordering in which the interior nodes F = F 1 . . . F m appear before the separators C = C 1 . . . C m . With this variable ordering, the Laplacian system (23) can be rewritten as,  For α ∈ [m], X α and B α denote the rows of X and B in (23) that correspond to robot α's interior variables F α . On the other hand, we treat separators from all robots as a single block C = C 1 . . . C m , and use X c and B c to denote the rows of X and B that correspond to the separator block C. (26)). In the multi-robot context, the overall Laplacian system (26) is stored distributedly across the robot team. Specifically, since each robot α knows the subgraph induced by its own vertices V α ( e.g., in Fig. 1a, robot 2 knows all edges incident to the blue vertices), it independently computes and stores its Laplacian blocks L αα and L αc . Similarly, each robot α also independently computes and stores the block B α . Meanwhile, we assume that the blocks defined over separators L cc and B c are handled by a special "server" that performs additional computations in the proposed methods. This role can be performed by a leading robot, or can be replicated on all robots as in DDF-SAM [10].

Remark 3 (Distributed computation of
In practice, the set of separators is usually smaller than the set of interior vertices. This fact together with the "arrowhead" sparsity pattern in (26) motivates us to first solve the reduced system defined over the separators,  In the following, let us define . Then, the matrix on the right-hand side of (27) can be written as, Meanwhile, the matrix S defined on the left-hand side of (27) is the Schur complement resulting from eliminating all interior nodes F from the full Laplacian matrix L, denoted as S = Sc(L, F). The next lemma shows S is the sum of multiple smaller Laplacian matrices.
as the graph induced by inter-robot loop closures E c . Then, the matrix S that appears in (27) can be written as, Since Laplacian matrices are closed under Schur complements [21,Fact 4.2], each S α defined in Lemma 1 is also a Laplacian matrix. 3 Furthermore, as a result of Remark 3, each robot α can independently compute S α = Sc(L(G α ), F α ) and U α = L cα L −1 αα B α . This observation motivates a distributed method in which robots first transmit their S α and U α to the server in parallel. Upon collecting S α and U α from all robots, the server can then form S using (29) and U using (28). It then solves the linear system SX c = U (27) and broadcasts the solution X c back to all robots. Finally, once robots receive the separator solutions X c , they can in parallel recover their interior solutions via back-substitution, The aforementioned method is a distributed implementation of domain decomposition. While it effectively exploits the separable structure in the problem, this method can incur significant communication cost as it requires each robot α to transmit its Schur complement matrix S α that is potentially dense. This issue is illustrated in Fig. 1b, where for robot 2 (blue) its S α corresponds to a dense graph over its separators.
In the following, we propose an approximate domain decomposition algorithm that is significantly more communication-efficient while providing provable approximation guarantees. Our method is based on the facts that (i) each local Schur complement S α is itself a graph Laplacian, and (ii) graph Laplacians admit sparse spectral sparsifications [18], i.e., for a given approximation threshold > 0, one can compute a sparse for each robot α in parallel do 3: Compute a sparse approximation Sα such that Sα ≈ Sα.

7:
Server solves SXc = U (where S is obtained from Algorithm 2), and broadcasts solution Xc to all robots.

8:
for each robot α in parallel do

9:
Compute interior solution via back-substitution Xα = L −1 αα (Bα − LαcXc) . Laplacian S α such that S α ≈ S α . In this work, we implement the sampling approach of [42], which produces S α with O(|C| log |C|) entries (compared to the original dense S α with O(|C| 2 ) entries). Fig. 1c illustrates a spectral sparsification for robot 2's dense reduced graph. In the proposed method, each robot transmits its sparse approximation S α instead of the original Schur complement S α . By summing together these S α matrices, the server can obtain a sparse approximation to the original dense Schur complement S; see Algorithm 2. Afterwards, we can follow the same procedure as standard domain decomposition to obtain an approximate solution to the Laplacian system (23); see Algorithm 3. Specifically, the server first solves an approximate reduced system using S obtained from Algorithm 2 (line 7). Then, the interior solution for each robot is recovered using back-substitution (line 9).
Together, Algorithms 2 and 3 provide a distributed and parallel procedure for computing an approximate solution to the original Laplacian system (23). The methods are communication-efficient, as each robot only needs to transmit a sparse matrix S α (line 4 in Algorithm 2) that has O(|C| log |C|) entries, together with a block vector U α ∈ R |C|×p (line 4 in Algorithm 3) that scales linearly with the number of separators |C|. Crucially, the use of spectral sparsifiers allows us to establish theoretical guarantees on the accuracy of the approximate solution, as summarized in the following theorem.
Theorem 2 (Approximation guarantees of Algorithms 2 and 3). Given a Laplacian system LX = B, Algorithms 2 and 3 together return a solution X ∈ R n×p such that L X = B, where L ∈ S n + satisfies, Furthermore, let X ∈ R n×p be an exact solution to the input linear system, i.e., LX = B. It holds that, where the constant c( ) is defined as, In Theorem 2, we show that the approximate solution X produced by Algorithms 2 and 3 remains close to the exact solution X when measured using the "norm" induced by the original Laplacian L. 4 Server computes B(R k )c that corresponds to all separator variables. 6: for each robot α in parallel do 7: Compute B(R k )α that corresponds to robot α's interior variables Fα. Solve V k = SparsifiedLaplacianSolver(L, B(R k ), S).

11:
// Distributed updates of all rotation variables 12: for each robot α in parallel do 13: Update iterates by R k+1 Furthermore, the quality of the approximation is controlled by the sparsification parameter through the function c( ) visualized in Fig. 2. Note that when = 0, sparsification is effectively skipped and robots transmit the original dense matrices S α . In this case, we have c( ) = 0 and the solution X produced by our methods is exact, i.e., L X = B. Meanwhile, by increasing , our methods smoothly trade off accuracy with communication efficiency.
Remark 4 (Connections with existing Laplacian solvers [21,22]). Our distributed Laplacian solver (Algorithms 2 and 3) is inspired by the centralized solvers developed in [21,22] for solving Laplacian systems in nearly linear time. However, our result differs from these works by focusing on the use of spectral sparsification in the multi-robot setting to achieve communication efficiency. Furthermore, in Sec. 4.2, we apply our distributed Laplacian solver to solve the non-convex Riemannian optimization problem underlying rotation averaging, and establish provable convergence guarantees for the resulting Riemannian optimization algorithm.

Collaborative Rotation Averaging
In this section, we present a fast and communication-efficient distributed solver for rotation averaging. Recall the centralized method in Algorithm 1, where each iteration solves a Laplacian system LV = B(R); see line 2. In the multi-robot setting, we can use the solver developed in Sec. 4.1 to obtain an approximate solution to this system. Algorithm 4 shows the pseudocode. First, an initial guess R 0 is computed (line 1) using existing distributed initialization techniques such as [5]. Then, at line 2, robots first form the approximate Schur complement S using SparsifiedSchurComplement (Algorithm 2). Each iteration consists of three main steps. At the first step (line 4-8), robots compute and store the right-hand side B(R) in a distributed fashion. Specifically, recall from Remark 3 that the overall B(R) is divided into multiple blocks, In our algorithm, each robot α ∈ [m] computes the block B(R) α corresponding to its interior variables F α , and the server computes the block B(R) c corresponding to all separators. At the second step (line 10), robots collaboratively solve for the update vector V k by calling SparsifiedLaplacianSolver (Algorithm 3). Finally, at the last step (line [11][12][13][14], we obtain the next iterate using the solutions V k , where robots in parallel update the rotation variables they own. L, ||X − X|| L ignores any component of X − X that lives on the kernel of L, which is spanned by the vector of all ones 1n. However, this does not create a problem for us since we only seek to compare X and X when considering both as solutions to the Laplacian system LX = B, and using · L naturally eliminates any difference on ker(L) that is inconsequential. Figure 3: Intuitions behind the convergence rate in Theorem 3. Recall from Theorem 1 that under bounded measurement noise, the original Newton system (left box) is locally δ-approximated by a linear system specified by a Laplacian L (middle box). In addition, in Theorem 2 we have shown that our distributed Laplacian solver approximates L with L where L ≈ L (right box). The composition of the two approximation relations thus gives H(R) ≈ δ+ ( L ⊗ Ip), which intuitively explains why the convergence rate in (36) depends on a function of δ + .
Remark 5 (Communication efficiency of Algorithm 4). In Algorithm 4, the majority of communication is induced by calls to SparsifiedSchurComplement and SparsifiedLaplacianSolver. However, note that we only require a single call to SparsifiedSchurComplement for constructing the sparse Schur complement S at the server (line 2). Consequently, each robot only needs to upload its sparse S α matrix once during the entire optimization process. Then, at each iteration, a single call to SparsifiedLaplacianSolver is made (line 10), which requires each robot to upload its block vector U α ∈ R |C|×p that scales linearly with the number of separators |C|. Together, these facts make Algorithm 4 very communication-efficient.
In the following, we proceed to establish theoretical guarantees for our collaborative rotation averaging algorithm. We will show that starting from a suitable initial guess, Algorithm 4 converges to a global minimizer at a linear rate. One might be tempted to state the linear convergence result on the total space, i.e., d(R k+1 , R ) ≤ γ d(R k , R ) where k is the iteration number, γ ∈ (0, 1) is a constant, and R is a global minimizer. However, it is challenging to prove this statement due to the gauge symmetry of rotation averaging. The iterates {R k } might converge to a solution R ∞ that is only equivalent to R up to a global rotation, i.e., and as a result d(R ∞ , R ) = 0 in general. Fortunately, this issue can be resolved using the machinery of Riemannian quotient manifolds. Instead of measuring the distance on the total space d(R k , R ), we will compute the distance between the underlying equivalence classes d ( (33). Under the assumptions of Theorem 1, suppose is selected such that γ(δ + ) < 1. In addition, suppose at each iteration k, the update vector v k ∈ R pn is orthogonal to the vertical space, i.e., v k ⊥ N . Let R be an optimal solution to Problem 1. There exists r > 0 such that for any R 0 where d([R 0 ], [R ]) < r , Algorithm 4 generates an infinite sequence {R k } where the corresponding sequence of equivalence classes [R k ] converges linearly to [R ]. Furthermore, the convergence rate factor is, Theorem 3 shows that using the distance metric on the quotient manifold, Algorithm 4 locally converges to the global minimizer at a linear rate. 5 Fig. 3 provides intuitions behind the convergence rate in (36).
Recall that δ appears in Theorem 1 where we show H(R) ≈ δ (L⊗I p ) under bounded measurement noise. On the other hand, is the parameter for spectral sparsification and is controlled by the user. In Theorem 2, we showed that our methods transform the input Laplacian L into an approximation L such that L ≈ L. The composition of the two approximation relations thus gives H(R) ≈ δ+ ( L ⊗ I p ), which intuitively explains why the convergence rate depends on a function of δ + . Lastly, we note that while Theorem 3 provides rigorous theoretical justification for Algorithm 4, the requirement γ(δ + ) < 1 still leads to conservative Server computes E k c that corresponds to all separator variables.

6:
for each robot α in parallel do 7: Compute E k α that corresponds to robot α's interior variables Fα. for each robot α in parallel do

Collaborative Translation Recovery
Similar to rotation averaging, we can develop a fast and communication-efficient method to solve translation recovery, which is equivalent to the Laplacian system (22) as shown in Sec. 3.2. Specifically, we employ our distributed Laplacian solver (Sec. 4.1) in an iterative refinement framework. Let M k t ∈ R n×d be our estimate for the translation variables at iteration k (in practice M 0 t can simply be initialized at zero). We attempt to find a correction D k to M k t by solving the residual system corresponding to (22): Observing that the system on the right-hand side of (37) is another Laplacian system in L ≡ L(G; τ ), we can deploy our distributed Laplacian solver to find an approximate solution D k . Algorithm 5 shows the pseudocode, which shares many similarities with the proposed collaborative rotation averaging method Algorithm 4. In particular, the computation of the right-hand side E k (line 4-8) and the update step (line [11][12][13][14] are performed in a distributed fashion. The two methods also share the same requirements in terms of communication; see Remark 5. The following theorem states the theoretical guarantees for Algorithm 5. Theorem 4 (Convergence rate of Algorithm 5). Suppose is selected such that the constant c( ) defined in (33) satisfies c( ) < 1. Let M t be an optimal solution to Problem 2 and let M k t denote the solution computed by Algorithm 5 at iteration k ≥ 1. It holds that, where L ≡ L(G; τ ).
Theorem 4 is simpler compared to its counterpart for rotation averaging (Theorem 3). The convergence rate (38) only depends on the sparsification parameter . Furthermore, since the translation recovery problem is convex, the convergence guarantee is global and holds for any initial guess. 6 While Theorem 4 requires c( ) < 1, our experiments in Sec. 5.2 show that Algorithm 5 is not sensitive to the choice of sparsification parameter and converges under a wide range of parameter settings.

Experimental Results
We perform extensive evaluations of our proposed collaborative rotation averaging and translation recovery methods using synthetic and benchmark datasets. Unless otherwise mentioned, we simulate scenarios with

Evaluation of Collaborative Rotation Averaging
In this subsection, we validate our theoretical results and proposed collaborative solver (Algorithm 4) for rotation averaging using benchmark pose graph SLAM datasets. Unless otherwise mentioned, we initialize Algorithm 4 using the distributed chordal initialization approach in [5], where the number of iterations is limited to 50. Our experiments mainly consider rotation averaging problems under the chordal distance metric. Appendix E provides additional results using the geodesic distance.
Numerical validation of Hessian approximation (Theorem 1). First, we perform numerical validation of the approximation relation established in Theorem 1, which serves as the basis of our proposed method. Recall that Theorem 1 states that the Hessian is well approximated by a Laplacian (in the sense of H(R) ≈ δ L ⊗ I p ) locally around a global minimizer. To verify this claim, we generate synthetic chordal rotation averaging problems over a 3D grid with 125 total rotations (Fig. 4a). With a probability of 0.3, we generate noisy relative measurements between pairs of nearby rotations. All rotation measurements are corrupted by increasing level of Langevin noise [49, Appendix A]. At each noise level, we generate 20 random problem instances. For each problem instance, we obtain the global minimizer R (global optimality is certified using the approach in [35]) and numerically compute the smallest constant δ such that H(R ) ≈ δ L ⊗ I p . Fig. 4b shows the evolution of the empirical approximation constant δ as a function of noise level. In the special case when there is no noise, δ is always zero. This result is expected as the global minimizer R coincides with the latent (ground truth) rotations and it can be shown that H(R ) = L ⊗ I p . In general, the empirical value of δ increases smoothly as the noise level increases. Since the Hessian H(R) varies smoothly with R, our results confirm that the Laplacian is a good approximation of the Hessian locally around R , as predicted by Theorem 1.
Evaluation on benchmark datasets. In this experiment, we evaluate Algorithm 4 on rotation averaging problems from benchmark SLAM datasets. Each dataset is divided to simulate a scenario with 5 robots. We first perform an in-depth evaluation of convergence rate and communication efficiency using the Cubicle dataset. Then, we report the comprehensive performance of Algorithm 4 on all benchmark datasets (Table 1). To start, we demonstrate the effectiveness of spectral sparsification on the Cubicle dataset. Recall   that Algorithm 4 calls the SparsifiedSchurComplement procedure (Algorithm 2), which requires each robot α to transmit its sparsified matrix S α . Fig. 5a shows the number of nonzero entries in S α as a function of the sparsification parameter . Note that when = 0, sparsification is effectively skipped and each robot transmits its exact S α matrix that is potentially large and dense. In Fig. 5a, this is reflected on robot 1 (blue curve) whose exact S α matrix has more than 2 × 10 4 nonzero entries and hence is expensive to transmit. However, spectral sparsification significantly reduces the density of the matrix and hence improves communication efficiency. In particular, for robot 1, applying sparsification with = 2 creates a sparse S α with 2300 nonzero entries, which is almost an order of magnitude sparser than the original S α . Next, we evaluate the convergence rate and communication efficiency of Algorithm 4 with varying sparsification parameter . We introduce three baseline methods for the purpose of comparison. The first baseline, called Newton in Fig. 5, implements the exact Newton update using domain decomposition, where each robot computes and transmits the Schur complement of its local Hessian matrix to the server. Inspired by existing works [14][15][16][17], we also implement two baselines that apply heuristic sparsification to Newton: in Block-Diagonal, each robot only transmits the diagonal blocks of its Hessian Schur complement (this strategy is also known as Jacobi preconditioning [25]), whereas in Block-Tree, each robot transmits both diagonal blocks and off-diagonal blocks that form a tree sparsity pattern. Fig. 5b shows the accuracy achieved by all methods (measured by norm of the Riemannian gradient) as a function of iterations. As expected, Newton achieves the best convergence speed and converges to a high-precision solution in two iterations. However, when combined with heuristic sparsifications in Block-Diagonal and Block-Tree, the resulting methods have very slow convergence. Intuitively, this result shows that a diagonal or tree sparsity pattern is not sufficient for preserving the spectrum of the original dense matrix. 7 In contrast, our proposed method achieves fast convergence under a wide range of sparsification parameter . Furthermore, by varying , the proposed method provides a principled way to trade off convergence speed with communication efficiency. Fig. 5c visualizes the accuracy as a function of total uploads to the server. Since both the Hessian and Laplacian matrices are symmetric, we only record the communication when uploading their upper triangular parts as sparse matrices. To convert the result to kilobyte (KB), we assume each scalar is transmitted in double precision. Our results show that the proposed method achieves the best communication efficiency under various settings of the sparsification parameter . Moreover, even without sparsification (i.e., = 0), the proposed method is still more communication-efficient than Newton. This result is due to the following reasons. First, since the Hessian matrix varies across iterations, Newton requires communication of the updated Hessian Schur complements at every iteration. In contrast, the proposed method works with a constant graph Laplacian, and hence only requires a one-time communication of its Schur complements; see line 2 in Algorithm 4. Second, Newton requires communication to form the Schur complement of the original Table 1: Evaluation on 5-robot chordal rotation averaging problems from benchmark pose graph SLAM datasets. |V| and |E| denote the total number of rotation variables and measurements, respectively. We run the baseline Newton method and the proposed method with varying sparsification parameter , and show the number of iterations and total uploads used by each method to reach a Riemannian gradient norm of 10 −5 . Results are averaged over 5 runs. pn-by-pn Hessian matrix, where n is the number of rotation variables and p = dim SO(d) is the intrinsic dimension of the rotation group (for the Cubicle dataset, n = 5750 and p = 3). In contrast, the proposed method operates on the smaller n-by-n Laplacian matrix, and the decrease in matrix size directly translates to communication reduction. Lastly, we extend our evaluation to rotation averaging problems from 12 benchmark pose graph SLAM datasets. For each problem, we simulate a scenario with 5 robots, and run the proposed method (Algorithm 4) with varying sparsification parameter as well as the baseline Newton method. All methods are terminated when the Riemannian gradient norm is smaller than 10 −5 . Since the spectral sparsification method we use [42] is randomized, we perform 5 random runs of our method for each setting of . Table 1 shows the average number of iterations and total uploads used by all methods to reach the desired precision. On all datasets, we are able to verify that all methods converge to the global minima of the considered rotation averaging problems. The proposed method achieves an empirical convergence speed that is close to Newton and typically converges in a few iterations. 8 We note that this is significantly faster than existing fully distributed methods (e.g., [6,8]) that often require hundreds of iterations to achieve a moderate precision. On all datasets, the proposed method is more communication-efficient than the baseline Newton method. We observe that the benefit of sparsification varies across datasets. For example, on Killian Court and INTEL, the effect of sparsification is limited because the exact Schur complement S is small or already sparse. Meanwhile, on datasets such as Grid and Rim, sparsification significantly improves communication efficiency. In particular, on Grid, setting = 1.5 achieves a 89% communication reduction compared to the case without sparsification ( = 0). Overall, we find that sparsification tends to be more effective when the input measurement graph has a large average degree (i.e., the ratio |E|/|V| between the number of measurements and variables is large), since, intuitively, the resulting Schur complements are more likely to be dense.

Datasets
Scalability with number of robots. In this experiment, we evaluate the scalability of Algorithm 4. For this purpose, we generate synthetic rotation averaging problems with increasing number of robots where each robot owns 125 rotation variables. All relative measurements are corrupted by Langevin noise with a standard deviation of 5 deg. Fig. 6a shows an example simulation with 9 robots. We run Algorithm 4 with sparsification parameter = 0.5 until the Riemannian gradient norm reaches 10 −5 . Fig. 6b shows the evolution of gradient norm as a function of iterations. Note that all curves in Fig. 6b have nearly constant slopes, which suggests that the empirical convergence rate of our method is not sensitive to the number     of robots. This observation is compatible with our theoretical analysis in Theorem 3, which shows that the local convergence rate of Algorithm 4 only depends on the Hessian approximation constant δ (which is mainly influenced by noise level as shown in Fig. 4b) and the user-defined sparsification parameter . This property makes our method more appealing than existing fully distributed methods, whose convergence speed typically degrades as the number of robots increases (e.g., see [6, Fig. 8]). On the other hand, the communication costs shown in Fig. 6c does increase, since increasing number of robots leads to larger total uploads at each iteration. Convergence basin. So far, we have used the distributed chordal initialization technique [5] to initialize Algorithm 4. In the next experiment, we test the robustness of our proposed method to poor initial guesses. For this purpose, we use a 9-robot simulation similar to Fig. 6a where each robot owns 512 rotation variables, and generate synthetic initial guesses by perturbing the global minimizer with increasing level of Langevin noise. Using the synthetic initialization, we run Algorithm 4 with sparsification parameter = 0.5 until the Riemannian gradient norm reaches 10 −5 or the number of iterations exceeds 50. At each noise level,   Sparsification runtime. We conclude this subsection by evaluating the runtime of our spectral sparsification implementation based on effective resistance sampling [42]. Recall that in the SparsifiedSchur-Complement step in Algorithm 4, each robot α sparsifies its S α matrix and transmits the result S α to the server. This step uses the majority of robots' local computation time. In Fig. 8, we evaluate the sparsification runtime on the 12 benchmark datasets shown in Table 1. For each dataset, we record the maximum sparsification time among all robots, and visualize the result as a function of the number of nonzero entries in the input matrix S α . On most datasets, the maximum runtime is below one second. On the Grid dataset, the input matrix has more than 5 × 10 5 nonzero entries and our implementation uses 10.2 seconds. Overall, we conclude that the runtime of our implementation is still reasonable. However, we believe that further improvements are possible, e.g., by approximately computing effective resistances during spectral sparsification as suggested in [42].

Evaluation of Collaborative Translation Recovery
In this subsection, we evaluate the proposed collaborative translation recovery algorithm (Algorithm 5) on benchmark pose graph SLAM datasets. In the following experiments, the rotation estimate R in (21) is fixed to be the solution of the rotation averaging subproblem obtained using Algorithm 4.
Evaluation on the Cubicle dataset. Similar to Sec. 5.1, we evaluate the convergence rate and communication efficiency of Algorithm 5 on a 5-robot problem simulated using the Cubicle dataset. Note that for translation recovery, setting = 0 (i.e., no sparsification) in Algorithm 5 recovers the exact solution in Table 2: Comparison to optimal pose graph SLAM solutions. |V| and |E| denote the total number of pose variables and measurements, respectively. We report the rotation and translation root-mean-square error (RMSE) of our approach with respect to optimal solutions computed by SE-Sync [49]. In addition, we also show the number of iterations and total uploads used by Algorithms 4 and 5. a single iteration. In Fig. 9b, the vertical line shows the total uploads incurred with this setting. We run Algorithm 5 under varying sparsification parameter . Similar to the rotation averaging experiment, we introduce two baselines in which each robot heuristically sparsifies its Schur complement matrix S α according to a diagonal sparsity pattern (Diagonal) or a tree sparsity pattern (Tree). Our results in Fig. 9 show that the proposed method dominates the baselines in terms of both iteration and communication complexity. Furthermore, spectral sparsification provides a way to trade off accuracy with communication efficiency. For instance, to obtain an approximate solution, we can run Algorithm 5 with = 0.4 and terminate when the gradient norm reaches 10 −1 , which would incur less communication compared to recovering the exact solution with = 0. Comparison to optimal PGO solutions. In this work, we considered a decoupled approach for solving PGO, where we first estimate the rotations by solving the rotation averaging problem (Algorithm 4), and then estimate the translations by solving translation recovery (Algorithm 5). This approach produces an approximate solution to the full PGO problem (6), since it disregards the coupling between the rotation and translation components. In this experiment, we present quantitative comparisons between solutions computed by this decoupled approach and optimal PGO solutions obtained from SE-Sync [49]. Specifically, we compute the rotation error between our solution R apx ∈ SO(d) n and the optimal solution R opt ∈ SO(d) n from SE-Sync as follows,

Datasets
Intuitively, (39) computes the root-mean-square error (RMSE) between two sets of rotations after alignment by a global rotation. The optimal alignment S in (39) has a closed-form expression; see [49,Appendix C.1].
Similarly, for translations, we report the RMSE between our solution t apx and the optimizer t opt after a global translation alignment. On each benchmark dataset, we simulate a scenario with 5 robots and run Algorithms 4 and 5 with sparsification parameter = 0.5. Both algorithms are terminated when the gradient norm reaches 10 −5 . Table 2 reports the final RMSE together with the number of iterations and total uploads used by our algorithms. We find that our approximate solutions remain very close to the optimal PGO solutions from SE-Sync on most datasets. This observation agrees with prior work [50], where the authors remarked that the translation measurements in pose graph SLAM usually have minor influence on rotation estimation. Meanwhile, on the Killian Court and Rim datasets where the errors are larger, our solutions could serve as high quality initial guesses for further optimization (e.g., using Newton or Gauss-Newton method). Overall, we believe the accuracy together with the computation and communication efficiency of our approach makes it an appealing alternative to state-of-the-art fully distributed solvers (e.g., [6,8]) that suffer from slow convergence.

Conclusion
In this work, we developed fast and communication-efficient methods for solving rotation averaging and translation recovery in multi-robot SLAM and distributed camera network applications. Our algorithms leverage theoretical relations between the Hessians of the optimization problems and the Laplacians of the underlying graphs. At each iteration, robots perform approximate second-order update by collaboratively solving a Laplacian system using the domain decomposition approach. Crucially, spectral sparsification is used to sparsify intermediate dense matrices during communication, which provides an effective way to trade off accuracy with communication efficiency. We perform rigorous analysis of our methods and prove that they achieve (local) linear rate of convergence. Extensive experiments validate our theoretical results and demonstrate the superior convergence rate and communication efficiency of our proposed methods. It remains an interesting open problem whether a similar communication-efficient solver can be developed for the full PGO problem. Our preliminary analysis shows that unlike rotation averaging, the Hessian of PGO is no longer well approximated by the corresponding graph Laplacian due to the coupling between rotation and translation terms. As a result, our current approach cannot be directly applied and additional techniques need to be considered. Lastly, in this work we have restricted our attention to standard least squares cost functions. Extending the theory and algorithm to problem formulations involving robust cost functions is another interesting direction for future research.

Appendix A Analysis of Riemannian Hessian of Rotation Averaging
In this appendix, we use d ∠ ≡ d to denote the geodesic distance on the rotation group, and use d chr to denote the chordal distance. Recall the definition of ϕ in Problem 1 as either the squared geodesic distance or the squared chordal distance, Using the notion of "reshaped distance" introduced in [52], we may express both cases as a function of the geodesic distance as follows, It can be shown that the scalar function ρ(·) is defined as, 2 − 2 cos(θ), for squared chordal distance (40b). (42b) The first case (42a) is readily verified. The second case (42b) makes uses of the relation between chordal and geodesic distances; see [53, Table 2]. To analyze the Hessian of rotation averaging (8), we start by considering the cost associated with a single relative rotation measurement, In addition to (43), we also consider its approximation defined on the Lie algebra, Note that (44) corresponds to a single term in the overall approximation defined in (11). Similar to (11), h ij depends on the current rotation estimates R, but we omit this from our notation for simplicity. Define the gradient and Hessian of h ij as follows, We prove the main theoretical results Theorem 1 and Corollary 1 for 3D rotation averaging. The case of d = 2 can be proved using the exact same arguments, and some steps would simplify due to the fact that 2D rotations commute. In the following, we first derive auxliary results that characterize g ij and H ij . Once we understand the properties of g ij and H ij , understanding the full rotation averaging problem becomes straightforward thanks to the additive structure in the cost function (8).

A.1 Auxiliary Results for 3D Rotation Averaging
Lemma 2. Consider a 3D rotation averaging problem. Let denote the angle-axis representation of the current rotation error. Then the gradient is given by, The Hessian is given by, where Proof. Introduce new rotation variables S i , S j ∈ SO(3) and consider the following function, Note that g ij and H ij correspond to the Riemannian gradient and Riemannian Hessian of (54) evaluated at (3) be the mapping such that, Then we have h ij (S i , S j ) = ϕ(F (S i , S j )). By chain rule, In (58), DF (S i , S j ) stands for the adjoint operator (transpose in matrix form) of the differential DF (S i , S j ). Using the standard basis for the Lie algebra so(3), Tron [52] showed that the Riemannian gradient of the geodesic distance is, c.f. [52,Equation (E.31)]. Substituting (59) into (58) and furthermore using the matrix form of the differential DF (S i , S j ) in the standard basis, we have For the Riemannian Hessian, differentiating (56) again yields, For the Hessian term in the middle of (61), we once again leverage existing results from [52, Proposition E.3.1]: where the inner matrix H ij is defined as, Using the expression for the differential of the logarithm map [52, Proposition E.2.1], the previous expression further simplifies to, To conclude, the Hessian is obtained by substituting above results into (61): The Hessian expression in Lemma 2 is complicated in general. However, we will show that as the angular error θ ij tends to zero, the Hessian H ij converges to a particular simple form.
Substituting into (51)- (53), Take the limit as θ tends to zero, Define the following matrix, From the definition of H ij in (50), Since lim θ→0 γ(θ) = lim θ→0 β(θ) = 0 and all matrices involved in (75) are bounded, we conclude that the last two terms in (75) vanish as θ ij converges to zero. For the first term in (75), notice that, As θ ij tends to zero, α(θ ij ) converges to 1 and the off-diagonal blocks in (76) tend to −I 3 . Hence the proof is completed.

A.2 Proof of Theorem 1
Proof. We will use Lemma 3 to prove the theorem for the case of squared geodesic cost. The case of squared chordal cost is analogous: instead of Lemma 3, we will use Lemma 4 and the remaining steps are the same.
Recall the approximation of the overall cost function defined in (11): Observe that the Hessian of h(v; R) is simply given by the sum of the Hessian matrices of h ij (v i , v j ), after "lifting" the latter to the dimension of the full optimization problem, i.e., In (86), W ij is formed by placing the p-by-p blocks of H ij defined in (50) in corresponding locations of the full matrix. For instance, H ij is the block of H ij that corresponds to vertex i. In the following, let θ ij (R) denote the residual of edge (i, j) ∈ E evaluated at R ∈ SO(d) n . From (86) and Lemma 3, we see that H(R) has the following limit point as all edge residuals tend to zero, Recall from Remark 2 that M 0 and ker(M ) = N where N is the vertical space defined in (12). In addition, recall the definition of H(R) in (14): Since P H is the (constant) orthogonal projection matrix onto the horizontal space H = N ⊥ , it holds that, Note that for singular symmetric matrices A and B, A ≈ δ B necessarily means that ker(A) = ker(B). Therefore, to prove the theorem, we must first show that ker(H(R)) = ker(M ) = N under our assumptions. Let λ 1 (A), λ 2 (A), . . . denote the eigenvalues of a symmetric matrix A sorted in increasing order. By construction, N is always contained in ker(H(R)), and thus dim(N ) = p eigenvalues of H(R) are always zero. Next, we will show that if all measurement residuals are sufficiently small, then the remaining pn − p eigenvalues of H(R) will be strictly positive. Define E(R) H(R) − M . Let x be any unit vector such that x ⊥ N . Note that, Since M = L(G; κ) ⊗ I p , it holds that λ p+1 (M ) = λ 2 (L(G; κ)). The latter is known as the algebraic connectivity which is always positive for a connected graph G. Thus λ p+1 (M ) > 0 and by (89), we also have lim θij (R)→0 E(R) = 0. Consequently, when all θ ij (R) are sufficiently small, the right-hand side of (90) is strictly positive, i.e., there existsθ 1 > 0 such that if θ ij (R) ≤θ 1 for all (i, j) ∈ E, we have, Under (91), the desired approximation H(R) ≈ δ M is equivalent to, where M Using the convergence result (89) and the eigenvalue perturbation bounds in [54,Corollary 6.3.8], we conclude that there existsθ 0 ∈ (0,θ 1 ] such that if R ∈ SO(d) n satisfies then (93) holds, i.e., we have the desired approximation, To conclude the proof, we need to show that there existθ, r > 0 such that condition (94) holds for all R ∈ B r (R ). Let us first consider residuals at the global minimizer R . Using assumption (15) and the cost function, we obtain the following simple bound: It can be verified that if,θ then (96) yields θ ij (R ) ≤θ 0 /2 for all edges (i, j) ∈ E. Finally, let us select r ∈ (0,θ 0 /4).
Using the triangle inequality and the fact that the geodesic distance d ∠ (·, ·) is bi-invariant, we can show that for any R ∈ B r (R ), In summary, we have shown that ifθ satisfies (97) and furthermore 0 < r <θ 0 /4, then the desired approximation M ≈ δ H(R) holds for all R ∈ B r (R ). This concludes the proof.

A.3 Proof of Corollary 1
Proof. To simplify notation, we use L to denote L(G; w). By (16), it holds that, Note that the eigenvalues of L⊗I p are given by the eigenvalues of L, repeated p times. Therefore, the desired result follows by noting that, Above, L cc is the block of the full Laplacian L that corresponds to the separators, denoted as L cc ≡ L(G) cc . Note that L cc can be decomposed as the sum, Intuitively, the first term in (107) accounts for inter-robot edges E c , and the second group of terms accounts for robots' local edges E α ; see Fig. 1a. Substitute (107) into (106), Sc(L(G α ), F α ). (108)

B.2 Proof of Theorem 2
Proof. Let us simplify the notations in the Laplacian system (26) by considering interior nodes from all robots as a single block: where By applying the Schur complement to (109), we obtain the following factorization for the input Laplacian system LX = B, where S = Sc(L, F) is the Schur complement that appears in (27). It can be verified that Algorithm 3 returns a solution to the following system, Recall from Lemma 1 that, Meanwhile, by construction, S is given by, where S α ≈ S α for all α ∈ [m]. Since spectral approximation is preserved under addition, it holds that S ≈ S. Furthermore, by comparing L defined in (114) and L defined in (115) and using [21,Fact 3.2], we conclude that L ≈ L and thus (31) is true. Lastly, (32) follows from Lemma 5.
Note that under the normal coordinates, we have x = 0, Exp x ( x) = x, and d(x, x ) = x . In addition, let us define, to be the coordinate expressions of vector fields, gradient, Hessian, and the retraction, respectively. Finally, let M denote the coordinate expression of the linear map M used in Algorithm 6: Let us express each iteration of Algorithm 6 in the chart, Using the triangle inequality, we can bound the distance between x k+1 and x , In the following, we will derive upper bounds for A and B as a function of x k − x = d(x k , x ).
Bounding A: Note that at x , we have D ϕ(x ) = I. Since ϕ is smooth, there exists r 1 > 0 such that for all It follows from (A2) and (126) that, Furthermore, Assumption (A1) implies that the gradient is Lipschitz continuous. Using the fact that g( x ) = 0 (since x is a critical point), we have the following upper bound for the norm of the Newton step, where L > 0 is a fixed constant. Using the local rigidity property of the retraction (e.g., see [23, Definition 4.1.1]), we have that, for all x in a neighborhood of x and all η sufficiently small; see also the discussions in [23, p. 115]. It follows from (132) and (131) that there exists r 2 , C 2 > 0 such that for all x k ∈ B r2 (x ).
Bounding B: To begin, we derive the following upper bound for B using triangle inequality, To bound B 1 , it follows from (130) that, Furthermore, in the proof of [23, Theorem 6.3.2], it is shown that there exist r 3 , C 3 > 0 such that, for x k ∈ B r3 (x ). Combining this result with (135), it holds that, It remains to establish an upper bound for the matrix that appears in B 2 : In the following, we first bound the norm of I − M (x k ) −1 Hess f (x k ). For any η ∈ T x k M, let us consider the following quantity, where ε H k = ε, Hess f (x k )[ε] denotes the norm induced by the Riemannian Hessian. Since M (x k ) ≈ Hess f (x k ) by Assumption (A3), we can use Lemma 5 to obtain an upper bound of (139), In addition, using Assumption (A1), it holds that, Combining (139)-(141) yields, Combining (142) and (129) in (138), we conclude that, Finishing the proof: We conclude the proof by combining the upper bounds (133), (137), and (143) in (128), and using the fact that The linear convergence factor in (120) is obtained by noting that the second term on the right-hand side vanishes at a quadratic rate.

C.2 Proof of Theorem 3
Proof. We will use the general linear convergence result established in Theorem 5. However, in order to properly account for the gauge symmetry in rotation averaging, we need to invoke Theorem 5 on the quotient manifold that underlies our optimization problem. In the following, we break the proof into three main parts (highlighted in bold). The proof makes heavy use of results regarding Riemannian quotient manifolds. The reader is referred to [24,Chapter 9] for a comprehensive review.
Rotation Averaging and Optimization on Quotient Manifolds. Following standard references, we denote a Riemannian quotient manifold as M = M/ ∼. For rotation averaging (Problem 1), the total space is given by M = SO(d) n . Let R = {R 1 , . . . , R n } and R = {R 1 , . . . , R n } be two points on M. We say that R and R are equivalent if they are related via a left group action: The equivalence class represented by R is defined as, Note that the cost function of Problem 1 is invariant within an equivalence class. In the following, we use T R M to denote the usual tangent space at R ∈ M, and T [R] M to denote the corresponding tangent space on the quotient manifold. Given a retraction Retr on M (see [24,Chapter 9.6]), we can execute the Riemannian Newton's method on M: Above, ξ k ∈ T [R k ] M is the solution of the following linear equation, where grad f ([R k ]) and Hess f ([R k ]) denote the Riemannian gradient and Hessian on the quotient manifold, respectively.
Expressing the iterates of Algorithm 4 on the quotient manifold. Recall that Algorithm 4 generates a sequence of iterates {R k } on the total space M. To prove convergence, we need to analyze the corresponding sequence of equivalence classes {[R k ]} on the quotient manifold M. Once we understand how [R k ] evolves, we can prove the desired result by invoking Theorem 5 on the quotient manifold M. To begin with, we write the update step in Algorithm 4 in the following general form: where v k is the vector corresponding to the matrix V k in line 10; see (19) for how v k and V k are related. Together, Retr R k (v k ) is the concise notation for the update steps in line 11-14, where we update each rotation R k i to Exp(v k i )R k i . Our notation Retr serves to emphasize that the retraction is performed on the total space M. Recall from (18)- (20) and Theorem 2 that v k is a solution to the linear system, where L ≈ L ≡ L(G; w) is defined in (31) and g(R k ) ∈ T R M is the Riemannian gradient in the total space defined in (13). In this proof, we will use the notations N R and H R to denote the vertical and horizontal spaces at R ∈ M. 9 By assumption, v k is the unique solution to (150) that satisfies v k ⊥ N R k , i.e., v k ∈ H R k . This implies that there is a unique tangent vector η k ∈ T [R k ] M on the tangent space of the quotient manifold such that v k is the horizontal lift [24, Definition 9.25] of η k at R k : v k = lift R k (η k ).
At any R ∈ M, define M (R) : H R → H R to be the linear map, where v ∈ H R is any vector from the horizontal space. Since ker By [24,Proposition 9.39], at iteration k, the Riemannian gradients on the total space and quotient space are related via, Combining (151)-(155), we see that (150) is equivalent to, Applying lift −1 R k to both sides of (156), Invoking Theorem 5 on the quotient manifold. To finish the proof, we will invoke Theorem 5 to show that the sequence of iterates [R k ] generated by (157)-(158) converges linearly to [R ]. This amounts to verifying that each condition in Assumption 1 holds on the quotient manifold M. To start, note that there exists r > 0 such that for any R ∈ SO(d) n , the condition d([R], [R ]) < r implies that R ∈ B r (R ) for some global minimizer R in the total space, where B r (R ) is the neighborhood within which Theorem 1 and Corollary 1 hold. Verifying (A1). We need to derive lower and upper bounds for the Hessian of the quotient optimization problem Hess f ([R]). For any η ∈ T [R] M, let v = lift R (η). By [24,Proposition 9.45], we have, where H(R) is defined in (14). Since v ∈ H R ≡ H belongs to the horizontal space, we conclude using Corollary 1 that v, H(R)v ≥ µ H v 2 where µ H is the constant defined in Corollary 1. Furthermore, since the quotient manifold M inherits the Riemannian metric from the total space M, we have η = v . We thus conclude that, The last inequality holds because v ⊥ N . Thus, we conclude that, Verifying (A3). Lastly, we need to show that the linear map M ([R]) is a spectral approximation of the Riemannian Hessian Hess f ([R]) on the quotient manifold. From Theorem 1, it holds that, In addition, from Theorem 2, we have Composing the two approximations yields, Note that the above result directly implies the following approximation relation on the quotient manifold, To see this, note that for any η ∈ T Proof. We note that the proof of a similar result can be found at [22,Claim 2.4]. In the following, we provide the proof for the case where L and L are singular. The non-singular case can be proved in the same way by replacing matrix psuedoinverse with the inverse. Observe that where X [:,i] denotes the i-th column of X . Therefore, we can first obtain an upper bound for the squared norm on a single column. To simplify notation, let x be a column of X , and let x and b be the corresponding columns of X and B. Let us expand the squared norm, Note that since L ≈ L, we have ker(L) = ker( L). In addition, any x where L x = b can be written as x = L † b + x ⊥ for some x ⊥ ∈ ker(L). Now, let us consider the middle term in (173), x L x = x L L † b = x L L † Lx = (L 1/2 x ) (L 1/2 L † L 1/2 )(L 1/2 x ).