Rapid Grassmannian Averaging with Chebyshev Polynomials

## Abstract

Abstract We propose new algorithms to efficiently average a collection of points on a Grassmannian manifold in both the centralized and decentralized settings. Grassmannian points are used ubiquitously in machine learning, computer vision, and signal processing to represent data through (often low-dimensional) subspaces. While averaging these points is crucial to many tasks (especially in the decentralized setting), existing methods unfortunately remain computationally expensive due to the non-Euclidean geometry of the manifold. Our proposed algorithms, Rapid Grassmannian Averaging (RGrAv) and Decentralized Rapid Grassmannian Averaging (DRGrAv), overcome this challenge by leveraging the spectral structure of the problem to rapidly compute an average using only small matrix multiplications and QR factorizations. We provide a theoretical guarantee of optimality and present numerical experiments which demonstrate that our algorithms outperform state-of-the-art methods in providing high accuracy solutions in minimal time. Additional experiments showcase the versatility of our algorithms to tasks such as K 𝐾 K italic_K -means clustering on video motion data, establishing RGrAv and DRGrAv as powerful tools for generic Grassmannian averaging. 1 1 1 All code available at https://github.com/brightonanc/rgrav

## 1 Introduction

Grassmannian manifolds, which represent sets of $K$-dimensional linear subspaces of $N$-dimensional spaces (Edelman et al., ), have been used extensively in machine learning (Huang et al., ; Zhang et al., ; Slama et al., ), computer vision (Harandi et al., ; Lui & Beveridge, ; Turaga et al., ), and signal processing (Gallivan et al., ; Mondal et al., ; Xu & Hassibi, ). Applications include Principal Component Analysis (PCA) (Jolliffe & Cadima, ), low-rank matrix completion (Keshavan et al., ), multi-task feature learning (Mishra et al., ), clustering (Gruber & Theis, ), array processing (Love et al., ), and distance metric learning (Meyer et al., ).

An essential primitive operation is finding an average of a collection of points on the manifold. There are several distinct yet reasonable definitions for an average of points on a Grassmannian (Marrinan et al., ). Arguably the most natural analog of the Euclidean mean for points on a Riemannian manifold (such as a Grassmannian) is the Fréchet (or Karcher) mean, defined as the point which minimizes the sum of squared distances to all sample points (Fréchet, ). Unfortunately, the Fréchet mean rarely admits a closed form solution, instead necessitating approximation via iterative algorithms (Jeuris et al., ). Such algorithms are often computationally expensive, scale poorly with dimension, and are not easily decentralized.

The induced arithmetic mean (IAM) is an alternative manifold average computed by first determining the Euclidean mean of the manifold sample points once embedded “naturally” in some Euclidean space and subsequently projecting this Euclidean mean “naturally” back onto the manifold. For Grassmannian manifolds, the standard embedding is the set of projection matrices and the standard projection operation is simply the closest matrix by Frobenius distance (Sarlette & Sepulchre, ). This manifold average may be computed much more efficiently in practice and lends itself well to decentralization as the Euclidean mean may be computed by average consensus (Nedic & Ozdaglar, ) (subsequent projections may be applied locally thereafter).

As the dimensionality of data grows, it becomes increasingly important to consider decentralized algorithms (Nedić et al., ) as data might be spread across many machines and only be accessible for processing via distributed algorithms (Beltrán et al., ). While a central server is sometimes employed in this regime, it is similarly common for the use of such a server to be infeasible or simply inefficient when compared to fully decentralized approaches (Sun et al., ; Feller et al., ).

We propose a novel method to efficiently compute the IAM of a collection of points on a Grassmannian manifold. Our method is highly amenable to decentralization, meaning it can be readily deployed to multi-agent systems or used in data centers operating on big data. Our algorithms operate similarly to the famous power method, with the distinction that Chebyshev polynomials are employed to leverage a “dual-banded” property of the problem in order to achieve never-before-seen efficiency in computation and communication. We demonstrate merit through a theoretical guarantee on the optimality of our approach among a class of algorithms, synthetic numerical experiments comparing our algorithms against state-of-the-art, and experiments on real-world problems showcasing the versatility of our algorithms.

## 2 Related Work

The problem of computing an appropriate average on specific manifolds has been investigated for many different manifolds, e.g., $S^{N}$, $\text{SO}\left(N\right)$, $\text{St}\left(N,K\right)$, even $\text{Gr}\left(N,K\right)$ (Downs, ; Buss & Fillmore, ; Galperin, ; Hueper & Manton, ; Absil et al., ; Moakher, ; Fiori et al., ; Yun, ). Focus is often given to the Fréchet mean (Chakraborty et al., ; Cheng et al., ; Le, ), however alternatives are becoming increasingly more popular (Fletcher et al., ; ; Arnaudon et al., ; Marrinan et al., ; Chakraborty & Vemuri, ; Lee & Jung, ). Similarly, the problem of consensus on a manifold in a multi-agent setting has been explored in works such as Sepulchre (); Tron et al. ().

There have been several algorithms proposed for decentralized optimization on manifolds such as Grassmannians. Sarlette & Sepulchre () proposes a decentralized gradient-based algorithm to solve the problem of computing the IAM for connected compact homogeneous manifolds, e.g. $\text{SO}\left(N\right)$ and $\text{Gr}\left(N,K\right)$. Deng & Hu () proposes two decentralized gradient-based algorithms for general optimization problems on Riemannian manifolds. Mishra et al. () proposes a decentralized gradient-based gossip algorithm for general optimization problems on a Grassmannian manifold. Similar works include Chen et al. (; ); Zhang & Sun ().

A problem which is closely related to Grassmannian averaging is that of PCA. Ye & Zhang () proposes the DeEPCA algorithm to solve the decentralized PCA problem. While computing a Grassmannian average is not the intended application of DeEPCA, it may be adapted to this task fairly naturally. Gang et al. (); Gang & Bajwa (); Froelicher et al. () similarly propose distributed algorithms for PCA; for a more comprehensive review of this field, see Wu et al. ().

## 3 Background

### 3.1 Averaging Subspaces

Given a collection of $M$ subspaces, our goal is to determine the average subspace as efficiently as possible. We choose to use the standard IAM definition of “average” (Sarlette & Sepulchre, ) as it leads to what we believe is the most efficient algorithm. Formally, let $\text{St}\left(N,K\right):=\left\{\boldsymbol{U}\in\mathbb{R}^{N\times K}~{}|~
{}\boldsymbol{U}^{\mathsf{T}}\boldsymbol{U}=\mathbf{I}_{K}\right\}$ be the set of $N\times K$ Stiefel matrices where $N\geq K$ and let $\text{Gr}\left(N,K\right):=\left\{\left[\boldsymbol{U}\right]~{}|~{}
\boldsymbol{U}\in\text{St}\left(N,K\right)\right\}$ (where the equivalence is defined as $\left[\boldsymbol{U}\right]:=\left\{\boldsymbol{U}\boldsymbol{Q}~{}|~{}
\boldsymbol{Q}\in\text{St}\left(K,K\right)\right\}$) be the Grassmannian representing the set of all $K$-dimensional subspaces of $\mathbb{R}^{N}$. The average of our collection $\left\{\left[\boldsymbol{U}_{m}\right]\right\}_{m=1}^{M}$ is then denoted $\left[\bar{\boldsymbol{U}}\right]$ and defined by the following optimization problem

$$ $\left[\bar{\boldsymbol{U}}\right]:=\operatorname*{argmin}_{\left[\boldsymbol{U }\right]\in\text{Gr}\left(N,K\right)}\left\lVert\frac{1}{M}\sum_{m=1}^{M} \boldsymbol{U}_{m}\boldsymbol{U}_{m}^{\mathsf{T}}-\boldsymbol{U}\boldsymbol{U} ^{\mathsf{T}}\right\rVert_{\mathrm{F}}^{2}$ (1) $$

Equation 1 may be manipulated algebraically to be interpreted equivalently in terms of the eigenvectors of $\bar{\boldsymbol{P}}:=\frac{1}{M}\sum_{m=1}^{M}\boldsymbol{U}_{m}\boldsymbol{U
}_{m}^{\mathsf{T}}$. Let

$$ $\bar{\boldsymbol{P}}=\tilde{\boldsymbol{V}}\tilde{\boldsymbol{\Lambda}}\tilde{ \boldsymbol{V}}^{\mathsf{T}}=\begin{bmatrix}\boldsymbol{V}&\boldsymbol{V}_{ \perp}\end{bmatrix}\begin{bmatrix}\boldsymbol{\Lambda}&\boldsymbol{0}\\ \boldsymbol{0}&\boldsymbol{\Lambda}_{\perp}\end{bmatrix}\begin{bmatrix} \boldsymbol{V}^{\mathsf{T}}\\ \boldsymbol{V}_{\perp}^{\mathsf{T}}\end{bmatrix}$ $$

denote an eigendecomposition where $\boldsymbol{V}\in\text{St}\left(N,K\right)$ and the entries of $\tilde{\boldsymbol{\Lambda}}$ are non-increasing. Assuming $\lambda_{K}>\lambda_{K+1}$ (where $\lambda_{k}$ denotes the $k$th largest eigenvalue of $\bar{\boldsymbol{P}}$), it can be shown that the solution to eq. 1 is precisely $\left[\bar{\boldsymbol{U}}\right]=\left[\boldsymbol{V}\right]$. Consequently, determining the span of the leading $K$ eigenvectors of $\bar{\boldsymbol{P}}$ is tantamount to solving eq. 1, which is the perspective we will later use to motivate our algorithms.

As we continue to discuss this problem, it is informative to keep in mind the following properties. The eigenvalues of $\bar{\boldsymbol{P}}$ are conveniently bounded by $\boldsymbol{0}\preceq\tilde{\boldsymbol{\Lambda}}\preceq\mathbf{I}_{N}$ and satisfy $\operatorname{tr}\left(\tilde{\boldsymbol{\Lambda}}\right)=K$, which may be determined by inspection. As a result, $\lambda_{K},\lambda_{K+1}$ are bounded as $\frac{1}{N-K+1}\leq\lambda_{K}\leq 1$ and $0\leq\lambda_{K+1}\leq\frac{K}{K+1}$. For convenience, we occasionally abuse notation to let $\left[\boldsymbol{X}\right]$ denote the Grassmannian equivalence class for the span of the columns of arbitrary (not necessarily Stiefel) matrix $\boldsymbol{X}\in\mathbb{R}^{N\times K}$.

### 3.2 Averaging Subspaces in a Decentralized Network

Decentralized or distributed optimization problems arise in numerous real-world scenarios where centralized approaches are impractical or undesirable. These problems are characterized by using information spread across multiple agents or nodes in a network. Motivating reasons include privacy, communication constraints, data storage limitations, scalability, etc.

In the context of this paper, consider the setting where there are $M$ agents, each holding a subspace $\left[\boldsymbol{U}_{m}\right]$, connected by a some undirected communication graph $\mathcal{G}$. We then want each agent to learn the solution $\left[\bar{\boldsymbol{U}}\right]$ to eq. 1 under the restriction that each agent only communicate with their neighbors in $\mathcal{G}$.

Average consensus (AC) is a useful primitive in decentralized optimization to quickly approximate the average of real numbers in a decentralized manner (Nedic & Ozdaglar, ). Unfortunately, the non-convex manifold structure of $\text{Gr}\left(N,K\right)$ precludes us from efficiently applying AC directly to solve eq. 1. While we could have each agent compute the matrices $\boldsymbol{U}_{m}\boldsymbol{U}_{m}^{\mathsf{T}}$ and then use AC to approximate $\bar{\boldsymbol{P}}$, this would incur a communication cost of $\mathcal{O}\left(N^{2}\right)$ (the size of $\bar{\boldsymbol{P}}$) which may be much larger than $\mathcal{O}\left(NK\right)$. We could instead use AC to average the matrices $\boldsymbol{U}_{m}$ with preferable communication cost $\mathcal{O}\left(NK\right)$, however the arbitrary choice of representative Stiefel matrix $\boldsymbol{U}_{m}$ from the Grassmannian equivalence class $\left[\boldsymbol{U}_{m}\right]$ makes this approach ill-posed.

In order to achieve the $\mathcal{O}\left(NK\right)$ communication cost without being ill-posed, one can have all agents compute $\boldsymbol{U}_{m}\boldsymbol{U}_{m}^{\mathsf{T}}\boldsymbol{X}$ and then use AC to approximate $\bar{\boldsymbol{P}}\boldsymbol{X}$, where $\boldsymbol{X}\in\mathbb{R}^{N\times K}$ is some matrix agreed upon by all agents a priori. In practice, the requirement that all agents agree upon $\boldsymbol{X}$ might be overly strict; in many scenarios, it often suffices to have each agent $m$ have instance $\boldsymbol{X}_{m}$ which are all approximately equal (i.e. there exists some $\boldsymbol{X}$ for which $\boldsymbol{X}_{m}\approx\boldsymbol{X}$ for all $m\in\left[M\right]$).

After one iteration of an algorithm, each agent will have some local approximation of the quantity $\bar{\boldsymbol{P}}\boldsymbol{X}$. If the matrix $\boldsymbol{X}$ is retained in memory for each agent, then AC may be applied to a linear combination of the $\bar{\boldsymbol{P}}\boldsymbol{X}$ approximations and $\boldsymbol{X}$ to approximate some quantity $\bar{\boldsymbol{P}}\left(a\bar{\boldsymbol{P}}\boldsymbol{X}+b\boldsymbol{X}
\right)=\left(a\bar{\boldsymbol{P}}^{2}+b\bar{\boldsymbol{P}}\right)
\boldsymbol{X}$. Applying this logic recursively reveals that such an algorithm can approximate $f_{t}\left(\bar{\boldsymbol{P}}\right)\boldsymbol{X}$ after $t$ iterations, where $f_{t}\in\mathcal{P}_{t}$ is a $t$th order polynomial; in Section 4 we will consider more thoroughly these polynomials and how they can translate to desirable algorithms.

Gradient tracking is a famous technique in decentralized optimization that improves upon the convergence rate of AC-based methods (Shi et al., ; Xu et al., ; Qu & Li, ; Deng & Hu, ; Ye & Zhang, ). In essence, it sets up a recursion using standard AC whose fixed point satisfies both a consensus condition (meaning all agents agree) and a stationarity condition (meaning the solution is locally optimal). We will later employ a form of gradient tracking over quantities of the form $f_{t}\left(\bar{\boldsymbol{P}}\right)\boldsymbol{X}$ for our decentralized algorithms.

### 3.3 The Power Method

The power method is a classical algorithm to estimate the leading eigenspace of a positive semidefinite matrix $\boldsymbol{A}\in\mathbb{R}^{N\times N}$. A single power iteration applies the matrix $\boldsymbol{A}$ and orthonormalizes (for numerical stability) the result, e.g.

$$ $\boldsymbol{U}^{(t)}=\text{QR}\left(\boldsymbol{A}\boldsymbol{U}^{(t-1)}\right)$ $$

The power method loop may be unrolled to reveal the following form

$$ $\boldsymbol{U}^{(t)}=\text{QR}\left(\underbrace{\boldsymbol{A}\boldsymbol{A} \cdots\boldsymbol{A}\boldsymbol{A}}_{\text{t times}}\boldsymbol{U}^{(0)}\right )=\text{QR}\left(\boldsymbol{A}^{t}\boldsymbol{U}^{(0)}\right)$ $$

Let $\boldsymbol{V}_{\star}\in\text{St}\left(N,K\right)$ be a basis for the leading $K$ eigenvectors of $\boldsymbol{A}$.
For random initialization $\boldsymbol{U}^{(0)}\in\mathbb{R}^{N\times K}$, the span of $\boldsymbol{U}^{(t)}$ converges to the span of the $\boldsymbol{V}_{\star}$ provided $\operatorname{rank}\left(\boldsymbol{V}_{\star}^{\mathsf{T}}\boldsymbol{U}^{(0
)}\right)=K$ (Golub & Van Loan, , Chapter 8.2).

At iteration $t$, the power method effectively applies the function $f_{t}(\lambda)=\lambda^{t}$ to each eigenvalue of $\boldsymbol{A}$. Consider for example the case where $\lambda_{K}\left(\boldsymbol{A}\right)=1$ and $\lambda_{K}\left(\boldsymbol{A}\right)-\lambda_{K+1}\left(\boldsymbol{A}\right
)>0$. As $t$ increases, the ratio between trailing and leading eigenvalues of $\boldsymbol{A}^{t}$ shrinks exponentially; formally, for $1\leq k\leq K$ and $K+1\leq\ell\leq N$ we have the following

$$ $\frac{\lambda_{\ell}\left(\boldsymbol{A}^{t}\right)}{\lambda_{k}\left( \boldsymbol{A}^{t}\right)}\leq\frac{\lambda_{K+1}\left(\boldsymbol{A}^{t} \right)}{\lambda_{K}\left(\boldsymbol{A}^{t}\right)}=\left(\frac{\lambda_{K+1} \left(\boldsymbol{A}\right)}{\lambda_{K}\left(\boldsymbol{A}\right)}\right)^{t }<1$ $$

While the power method works well, large values of $\lambda_{K+1}$ can slow convergence. We will later show how polynomials other than $\lambda^{t}$ can overcome this shortcoming and use them in our algorithms.

## 4 Methods

### 4.1 Motivation

The goal of the RGrAv algorithms is to solve eq. 1 as efficiently as possible. Recall from Section 3 that determining the span of the $K$ leading eigenvectors of $\bar{\boldsymbol{P}}$ is tantamount to solving eq. 1. Motivated by the decentralized setting and the power method, we restrict our consideration to iterative algorithms which after $t$ iterations compute $f_{t}\left(\bar{\boldsymbol{P}}\right)\bar{\boldsymbol{U}}^{(0)}$ for some $t$th-order polynomial $f_{t}$ and initial estimate $\bar{\boldsymbol{U}}^{(0)}$ (see Section 3.2).

We consider the case where the spectrum of $\bar{\boldsymbol{P}}$ is dual-banded, i.e. there exists some $0<\alpha<\beta<1$ such that $\boldsymbol{\Lambda}_{\perp}\preceq\alpha\mathbf{I}_{N-K}$, $\beta\mathbf{I}_{K}\preceq\boldsymbol{\Lambda}$, and $\beta-\alpha\gg 0$ (e.g. $\beta-\alpha=\frac{1}{3}$). This situation can arise, for instance, when points are normally distributed on the manifold (see Section 5.1). For simplicity, we refer to the intervals $\left[0,\alpha\right]$ and $\left[\beta,1\right]$ as the “stop-band” and “pass-band”, respectively. Similar to the power method, we would like our polynomial $f_{t}$ to decrease the ratio between eigenvalues in the stop-band relative to eigenvalues in the pass-band. However, unlike the power method, knowledge of this dual-banded structure (even heuristically) allows us to choose polynomials which optimize the worst case value of this ratio criteria. Leveraging this spectral structure is how we will choose our optimal polynomials $f_{t}^{\star}$ (see Theorem 1).

There is, however, an important consideration for high-dimensional data. In the case where $N\gg MK$ we are guaranteed that the nullspace of $\bar{\boldsymbol{P}}$ is (at least $N-MK$) high-dimensional by rank subadditivity, meaning there will be a cluster of eigenvalues at 0. For this reason, we will constrain our polynomials $f_{t}$ to always satisfy $f_{t}(0)=0$. Since the ratio criteria of Theorem 1 is invariant to scaling of $f_{t}$, we provide one final constraint of $f_{t}(1)=1$ simply for uniqueness of solution and numerical stability.

###### Theorem 1 .

For $t\geq 1$, the minimization problem

$$ $\operatorname*{minimize}_{f_{t}\in\mathcal{P}_{t}^{\prime}}\frac{\max_{\lambda \in\left[0,\alpha\right]}\left\lvert f_{t}\left(\lambda\right)\right\rvert}{ \min_{\lambda\in\left[\beta,1\right]}\left\lvert f_{t}\left(\lambda\right) \right\rvert},$ (2) $$

where $\mathcal{P}_{t}^{\prime}$ is the set of $t$th order polynomials such that $f_{t}(0)=0$ and $f_{t}(1)=1$, is solved by

$$ $f_{t}^{\star}\left(\lambda\right)=\prod_{s=0}^{t-1}\frac{\lambda-r_{s,t}}{1-r_ {s,t}},\qquad r_{s,t}:=\alpha\frac{\cos\left(\frac{\pi\left(s+1/2\right)}{t} \right)+\cos\left(\frac{\pi}{2t}\right)}{1+\cos\left(\frac{\pi}{2t}\right)}$ $$

which is a modification of a Chebyshev polynomial of the first kind.(^2^22A proof may be found in Section B.1)

Figure: Figure 1: A visual comparison between the power method and our method. Each method’s corresponding $t$th-order polynomial is applied to the eigenvalues $\lambda$ in the domain $\left[0,1\right]$. The Chebyshev recursion with threshold parameter $\alpha=0.5$ results in the polynomial oscillations being reduced and flattened in the range $[0,\alpha]$.
Refer to caption: extracted/5920205/eig_polynomial.png

###### Theorem 1 .

### 4.2 RGrAv Algorithms

The polynomial $f_{t}^{\star}$ may be exactly implemented by iteratively multiplying each factor in the product given in Theorem 1; the “finite” variants of the RGrAv algorithms (Algorithms 5 and 6) do precisely this. While this approach may be acceptable when the number of iterations $t$ is known in advance, intermediate solutions can be quite sub-optimal.

Ideally, one would be able to describe $f_{t}^{\star}$ in terms of only a constant number of previous $f_{s}^{\star}$, e.g. $f_{t-1}^{\star},f_{t-2}^{\star}$. This would yield an efficient algorithm with optimal intermediate solutions whose memory/compute costs do not grow as $t\to\infty$. Unfortunately, this is not the case; fortunately, $f_{t}^{\star}$ is well-approximated in terms of $f_{t-1}^{\star},f_{t-2}^{\star}$. For $t\geq 2$, coefficients $a_{t},b_{t},c_{t}$ are chosen such that

$$ $\tilde{f}_{t}^{\star}\left(\lambda\right)=a_{t}\left(\left(\lambda+b_{t}\right )f_{t-1}\left(\lambda\right)+c_{t}f_{t-2}\left(\lambda\right)\right)$ $$

matches $f_{t}^{\star}\left(\lambda\right)$ in its leading three terms, i.e. $f_{t}^{\star}\left(\lambda\right)-\tilde{f}_{t}^{\star}\left(\lambda\right)\in
\mathcal{P}_{t-3}$ (see Algorithm 7); the “asymptotic” variants of the RGrAv algorithms (Algorithms 1 and 2) use this $\tilde{f}_{t}^{\star}$.

As discussed in Section 3.3, orthonormalization must occur periodically for numerical stability. To minimize the frequency of the orthonormalization schedule, our algorithms effectively cache the operation of the most recent exact orthonormalization to a matrix $\boldsymbol{S}$ and then efficiently approximate the orthonormalization procedure for iterations between the schedule by left-application of $\boldsymbol{S}$. We choose the QR factorization for our orthonormalization method, however alternative methods would be acceptable.

In the centralized setting, the RGrAv algorithms need not worry about inaccuracy from average consensus and so the order of operations focuses on minimizing the number of computations performed. In the decentralized setting, the order of operations focuses on minimizing the error in the gradient tracking procedure. Additionally, in the decentralized setting one must take care to use a “stable” orthonormalization method which will not change drastically with small perturbations in the input. We omit the nuances of numerical linear algebra that lead to this problem (see Golub & Van Loan (, Chapter 5) for more information) and simply present Algorithm 3, which is a stable wrapper for any implementation of a possibly unstable QR factorization.

Figure: Algorithm 1 Asymptotic RGrAv (Rapid Grassmannian Averaging)

Figure: Algorithm 2 Asymptotic DRGrAv (Decentralized Rapid Grassmannian Averaging)

Figure: Algorithm 3 StableQR

## 5 Experiments

### 5.1 Decentralized Grassmannian Averaging

In these experiments, we consider the problem where a network of $M$ connected agents each has a local instance of a Grassmannian basis $\boldsymbol{U}_{m}\in\text{St}\left(N,K\right)$ and we would like for all agents to learn an average Grassmannian basis of all $\boldsymbol{U}_{m}$ in a strictly decentralized manner (i.e. there is no central server, all communication is neighbor-to-neighbor). Our experiments had parameters $M=64$, $N=150$, $K=30$. To demonstrate the practicality of DRGrAv in both well-connected and sparse communication graphs, we performed experiments for two communication graphs: the hypercube graph and the cycle graph.

We compared DRGrAv to several alternative methods for Grassmannian averaging. Given below are the algorithms, their sources, and considerations for their tuning such that the comparison would be fair.

DRGrAv (this paper): Contrary to the following algorithms for which hyperparameters were chosen over large ranges to be empirically optimal, the hyperparameter $\alpha$ is chosen here heuristically as $0.15$. We also choose to use the approximate asymptotic variant of DRGrAv, and set the orthonormalization schedule to orthonormalize at every iteration (to match DeEPCA). These choices were made to demonstrate DRGrAv’s competitive performance even when sub-optimally tuned.

DeEPCA (Ye & Zhang, ): While this algorithm is not intended for decentralized Grassmannian averaging, we found that it could be easily adapted to this setting and gave competitive results: one simply substitutes $\boldsymbol{U}_{m}\boldsymbol{U}_{m}^{\mathsf{T}}$ for the paper’s $\boldsymbol{A}_{j}$. Also, for numerical stability, the paper’s QR + SignAdjust procedure is replaced with the StableQR procedure.

DPRGD/DPRGT (Deng & Hu, ): These algorithms are adapted to the decentralized Grassmannian averaging problem by using $-\frac{1}{2}\left\lVert\boldsymbol{U}_{m}^{\mathsf{T}}\bar{\boldsymbol{U}}_{m}
^{(t)}\right\rVert_{\mathrm{F}}^{2}$ for the paper’s $f_{i}$. These algorithms each have a single hyperparameter for step size $\alpha$, which was chosen for each algorithm (up to precision of 2 significant figures) by searching for the value $\alpha\in\left[10^{-4},10^{3}\right]$ which approximately minimized the MSE; the solutions were $\alpha$ on the order of $10^{0}$.

COM (Consensus Optimization on Manifolds) (Sarlette & Sepulchre, ): This algorithm is the discrete-time variant of the continuous-time dynamics presented in Equation 20 of Sarlette & Sepulchre (). There is a single hyperparameter for step size $\alpha$, which was chosen (up to precision of 2 significant figures) by searching for the value in $\alpha\in\left[10^{-4},10^{3}\right]$ which approximately minimized the MSE; the solutions were $\alpha$ on the order of $10^{-1}$.

Gossip (Mishra et al., ): This is Algorithm 1 of Mishra et al. () where $\frac{1}{4}\sum_{m=1}^{M}\sum_{n\in\mathcal{N}(m)}d^{2}\left(\left[\bar{U}_{m}
^{(t)}\right],\left[\bar{U}_{n}^{(t)}\right]\right)$ is used as the paper’s $g$ (where $\mathcal{N}(m)$ is the set of neighbors of $m$). This algorithm is a gossip algorithm, not an average-consensus-based algorithm. As a result, to keep the comparison fair we let each agent perform a gradient step in parallel for each round of consensus the other algorithms perform. In precise terms, during each round of consensus this algorithm will select edges from the graph uniformly at random until there no longer remain any 2 neighboring agents who both have not yet been selected; each of these agents then performs a gradient step and the process repeats. There are 2 hyperparameters (w.l.o.g. $\rho=1$) for step size $a$ and $b$, which were chosen (up to precision of 2 significant figures) by searching for values in $a\in\left[10^{-4},10^{3}\right],b\in\left[10^{-8},10^{0}\right]$ which approximately minimized the MSD; the solutions were $a$ on the order of $10^{0}$ and $b$ on the order of $10^{-4}$.

Figure: (a) (Hypercube) Mean Squared Error
Refer to caption: x1.png

**Table 1: Comparison of runtimes for various algorithms. The first five data columns display time (in milliseconds) until the MSE across agents goes below the given tolerance. COM and Gossip do not in general converge to any specific point, so their metric for tolerance is instead MSD. The minimal quantity in each column is bolded. The final column represents time per algorithm iteration. These decentralized algorithms are not truly run on separate devices, only simulated as such, so these runtimes should be interpreted broadly as general evidence that DRGrAv would perform well in a true decentralized setting.**
| Time (ms) until tolerance… | 1e-3 | 1e-6 | 1e-9 | 1e-12 | 1e-15 | Per Iter. |
| --- | --- | --- | --- | --- | --- | --- |
| DRGrAv (This Paper) | 35.4 | 47.4 | 47.4 | 56.7 | 66.1 | 11.8 |
| DeEPCA Ye & Zhang () | 35.8 | 45.7 | 53.0 | 61.1 | 80.6 | 9.13 |
| DPRGD Deng & Hu () | 1860 | 61200 | $>$100000 | $>$100000 | $>$100000 | 9.24 |
| DPRGT Deng & Hu () | 2270 | 2910 | 3470 | 4200 | 4780 | 14.5 |
| COM* Sarlette & Sepulchre () | 5050 | 7290 | 9260 | 11000 | 13200 | 16.3 |
| Gossip* Mishra et al. () | 1280 | 1730 | 2150 | 2590 | 3390 | 427 |

In order to have a fair comparison, all average-consensus-based algorithms mentioned above used the same consensus protocol. Both graphs used the optimal Laplacian-based communication matrix (i.e. $\boldsymbol{W}=\mathbf{I}-\frac{1}{7}\boldsymbol{L}$ for the hypercube graph, $\boldsymbol{W}\approx\mathbf{I}-\frac{1}{2}\boldsymbol{L}$ for the cycle graph, where $\boldsymbol{L}$ is the corresponding graph Laplacian matrix) for 10 rounds of communication in the hypercube graph case and 50 rounds of communication in the cycle graph case.

A single synthetic dataset $\left\{\boldsymbol{U}_{m}\right\}_{m=1}^{M}$ of “normally distributed points with standard deviation $\frac{\pi}{4}$” was used for all experiments. In precise terms, said dataset was generated by sampling a center point $\boldsymbol{U}_{\mathrm{C}}$ uniformly at random on $\text{Gr}\left(N,K\right)$ and then computing $\boldsymbol{U}_{m}:=\exp_{\boldsymbol{U}_{\mathrm{C}}}\left(\boldsymbol{T}_{m}\right)$ for all $m\in\left[M\right]$, where $\boldsymbol{T}_{m}:=\tilde{\boldsymbol{U}}_{m}\tilde{\boldsymbol{\Sigma}}_{m}
\tilde{\boldsymbol{V}}_{m}^{\mathsf{T}}$ is a random tangent vector at $\boldsymbol{U}_{\mathrm{C}}$ such that $\tilde{\boldsymbol{U}}_{m},\tilde{\boldsymbol{V}}_{m}$ are sampled uniformly at random from sets $\left\{\tilde{\boldsymbol{U}}~{}|~{}\tilde{\boldsymbol{U}}\in\text{St}\left(N,
K\right),\boldsymbol{U}_{\mathrm{C}}^{\mathsf{T}}\tilde{\boldsymbol{U}}=
\boldsymbol{0}\right\},\text{St}\left(K,K\right)$ respectively and $\tilde{\boldsymbol{\Sigma}}:=\operatorname{Diag}\left(\frac{\pi}{4}\boldsymbol
{z}\right)$ where $\boldsymbol{z}\sim\mathcal{N}\left(\boldsymbol{0}_{K},\mathbf{I}_{K}\right)$ is a vector draw of $K$ i.i.d. standard normal random variables.

The Mean Squared Error quantity at time $t$ was computed as $\frac{1}{M}\sum_{m=1}^{M}d^{2}\left(\left[\bar{\boldsymbol{U}}_{m}^{(t)}\right
],\left[\bar{\boldsymbol{U}}\right]\right)$ where $d$ is the extrinsic (or chordal) distance on the Grassmannian defined as $d\left(\left[\boldsymbol{U}_{1}\right],\left[\boldsymbol{U}_{2}\right]\right):
=2^{-1/2}\left\lVert\boldsymbol{U}_{1}\boldsymbol{U}_{1}^{\mathsf{T}}-
\boldsymbol{U}_{2}\boldsymbol{U}_{2}^{\mathsf{T}}\right\rVert_{\mathrm{F}}$ and $\bar{\boldsymbol{U}}$ is the true IAM average (see Section 3). Similarly, the Mean Squared Disagreement quantity represents the amount to which the agents’ estimates vary at time $t$ and was computed as $\frac{2}{M(M-1)}\sum_{m=1}^{M}\sum_{n=m+1}^{M}d^{2}\left(\left[\bar{
\boldsymbol{U}}_{m}^{(t)}\right],\left[\bar{\boldsymbol{U}}_{n}^{(t)}\right]\right)$.

The results of these experiments are shown in Figure 2 and Table 1. Six iterations were chosen for the hypercube graph case because at this point DRGrAv reaches floating point tolerance (FPtol), demonstrating that such precision is achievable by an algorithm in such time; Nine iterations were chosen for the cycle graph in order to demonstrate the effective consensus-permitted tolerance (ECPtol) of around $10^{-10}$. DRGrAv performs the best out of all algorithms, converging in the hypercube graph case to FPtol in only 6 iterations and converging in the cycle graph case to ECPtol in only 5 iterations. The adapted DeEPCA method performs most closely to DRGrAv, however still lags behind several orders of magnitude in MSE. In the cycle graph, DeEPCA manages to barely beat DRGrAv in at the end, however given both are more or less at ECPtol we do not think this provides strong evidence to prefer DeEPCA to DRGrAv. Since both COM and Gossip begin with $\bar{\boldsymbol{U}}_{m}^{(0)}\leftarrow\boldsymbol{U}_{m}$ instead of some pre-agreed upon starting $\boldsymbol{U}^{(0)}$, they are able to have superior performance to DRGrAv in the short term; however after only 2 iterations this short term behavior ends. DPRGD and DPRGT, being generic algorithms for use on any compact submanifold, do not leverage any of the structure specific to the Grassmannian problem and consequently are not empirically competitive to algorithms which do, e.g. DRGrAv. All algorithms presented have their specific ideal use cases, and we claim that the problem of decentralized Grassmannian averaging is the ideal use case of DRGrAv.

### 5.2 K-Means for Video Motion Clustering

We consider the application of the RGrAv algorithm to the problem of video motion analysis, extending the work of Marrinan et al. (). Their study applied centralized subspace averaging methods to multiple tasks on the DARPA Mind’s Eye video dataset. In our work, we focus specifically on the task of $K$-means clustering and compare against the best algorithm presented in their work.

The Mind’s Eye dataset consists of a set of “tracklets” — short grayscale videos sequences of moving objects, primarily people. Each tracklet consists of 48 frames of size $32\times 32$ pixels. To prepare these tracklets for subspace analysis, they are flattened into matrices $\boldsymbol{X}_{t}\in\mathbb{R}^{1024\times 48}$, where each column represents a vectorized frame. The subspaces $\boldsymbol{U}_{t}=\rm{span}(\boldsymbol{X}_{t})$ spanned by these columns are treated as points on the Grassmannian, effectively encoding the essential motion patterns in the video.

Each tracklet is annotated with a label describing the type of motion it contains, such as “walk” or “ride-bike.” These labels provide ground truth for evaluating the effectiveness of clustering algorithms; clusters are considered high-quality if most tracklets in the cluster share the same label. In Marrinan et al. (), the authors find that the choice of averaging algorithm does not significantly affect cluster quality as the number of clusters increases, whereas runtime can differ significantly. As a result, we compare to the fastest averaging algorithm from their work – the flag mean – and show that applying RGrAv can reduce runtime for the clustering task.

Figure: Figure 3: A comparison of runtime for K-means with various averaging algorithms and numbers of clusters K. The four colors represent the averaging algorithm as RGrAv (green), flag mean (orange), Fréchet mean (blue), and power method (pink). The four algorithms produce clusters with similar quality (excluded for brevity), but the RGrAv algorithm is significantly faster, showing $2\times$-$10\times$ speedup over the other averaging algorithms.
Refer to caption: extracted/5920205/runtime_comparison_power.png

The standard $K$-means algorithm can be extended to cluster points on the Grassmannian by defining two primitives: a distance metric and an averaging operation. The standard $K$-means algorithm with these operations is shown in Algorithm 4. The centers $\bar{\boldsymbol{U}}_{c}$ are initialized randomly. At each iteration, the points $\boldsymbol{U}_{t}$ in the dataset are each assigned to their closest mean using the metric to form clusters. The means are then updated to the average of their respective clusters, and the steps are repeated until the means converge.

Figure: Algorithm 4 Grassmannian K-Means

For our clustering experiments, we test on the first 200 tracklets in the Mind’s eye dataset, which have a total of 24 unique labels. Given the success of DeEPCA (Ye & Zhang, ) in our benchmarks as a runner-up, we use the centralized version (the block power method) in these experiments as well. We test the performance of $K$-means with four different averaging algorithms, namely RGrAv, the power method, the Fréchet mean, and the flag mean. The distance operation is chosen to be the chordal distance for computational efficiency.
The Fréchet mean is computed via iterated gradient descent on the sum of squared distances cost function. Similar to Marrinan et al. (), we find that the four averaging algorithms produce clusters with similar quality across various values for $K$ (these results are not shown for brevity). However as can be seen in Figure 3, the runtime varies significantly between algorithms. The Fréchet mean has the slowest runtime regardless of number of clusters, while RGrAv offers a $2\times$-$10\times$ speedup over the other averaging algorithms for all $K$ values.

## 6 Reproducibility Statement

Section 5.1 is intended to be a comprehensive description sufficient for reproducibility; however, in addition all experiments from this section may be reproduced by running the scripts/decentralized_grassmannian_averaging.py script in the supplemental code; precise instructions and data formatting are described in the header of this file. The $K$-means experiment can be reproduced in three steps. First is by downloading the SUMMET dataset and putting it in the right subdirectory. Simply follow instructions in data_sources/video_separation.py. Next to run the experiment, return to the base directory and run the command python -m scripts.tracklet_clustering to run it as a module. Finally once this completes, there should be a .pkl file with the results. To create the visualization seen in this paper, run python -m vis.visualize_tracklet and look in the new plots/ subdirectory for the results.

## Appendix A Additional Algorithms

Figure: Algorithm 5 Finite RGrAv (Rapid Grassmannian Averaging)

Figure: Algorithm 6 Finite DRGrAv (Decentralized Rapid Grassmannian Averaging)

Figure: Algorithm 7 ChebyshevCoefficients

Figure: Algorithm 8 ChebyshevRoot

## Appendix B Auxiliary Theorems

###### Lemma 2 .

Suppose $f^{\star}\in\mathcal{P}_{t}^{\prime}$ is a solution to eq. 2. Then the $f^{\star}$ has no roots in $[\alpha,1]$.

###### Proof.

If $f^{\star}$ had a root in $[\beta,1]$, then $\min_{\lambda\in\left[\beta,1\right]}\left\lvert f^{\star}\left(\lambda\right)
\right\rvert=0$ and the objective function is unbounded, so $f^{\star}$ cannot have any roots in $[\beta,1]$. This also allows us to conclude $f^{\star}$ is not the 0 0 polynomial, which is used in what follows.
Now, we show that if $f^{\star}$ has a root in $(\alpha,\beta)$, moving this root to $\alpha$ strictly decreases the objective function’s value showing that $f^{\star}$ could not have been a solution to the minization problem.
If $f^{\star}$ had a root in $(\alpha,\beta)$, then we can write $f^{\star}(\lambda)=(\lambda-\lambda_{1})g(\lambda)$ for some $\lambda_{1}\in(\alpha,\beta)$ and $g\in\mathcal{P}_{t-1}^{\prime}$ where $g$ is not the 0 0 polynomial. It follows that

$$ $\displaystyle\frac{\max_{\lambda\in\left[0,\alpha\right]}\left\lvert f^{\star} \left(\lambda\right)\right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left \lvert f^{\star}\left(\lambda\right)\right\rvert}$ $\displaystyle=\frac{\max_{\lambda\in\left[0,\alpha\right]}\left\lvert(\lambda- \lambda_{1})g(\lambda)\right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left \lvert(\lambda-\lambda_{1})g(\lambda)\right\rvert}$ $\displaystyle=\frac{\sup_{\lambda\in[0,\alpha)}\left\lvert(\lambda-\alpha)g( \lambda)\right\rvert\left\lvert\frac{\lambda-\lambda_{1}}{\lambda-\alpha} \right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left\lvert(\lambda-\alpha) g(\lambda)\right\rvert\left\lvert\frac{\lambda-\lambda_{1}}{\lambda-\alpha} \right\rvert}$ $\displaystyle>\frac{\max_{\lambda\in\left[0,\alpha\right]}\left\lvert(\lambda- \alpha)g(\lambda)\right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left \lvert(\lambda-\alpha)g(\lambda)\right\rvert}$ $$

where the final line results from the fact that $\left\lvert\frac{\lambda-\lambda_{1}}{\lambda-\alpha}\right\rvert<1$ for $\lambda\in[\beta,1]$ and $\left\lvert\frac{\lambda-\lambda_{1}}{\lambda-\alpha}\right\rvert>1$ for $\lambda\in[0,\alpha)$. Now we observe that the function $h(\lambda):=(\lambda-\alpha)g(\lambda)\in\mathcal{P}_{t}^{\prime}$ achieves a strictly lower value for the objective function via moving the root $\lambda_{1}\in(\alpha,\beta)$ to $\lambda=\alpha$. Thus, $f^{\star}$ would not be a solution to eq. 2 and therefore $f^{\star}$ cannot have roots in $(\alpha,\beta)$.
Now we show that if $f^{\star}$ had a root at $\lambda=\alpha$, we could slightly move the root to some point to the left of $\alpha$ and decrease the objective function’s value.
Suppose $f^{\star}$ has a root at $\lambda=\alpha$. We can write $f^{\star}(\lambda)=(\lambda-\alpha)g(\lambda)$ for some $g\in\mathcal{P}_{t-1}^{\prime}$, where $g$ is not the 0 0 polynomial. Let $\delta>0$ be such that

$$ $\displaystyle\max_{\lambda\in[\alpha-\frac{\delta}{2},\alpha]}\lvert f^{\star} (\lambda)\rvert<\min_{\epsilon\in[0,\alpha]}\max_{\lambda\in[0,\alpha]}\lvert( \lambda-\epsilon)g(\lambda)\rvert$ $$

which exists since $\min_{\epsilon\in[0,\alpha]}\max_{\lambda\in[0,\alpha]}\lvert(\lambda-\epsilon
)g(\lambda)\rvert>0$ (as $g\not\equiv 0$ from the beginning of the proof) and $f^{\star}$ is continuous. It follows that

$$ $\displaystyle\frac{\max_{\lambda\in\left[0,\alpha\right]}\left\lvert f^{\star} \left(\lambda\right)\right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left \lvert f^{\star}\left(\lambda\right)\right\rvert}$ $\displaystyle=\frac{\max_{\lambda\in\left[0,\alpha\right]}\left\lvert(\lambda- \alpha)g(\lambda)\right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left \lvert(\lambda-\alpha)g(\lambda)\right\rvert}$ $\displaystyle=\frac{\sup_{\lambda\in[0,\alpha]\setminus\{\alpha-\delta\}}\left \lvert(\lambda-(\alpha-\delta))g(\lambda)\right\rvert\left\lvert\frac{\lambda- \alpha}{\lambda-(\alpha-\delta)}\right\rvert}{\min_{\lambda\in\left[\beta,1 \right]}\left\lvert(\lambda-(\alpha-\delta))g(\lambda)\right\rvert\left\lvert \frac{\lambda-\alpha}{\lambda-(\alpha-\delta)}\right\rvert}$ $\displaystyle=\frac{\max\{\sup_{\lambda\in[0,\alpha-\frac{\delta}{2})\setminus \{\alpha-\delta\}}\left\lvert(\lambda-(\alpha-\delta))g(\lambda)\right\rvert \left\lvert\frac{\lambda-\alpha}{\lambda-(\alpha-\delta)}\right\rvert,\Phi\}}{ \min_{\lambda\in\left[\beta,1\right]}\left\lvert(\lambda-(\alpha-\delta))g( \lambda)\right\rvert\left\lvert\frac{\lambda-\alpha}{\lambda-(\alpha-\delta)} \right\rvert}$ $\displaystyle>\frac{\max_{\lambda\in[0,\alpha]}\left\lvert(\lambda-(\alpha- \delta))g(\lambda)\right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left \lvert(\lambda-(\alpha-\delta))g(\lambda)\right\rvert}$ $$

with $\Phi=\max_{\lambda\in[\alpha-\frac{\delta}{2},\alpha]}\left\lvert(\lambda-(
\alpha-\delta))g(\lambda)\right\rvert\left\lvert\frac{\lambda-\alpha}{\lambda-
(\alpha-\delta)}\right\rvert$. The final line results from the fact that $\left\lvert\frac{\lambda-\alpha}{\lambda-(\alpha-\delta)}\right\rvert<1$ for $\lambda\in[\beta,1]$, $\left\lvert\frac{\lambda-\alpha}{\lambda-(\alpha-\delta)}\right\rvert\geq 1$ for $\lambda\in[0,\alpha-\frac{\delta}{2})$, and since $\left\lvert\frac{\lambda-\alpha}{\lambda-(\alpha-\delta)}\right\rvert\leq 1$ for $\lambda\in[\alpha-\frac{\delta}{2},\alpha]$,

$$ $\displaystyle\max_{\lambda\in[\alpha-\frac{\delta}{2},\alpha]}\left\lvert( \lambda-(\alpha-\delta))g(\lambda)\right\rvert\left\lvert\frac{\lambda-\alpha} {\lambda-(\alpha-\delta)}\right\rvert$ $\displaystyle<\min_{\epsilon\in[0,\alpha]}\max_{\lambda\in[0,\alpha]}\lvert( \lambda-\epsilon)g(\lambda)\rvert*1$ $\displaystyle\leq\max_{\lambda\in[0,\alpha]}\left\lvert(\lambda-(\alpha-\delta ))g(\lambda)\right\rvert$ $$

Now we observe that the function $h(\lambda):=(\lambda-(\alpha-\delta))g(\lambda)\in\mathcal{P}_{t}^{\prime}$ achieves a strictly lower value for the objective function via moving the root from $\alpha$ to $\lambda=\alpha-\delta$. Thus, $f^{\star}$ would not be a solution to eq. 2 and therefore $f^{\star}$ cannot have a root at $\lambda=\alpha$.
∎

###### Definition 1 .

A polynomial $f$ is $t$-equioscillatory on an interval
$\left[a,b\right]$ if there exists some $a\leq\gamma_{0}<\gamma_{1}<\dots<\gamma_{t-1}\leq b$ such that $\left\lvert f\left(\gamma_{s}\right)\right\rvert=\max_{\lambda\in\left[a,b
\right]}\left\lvert f\left(\lambda\right)\right\rvert$ for all $0\leq s\leq t-1$ and $f\left(\gamma_{0}\right)=-f\left(\gamma_{1}\right)=f\left(\gamma_{2}\right)=-f
\left(\gamma_{3}\right)=\dots$.

###### Lemma 3 .

A solution to the minimization problem

$$ $\displaystyle\operatorname*{minimize}_{f_{t}\in\mathcal{P}_{t}^{\prime}}\frac{ \max_{\lambda\in\left[0,\alpha\right]}\left\lvert f_{t}\left(\lambda\right) \right\rvert}{\min_{\lambda\in\left[\beta,1\right]}\left\lvert f_{t}\left( \lambda\right)\right\rvert}$ $$

must be $t$-equioscillatory.

###### Proof.

Suppose $f^{\star}$ is a minimizer in to eq. 2 that is not $t$-equioscillatory. First, we assume that without loss of generality, $f^{\star}(\alpha)>0$ since if $f^{\star}(\alpha)<0$, we could replace $f^{\star}$ with $-f^{\star}$ and it would still be a minimizer to the eq. 2, and we cannot have $f^{\star}(\alpha)=0$ by Lemma 2.
Let $\{\lambda_{1}=0,\lambda_{2},\cdots,\lambda_{m}\}\in[0,\alpha]$ for some $m\leq t$ denote the distinct roots of $f^{\star}$ in increasing order that lie in $[0,\alpha]$. We have $m$-intervals, $\mathcal{I}=\{I_{i}\}_{i=1}^{m}$ where $I_{i}=[\lambda_{i},\lambda_{i+1}]$ and $\lambda_{m+1}=\alpha$. We have that for any $i\in[m]$, $\operatorname{sgn}\left(f^{\star}(\lambda)\right)=c$ for every $\lambda\in\text{int}(I_{i})$ where $c\in\{-1,1\}$. In other words, in the interior of any intervals from $\mathcal{I}$, $f^{\star}$ takes strictly all positive values or strictly all negative values. Define, for any interval $I\subseteq\mathbb{R}$ and function $g$

$$ $\displaystyle\operatorname{sgn}\left(g|_{I}\right)=c$ $$

where $c\in\{-1,1\}$ is the value such that $\operatorname{sgn}\left(f^{\star}(\lambda)\right)=c$ for every $\lambda\in\text{int}(I)$. If such a value does not exist, $\operatorname{sgn}\left(g|_{I}\right)$ is left undefined.
We also define for any compact $A\subseteq\mathbb{R}$ and $g\in C(\mathbb{R})$,

$$ $\displaystyle\|g\|_{A}:=\max_{x\in A}\lvert g(x)\rvert$ $$

Let

$$ $\displaystyle L$ $\displaystyle=\{I\in\mathcal{I}:\max_{\lambda\in I}\lvert f^{\star}(\lambda) \rvert<\|f^{\star}\|_{[0,\alpha]}\}$ $\displaystyle M$ $\displaystyle=\{I\in\mathcal{I}:\max_{\lambda\in I}\lvert f^{\star}(\lambda) \rvert=\|f^{\star}\|_{[0,\alpha]}\}=\mathcal{I}\setminus L$ $\displaystyle J$ $\displaystyle=\bigcup_{I\in L}I$ $$

Define $g(\lambda)=\lvert f^{\star}(\lambda)\rvert-\|f^{\star}\|_{[0,\alpha]}$ and denote $\epsilon=\min_{\lambda\in J}\lvert g(\lambda)\rvert$, i.e the minimum distance by which any point in $J$ misses one of $\pm\|f^{\star}\|_{[0,\alpha]}$.
Let $k=\lvert\{\lambda\in[0,\alpha]:f^{\star}(\lambda)=\|f\|_{[0,\alpha]}\}\rvert$. Note that $k<t$ by assumption. Let $M=\{M_{1},\cdots,M_{k}\}$ be intervals listed from left to right where $i_{j}\in[m]$ are indices such that $M_{j}=I_{i_{j}}=[\lambda_{i_{j}},\lambda_{i_{j}+1}]$ for $j\in[k]$ and $\lambda_{i_{j}}<\lambda_{i_{m}}$ for $1\leq j<m\leq k$.
Now define

$$ $\displaystyle R=\{\lambda_{i_{j}}:\operatorname{sgn}\left(f^{\star}|_{M_{j}} \right)\neq\operatorname{sgn}\left(f^{\star}|_{M_{j-1}}\right)\text{ for some }2\leq j\leq k\}$ $$

We have that since $k\leq t-1$ by assumption, then $\lvert R\rvert\leq k-1\leq t-2$.

With $R=\{r_{1},\cdots,r_{q}\}$ for some $q\leq t-2$, we define a polynomial $r:\mathbb{R}\to\mathbb{R}$ based on the sign of $f^{\star}$ on $M_{k}$.
Case A: If $\operatorname{sgn}\left(f^{\star}|_{M_{k}}\right)=1$, we define

$$ $\displaystyle r(\lambda):=c_{r}\lambda(\lambda-\beta)\Pi_{i=1}^{q}(\lambda-r_{ i})$ $$

Case B: Otherwise (when $\operatorname{sgn}\left(f^{\star}|_{M_{k}}\right)=-1$), define

$$ $\displaystyle r(\lambda):=c_{r}\lambda(\lambda-1)\Pi_{i=1}^{q}(\lambda-r_{i})$ $$

In either case, $r$ be a polynomial of degree at most $t$. We now select $c_{r}\in\mathbb{R}$ so that $\|r(\lambda)\|_{[0,\alpha]}<\epsilon$ and set the $\operatorname{sgn}\left(c_{r}\right)$ so that $\operatorname{sgn}\left(r|_{M_{1}}\right)=-\operatorname{sgn}\left(f^{\star}|_
{M_{1}}\right)$. We now show $\operatorname{sgn}\left(r|_{M_{j}}\right)=-\operatorname{sgn}\left(f^{\star}|_
{M_{j}}\right)$ for every $j\in[k]$ via induction. Our base case holds via how we set $c_{r}$ above.
Suppose inductively that $\operatorname{sgn}\left(r|_{M_{j-1}}\right)=-\operatorname{sgn}\left(f^{\star}
|_{M_{j-1}}\right)$ for some $j\leq k$. Then, if $\operatorname{sgn}\left(f^{\star}|_{M_{j}}\right)=\operatorname{sgn}\left(f^{
\star}|_{M_{j-1}}\right)$, we have that as no root was added to $r$ between these intervals and so

$$ $\displaystyle\operatorname{sgn}\left(r|_{M_{j}}\right)=\operatorname{sgn}\left (r|_{M_{j-1}}\right)=-\operatorname{sgn}\left(f^{\star}|_{M_{j-1}}\right)$ $$

with the last equality by the induction hypothesis. Otherwise, if $\operatorname{sgn}\left(f^{\star}|_{M_{j}}\right)\neq\operatorname{sgn}\left(f
^{\star}|_{M_{j-1}}\right)$, then since a root is added between them

$$ $\displaystyle\operatorname{sgn}\left(r|_{M_{j}}\right)=-\operatorname{sgn} \left(r|_{M_{j-1}}\right)=\operatorname{sgn}\left(f^{\star}|_{M_{j-1}}\right)= -\operatorname{sgn}\left(f^{\star}|_{M_{j}}\right)$ $$

So in both cases, we have $\operatorname{sgn}\left(r|_{M_{j}}\right)=-\operatorname{sgn}\left(f^{\star}|_
{M_{j}}\right)$ which closes the induction.
Now we obtain that for any $j\in[k]$,

$$ $\displaystyle\|f^{\star}+r\|_{M_{j}}<\|f^{\star}\|_{[0,\alpha]}$ $$

as, for $M_{j}\neq[\lambda_{m},\alpha]$, the maximum $\|f^{\star}\|_{[0,\alpha]}$ is achieved by $\lvert f^{\star}\rvert$ in the interior of $M_{j}$ where $\lvert r\rvert>0$ as $r$ has no roots in this interior and $\operatorname{sgn}\left(r|_{M_{j}}\right)=-\operatorname{sgn}\left(f^{\star}|_
{M_{j}}\right)$. Additionally, one recognizes that if $M_{j}=[\lambda_{m},\alpha]$, then $r(\lambda)<0$ for $x\in(\lambda_{m},\alpha]$ and so the above bound still holds even when the maximum is attained on the boundary of the interval, i.e when$f^{\star}(\alpha)=\|f^{\star}\|_{[0,\alpha]}$. Furthermore, we have that

$$ $\displaystyle\|f^{\star}+r\|_{J}$ $\displaystyle\leq\max_{\lambda\in J}\lvert f^{\star}(\lambda)\rvert+\|r\|_{[0, \alpha]}$ $\displaystyle<\max_{\lambda\in J}\lvert f^{\star}(\lambda)\rvert+\epsilon=\|f^ {\star}\|_{[0,\alpha]}$ $$

by $\epsilon$’s definition. Now we have obtained

$$ $\displaystyle\|f^{\star}+r\|_{[0,\alpha]}<\|f^{\star}\|_{[0,\alpha]}$ (3) $$

Now we turn to bounding the denominator of the objective function using $f^{\star}$ from above.
Note that $\operatorname{sgn}\left(f^{\star}|_{[\alpha,1]}\right)=1$ as $f^{\star}(\alpha)>0$ and $f^{\star}$ has no roots in $[\alpha,1]$ by Lemma 2. We will now show that in either case,

$$ $\displaystyle\min_{\lambda\in[\beta,1]}\lvert(f^{\star}+r)(\lambda)\rvert\geq \min_{\lambda\in[\beta,1]}\lvert f^{\star}(\lambda)\rvert$ (4) $$

Case A: Since $\operatorname{sgn}\left(f^{\star}|_{M_{k}}\right)=1$, $\operatorname{sgn}\left(r|_{M_{k}}\right)=-1$ which gives $\operatorname{sgn}\left(r|_{[\lambda_{i_{k}},\beta]}\right)=-1$ as $r$ has no roots in $(\lambda_{i_{k}},\beta)$. Since $r$ has a simple root at $\lambda=\beta$ with no other roots greater than $\lambda=\beta$, we have that $\operatorname{sgn}\left(r|_{[\beta,1]}\right)=1$. It immediately follows that

$$ $\displaystyle\min_{\lambda\in[\beta,1]}\lvert(f^{\star}+r)(\lambda)\rvert\geq \min_{\lambda\in[\beta,1]}\lvert f^{\star}(\lambda)\rvert$ $$

Case B: Since $\operatorname{sgn}\left(f^{\star}|_{M_{k}}\right)=-1$, $\operatorname{sgn}\left(r|_{M_{k}}\right)=1$ which gives $\operatorname{sgn}\left(r|_{[\lambda_{i_{k}},1]}\right)=1$ as $r$ has no roots in $(\lambda_{i_{k}},1)$. So $\operatorname{sgn}\left(r|_{[\beta,1]}\right)=1$ and it follows that

$$ $\displaystyle\min_{\lambda\in[\beta,1]}\lvert(f^{\star}+r)(\lambda)\rvert\geq \min_{\lambda\in[\beta,1]}\lvert f^{\star}(\lambda)\rvert$ $$

Combining eq. 3 and eq. 4 gives

$$ $\displaystyle\frac{\|f^{\star}\|_{[0,\alpha]}}{\min_{\lambda\in[\beta,1]} \lvert f^{\star}(\lambda)\rvert}\geq\frac{\|f^{\star}+r\|_{[0,\alpha]}}{\min_{ \lambda\in[\beta,1]}\lvert(f^{\star}+r)(\lambda)\rvert}$ $$

We note that $(f^{\star}+r)(0)=0$ and $f^{\star}+r$ is an at most $t$-degree polynomial and therefore $f^{\star}+r$ is a feasible function for the minimization problem eq. 2. Thus, $f^{\star}$ is not a solution to eq. 2 as $f^{\star}+r$ is feasible and achieves a strictly smaller value for the objective function, which proves the lemma.
∎

###### Lemma 4 .

Define $\mathcal{P}_{t}^{\prime\prime}:=\left\{f_{t}\in\mathcal{P}_{t}~{}|~{}f_{t}
\left(0\right)=0,f_{t}\text{ is
$t$-equioscillatory on }\left[0,\alpha\right],f_{t}\left(\beta\right)=1\right\}$ for $\beta>\alpha$. The problem

$$ $\operatorname*{minimize}_{f_{t}\in\mathcal{P}_{t}^{\prime\prime}}\max_{\lambda \in\left[0,\alpha\right]}\left\lvert f_{t}\left(\lambda\right)\right\rvert$ $$

is solved by

$$ $\displaystyle f_{t}^{\star}\left(\lambda\right)$ $\displaystyle\phantom{:}=\prod_{s=0}^{t-1}\frac{\lambda-r_{s,t}}{\beta-r_{s,t}}$ $\displaystyle r_{s,t}$ $\displaystyle:=\alpha\frac{\cos\left(\frac{\pi\left(s+1/2\right)}{t}\right)+ \cos\left(\frac{\pi}{2t}\right)}{1+\cos\left(\frac{\pi}{2t}\right)}$ $$

###### Proof.

First note that $f_{t}^{\star}$ is indeed feasible; from inspection one
realizes $f_{t}^{\star}\left(\beta\right)=1$ and $r_{t-1,t}=0$ (consequently
$f_{t}^{\star}\left(0\right)=0$) and as a variant of a Chebyshev polynomial (of the
first kind), $f_{t}^{\star}$ is $t$-equioscillatory on $\left[0,\alpha\right]$ with extremal points

$$ $\gamma_{s}=\alpha\frac{\cos\left(\frac{\pi s}{t}\right)+\cos\left(\frac{\pi}{2 t}\right)}{1+\cos\left(\frac{\pi}{2t}\right)},\quad 0\leq s\leq t-1$ $$

Assume for the sake of contradiction that there exists some $g\in\mathcal{P}_{t}^{\prime\prime}$
such that $\max_{\lambda\in\left[0,\alpha\right]}\left\lvert g\left(\lambda\right)\right
\rvert<\max_{\lambda\in\left[0,\alpha\right]}\left\lvert f_{t}^{\star}\left(
\lambda\right)\right\rvert$. This would
then imply that the difference polynomial $h\left(\lambda\right):=f_{t}^{\star}\left(\lambda\right)-g\left(\lambda\right)$ satisfies the following

$$ $\displaystyle\forall s=0,2,\dots~{}:\quad h\left(\gamma_{s}\right)$ $\displaystyle=f_{t}^{\star}\left(\gamma_{s}\right)-g\left(\gamma_{s}\right)$ $\displaystyle>f_{t}^{\star}\left(\gamma_{s}\right)-\max_{\lambda\in\left[0, \alpha\right]}\left\lvert f_{t}^{\star}\left(\lambda\right)\right\rvert$ $\displaystyle=0$ $\displaystyle\forall s=1,3,\dots~{}:\quad h\left(\gamma_{s}\right)$ $\displaystyle=f_{t}^{\star}\left(\gamma_{s}\right)-g\left(\gamma_{s}\right)$ $\displaystyle<f_{t}^{\star}\left(\gamma_{s}\right)+\max_{\lambda\in\left[0, \alpha\right]}\left\lvert f_{t}^{\star}\left(\lambda\right)\right\rvert$ $\displaystyle=0$ $\displaystyle h\left(0\right)$ $\displaystyle=f_{t}^{\star}\left(0\right)-g\left(0\right)$ $\displaystyle=0$ $\displaystyle h\left(\beta\right)$ $\displaystyle=f_{t}^{\star}\left(\beta\right)-g\left(\beta\right)$ $\displaystyle=0$ $$

Since $h$ changes sign on every $\gamma_{s}$, it must have at least one root
in each of the $t-1$ intervals between consecutive $\gamma_{s}$. By
construction, $h$ also has a 2 more roots at 0 0 and $\beta$, meaning $h$
has at minimum $t+1$ distinct roots. However, $h$ is a $t$th order
polynomial, leading to a contradiction.
∎

###### Lemma 2 .

###### Proof.

###### Definition 1 .

###### Lemma 3 .

###### Proof.

###### Lemma 4 .

###### Proof.

### B.1 Proof of Theorem 1

Proof. First, by Lemma 3, since any minimizer $f^{\star}$ of eq. 2 is $t$-equiosciallatory, all $t$ of $f^{\star}$’s roots lie in $[0,\alpha]$. This implies that $f^{\star}$ is increasing in $[\beta,1]$ and therefore we have

$$ $\min_{\lambda\in[\beta,1]}\lvert f^{\star}(\lambda)\rvert=f^{\star}(\beta)$ $$

which is assumed to be positive without loss of generality also as in Lemma 3. Now we can scale $f^{\star}$ so that $f(\beta)=1$ without changing the value of the objective function, i.e

$$ $\displaystyle\frac{\max_{\lambda\in\left[0,\alpha\right]}\left\lvert f^{\star} \left(\lambda\right)\right\rvert}{f^{\star}(\beta)}=\frac{\max_{\lambda\in \left[0,\alpha\right]}\left\lvert Cf^{\star}\left(\lambda\right)\right\rvert}{ Cf^{\star}(\beta)}=\max_{\lambda\in\left[0,\alpha\right]}\left\lvert Cf^{\star }\left(\lambda\right)\right\rvert$ $$

where $C=\frac{1}{f^{\star}(\beta)}$. We have now transformed the problem into

$$ $\operatorname*{minimize}_{f_{t}\in\mathcal{P}_{t}^{\prime\prime}}\max_{\lambda \in\left[0,\alpha\right]}\left\lvert f_{t}\left(\lambda\right)\right\rvert$ $$

as in Lemma 4 which is solved by

$$ $\displaystyle f_{t}^{\star}\left(\lambda\right)$ $\displaystyle\phantom{:}=\prod_{s=0}^{t-1}\frac{\lambda-r_{s,t}}{\beta-r_{s,t}}$ $\displaystyle r_{s,t}$ $\displaystyle:=\alpha\frac{\cos\left(\frac{\pi\left(s+1/2\right)}{t}\right)+ \cos\left(\frac{\pi}{2t}\right)}{1+\cos\left(\frac{\pi}{2t}\right)}$ $$

as desired.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: Additional issues and questions:
- Eq. 1: It would be good to add parantheses to clarify whether summation includes $UU^T$ to avoid potential confusion.
- Could you elaborate on the connection between Theorem 1 and the equioscillation theorem? Would this reproduce the best polynomial approximation of the step function?
- Please explain how the DeEPCA algorithm is adapted for decentralized Grassmannian averaging.
- Line 172: Unrolling the power method loop—could you clarify if this is straightforward or requires further justification?
- Line 198: Are you assuming that $\alpha$ and $\beta$ are known?
- Figure 1: Please elaborate on what is shown and its significance. Additionally, Figure 1 is not referenced in the main text.
- Line 264: "focuses on minimizing the error in the gradient tracking procedure" — please provide more details.
- Table 1: Are these results obtained under the hypercube or cycle setup?
- Line 438: How is the true IAM average computed? Could you include timing for this computation as well?
- Section 5.2: Please clarify that this experiment does not test the decentralized algorithm.
- $\alpha$ is used as a bound on eigenvalues in Section 4 and later as a step size parameter in Section 5, which could be confusing. Please consider revising this notation.

Additional minor issues:
- Line 51: "subsequent projections may be applied locally thereafter)." - could you please elaborate?
- Line 62: "We demonstrate merit through a theoretical guarantee on the optimality" - Theorem 1 proves optimality of the choice for the polynomial $f_t$. Does this guarantee the optimality of your approach? In what sense?
- Line 72: Define Grassmannian and Stiefel manifolds, as these terms and notations are introduced here for the first time.
- line 377: What is the parameter $\rho$?
soundness: 3
strengths: - The paper addresses an interesting problem and is generally well-written, presenting the proposed approach clearly, though it occasionally lacks detail.
- The authors prove that Chebyshev polynomials are optimal for their modified power method.
- Their experiments demonstrate that the algorithm outperforms competing methods in synthetic Grassmannian averaging and K-means clustering of subspaces spanned by video sequences for motion clustering. They have also successfully adapted related methods (such as DeEPCA) to the problem of Grassmannian averaging, which can be seen as part of their contribution.
- The authors provide code that implements the proposed approach and reproduces the experiments presented in the paper (although I have not tested the code).
summary: This paper addresses the problem of computing the induced arithmetic mean of subspaces (points on the Grassmannian manifold). The authors focus on a decentralized setting where each agent has limited data and restricted communication with its neighbors. Their approach is based on spectral analysis of the associated optimization problem, resulting in a power-method-like algorithm. The utility of the method is demonstrated through synthetic experiments and an application to clustering video motion data via Grassmannian K-means clustering.
weaknesses: 1. Importance of Decentralized Grassmannian Averaging: The paper emphasizes the importance of decentralized averaging, which the authors claim is significant. However, the paper does not provide concrete examples where decentralization is essential and demonstrates it only through a single set of synthetic experiments (emulating a decentralized setup). Could you provide more compelling examples to illustrate the importance of this problem?

2. Comparison and Relation to Existing Literature: How does RGrAv compare or relate to existing methods, such as block power methods (e.g., Chebyshev-Davidson, spectral filtering) or other similar approaches? Given that the problem in Eq. (1) is reduced to a spectral problem in Section 3.1, a more thorough exploration of relevant methods, both centralized and decentralized, would be beneficial. (I am not an expert in this area, but there appears to be relevant literature available.)

3. Need for Additional Explanations: While the paper is overall well-written, certain choices or statements would benefit from further explanation or justification. I was able to understand some of these upon reflection, while others remain unclear. For example:
  - Section 3.1: The definition of the optimization problem (1) and the equivalent spectral perspective are presented very densley, without references, which could be difficult for readers unfamiliar with the topic. For instance, why is (1) invariant to the choice of representatives $[U_m]$? Why is the projection on the K largest eigenvectors $[V]$ the minimizer of (1)? Additional explanations or citations would be helpful.
  - Section 3.2: What is average consensus (AC)? Without prior knowledge, the proposed procedure is unclear. Why can’t AC be used directly to approximate (1)? Most importantly, why is it sufficient to approximate the lower dimensional $PX$? In what way does this guarantee an approximation of $P$?
  - Section 4.1: How are the optimal polynomials $f_t$ used to recover the space $U$? This is again central to the approach, but is challenging to understand. Theorem 1 proves the optimality of $f_t$ in the sense of Eq. (2), but in what sense does this translate to the optimality of the approximate average space $[U]$?
  - Section 4.2: "fortunately, $f_t$ is well-approximated in terms of $f_{t-1}^*$ and $f_{t-2}^*$, in what sense does this hold? Is it an experimental observation or is there a theoretical justification?

4. Experimental Results Are Limited:
  - The experiments in Section 5.1 focus on the decentralized setting and raise a few questions:
    - In these experiments, DRGrAv is "sub-optimally tuned" (Line 353) while parameters for other methods are carefully tuned through parameter search (commendable effort). How would DRGrAv perform if optimally tuned?
    - The experiments are conducted for a single setup (M,N,K)=(64,150,30) with two types of agent connectivity (hypercube and cycle). Is this setup representative?
    - Are similar results obtained against competitive methods in different setups?
    - Is this a setup where decentralization is necessary or preferable? How would DRGrAv perform compared to competitors in a scenario with very large M or N, where decentralization is crucial?
    - Additional experiments supporting the importance and benefits of the proposed approach for decentralized averaging would significantly strengthen the manuscript, given that this is a primary focus of the paper.

  - The experiment in Section 5.2 uses Grassmannian averaging to define Grassmannian K-means for motion clustering in videos. The authors present only timing results but do not include any assessment of clustering quality for each method. The statement that "the four averaging algorithms produce clusters with similar quality across various values for K  (these results are not shown for brevity)" is unclear. Showing these results alongside runtime comparisons seems essential for evaluating each method's overall effectiveness. These results could be briefly mentioned in the main text and detailed in the appendix if space is limited. Does RGrAv produce the same results as Flag? How does RGrAv compare to the block power method, ground truth IAM? Does the spectrum gap assumption (Section 4.1) hold in these experiments?

  - The evidence for a performance advantage is somewhat limited. The decentralized results in Section 5.1 are comparable to DeEPCA, while the centralized results in Section 5.2 are closely followed by the block power method.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: omitted.
soundness: 2
strengths: omitted.
summary: The paper proposes a method for computing the average of some subspaces in the standard and decentralized setting.
weaknesses: - The contribution at the technical level is incremental, in my opinion. Theorem 1 seems to be a basic variant of the known result for Chebyshev polynomials. The only difference seems to be that the paper adds an extra constraint $f_t(0)=0$.

- The assumption on the dual-banded spectrum is also not new. This assumption is used in the AISTATS 2022 paper: Super-Acceleration with Cyclical Step-sizes (https://proceedings.mlr.press/v151/goujaud22a/goujaud22a.pdf); see Eq. 3 there. That paper has a much deeper exploration of related aspects than the present manuscript. If NeurIPS holds the same high standard, the lack of technical depth forms a ground for rejection.

- From a practical viewpoint, the proposed method performs similarly to the baselines it compares, and there is no significant advantage in terms of acceleration or efficiency, as shown in Section 5.1

The above two points suggest Lines 60 - 62 are never-before-seen overstatements: "*Our algorithms operate similarly to the famous power method, with the distinction that Chebyshev polynomials are employed to leverage a “dual-banded” property of the problem in order to achieve never-before-seen efficiency in computation and communication.*"

- The algorithm listings (Algorithms 1 and 2) are highly repeated.
- Section 5.2 is disconnected from the prior sections of the paper, which puts emphasis on the decentralized setting. From Section 1 to Section 5.1 there is no mention of k-means clustering.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. Can you clarify what is meant by "QR" in the first equation of Sec. 3.3? I know the QR decomposition, but this generally returns *two* values (the Q and the R matrices). In the text, only a single matrix is returned, so it would be good to clarify if this was Q or R.
2. Can you provide a citation for "StableQR"?
3. In the last experiment, you compute a Frechet mean, but use the chordal distance. Isn't the Frechet mean defined according to the geodesic distance? Can you clarify what is actually being done here?
soundness: 4
strengths: The paper observes that the eigenvalues of projection matrices have a particular structure, which can be exploited to build a more efficient variant of the power iteration to compute eigenvalues/vectors. As far as I can tell, this contribution is novel and is potentially more widely applicable than explored in the submitted paper.

The paper is generally well-written and easy to understand.

The developed algorithm is shown to be easy to extend to a decentralized implementation. This is a quite nice practical property.
summary: The paper provides a fast algorithm for computing a particular average on the Grassmann manifold. The algorithm is, in essence, a variant of the classic power iteration, which is tailored for projection matrices (which can, in turn, be used to represent points on the Grassmannian). The algorithm is suitable for decentralized computations, making it easy to scale across machines. Empirical results show that the algorithm is faster than a selection of baseline algorithms while providing comparable results.
weaknesses: My main concern with the paper is the studied problem. The paper states that points on the Grassmann manifold are "used ubiquitously in machine learning, computer vision, and signal processing". It is true that some years ago (before deep learning took off), some work went into using subspaces (and hence Grassmannians) to represent batches of data. This approach has, however, not caught on and is quite rare today.

I do not want to argue that a paper is uninteresting because it explores a (current) niche approach. However, I consider it problematic that the paper does not demonstrate that the problem of study is relevant to solve. The paper demonstrates that the proposed algorithm is faster than alternatives, but it fails to demonstrate that the algorithm computes a quantity of interest. I would consider the paper significantly more interesting if it was demonstrated that Grassmann Averaging is useful for solving real tasks.

Some further minor comments:
* It seems to me that Theorem 1 is more generally applicable than just computing Grassmann averages. If this understanding is correct, it would have been valuable to have pointers in such directions.
* I miss references to early work on using Grassmann Averages to define principal components of data: "Grassmann Averages for Scalable Robust PCA", CVPR 2014. This relies on a crude chordal-like metric.
* While I appreciate that many algorithms are provided in the paper, I miss general explanations of what these algorithms do. On a first reading, it can be difficult to read an algorithm, so I would have appreciated more hand-holding.

Meta-review: This paper introduces an efficient approach for performing subspace averaging on the Grassmannian manifold in a decentralized setting. The authors provide theoretical guarantees of optimality and empirical results showcasing the proposed algorithm’s advantages in terms of accuracy and runtime compared to relevant approaches. Although the paper is well-organized and easy to follow, the authors have not sufficiently emphasized the significance of their contributions. The reviewers raised concerns regarding the practical relevance of the proposed method, the novelty of its technical ideas, as well as the thoroughness of its experimental validation and comparisons with existing solutions. In my view, the paper requires major revisions before being ready for acceptance at ICLR. Specifically, the authors should more clearly motivate their approach, emphasize the importance of the problem, improve the clarity of the technical contributions, and strengthen the empirical validation section. That being said, I recommend rejecting this paper in its current form.

Decision: Reject