Exact Minimax Optimality of Spectral Methods in Phase Synchronization and Orthogonal Group Synchronization

We consider the phase synchronization problem with additive Gaussian noises and incomplete data [2, 5, 39, 18]. Let $z^{*}_{1},\ldots,z^{*}_{n}\in\mathbb{C}_{1}$ where $\mathbb{C}_{1}:=\{x\in\mathbb{C}:\left|z\right|=1\}$, the set of all unit complex numbers. Then each $z^{*}_{j}$ can be written equivalently as $e^{i\theta^{*}_{j}}$ for some phase (or angle) $\theta^{*}_{j}\in[0,2\pi)$ . For each $1\leq j<k\leq n$, the observation $X_{jk}\in\mathbb{C}$ is missing at random. Let $A_{jk}\in\{0,1\}$ and $X_{jk}$ satisfy

$\displaystyle X_{jk}:=\begin{cases}{z_{j}^{*}\overline{z_{k}^{*}}+\sigma W_{jk}},\text{ if }A_{jk}=1,\\ 0,\quad\quad\quad\quad\quad\text{ if }A_{jk}=0,\end{cases}$

(1)

where $A_{jk}\sim\text{Bernoulli}(p)$ and $W_{jk}\sim\mathcal{CN}(0,1)$. That is, each $X_{jk}$ is missing with probability $1-p$ and is denoted as 0. If it is not missing, it is equal to $z_{j}^{*}\overline{z_{k}^{*}}$ with an additive noise $\sigma W_{jk}$ where $W_{jk}$ follows the standard complex Gaussian distribution: ${\rm Re}(W_{jk}),{\rm Im}(W_{jk})\sim\mathcal{N}(0,1/2)$ independently. Each $A_{jk}$ is the indicator of whether $X_{jk}$ is observed or not. We assume all random variables $\{A_{jk}\}_{1\leq j<k\leq n},\{W_{j,k}\}_{1\leq j<k\leq n}$ are independent of each other. The goal is to estimate the phase vector $z^{*}:=(z_{1}^{*},\ldots,z_{n}^{*})\in\mathbb{C}_{1}^{n}$ from $\{A_{jk}\}_{1\leq j<k\leq n}$ and $\{X_{jk}\}_{1\leq j<k\leq n}$.

The observations can be seen as entries of a matrix $X\in\mathbb{C}^{n\times n}$ with $X_{jj}:=0$ and $X_{kj}:=\overline{X_{jk}}$ for any $1\leq j<k\leq n$. Define $A_{jj}:=0$ and $A_{kj}:=A_{jk}$ for all $1\leq j<k\leq n$. Then the matrix $A\in\{0,1\}^{n\times n}$ can be interpreted as the adjacency matrix of an Erdös-Rényi random graph with edge probability $p$. Define $W\in\mathbb{C}^{n\times n}$ such that $W_{jj}:=0$ and $W_{kj}:=\overline{W_{jk}}$ for all $1\leq j<k\leq n$. Then all the matrices $A,W,X$ are Hermitian and $X$ can be written equivalently as

$\displaystyle X=A\circ\left(z^{*}z^{*{\mathrm{\scriptscriptstyle H}}}+\sigma W\right)=A\circ\left(z^{*}z^{*{\mathrm{\scriptscriptstyle H}}}\right)+\sigma A\circ W.$

(2)

Note that $X$ can be seen as a noisy version of

$\displaystyle\mathbb{E}X=pz^{*}z^{*{\mathrm{\scriptscriptstyle H}}}-pI_{n}$

(3)

whose leading eigenvector is $z^{*}/\sqrt{n}$. This motivates the following spectral method [36, 18, 13]. Let $u\in\mathbb{C}^{n}$ be the leading eigenvector of $X$. Then the spectral estimator $\widehat{z}\in\mathbb{C}_{1}^{n}$ is defined as

$\displaystyle\widehat{z}_{j}:=\begin{cases}\frac{u_{j}}{\left|u_{j}\right|},\text{ if }u_{j}\neq 0,\\ 1,\quad\text{ if }u_{j}=0,\end{cases}$

(4)

for each $j\in[n]$, where each $u_{j}$ is normalized so that $\widehat{z}_{j}\in\mathbb{C}_{1}$. The performance of the spectral estimator can be quantified by a normalized squared $\ell_{2}$ loss

$\displaystyle\ell(\widehat{z},z^{*}):=\min_{a\in\mathbb{C}_{1}}\frac{1}{n}\sum_{j=1}^{n}\left|\widehat{z}_{j}-z_{j}^{*}a\right|^{2},$

(5)

where the minimum over $\mathbb{C}_{1}$ is due to the fact that $z^{*}_{1},\ldots,z^{*}_{n}$ are identifiable only up to a phase.

The spectral estimator $\widehat{z}$ is simple and easy to implement. Regarding its theoretical performance, it was suggested in [18] that an upper bound $\ell(\widehat{z},z^{*})\leq C(\sigma^{2}+1)/(np)$ holds with high probability for some constant $C$. However, the minimax risk of the phase synchronization was established in [18] and has the following lower bound:

$\displaystyle\inf_{z\in\mathbb{C}^{n}}\sup_{z^{*}\in\mathbb{C}_{1}^{n}}\mathbb{E}\ell(z,z^{*})\geq\left(1-o(1)\right)\frac{\sigma^{2}}{2np}.$

(6)

To provably achieve the minimax risk, the spectral method is often used as an initialization for some more sophisticated procedures. For example, it was used to initialize a generalized power method (GPM) [7, 29, 31] in [18]. Nevertheless, numerically the performance of the spectral method is already very good and the improvement from GPM is often marginal. This raises the following questions about the performance of the spectral method: Can we derive a sharp upper bound? Does the spectral method already achieve the minimax risk or not?

In this paper, we provide complete answers to these questions. We carry out a sharp $\ell_{2}$ analysis of the performance of the spectral estimator $\widehat{z}$ and further show it achieves the minimax risk with the correct constant. Our main result is summarized below in Theorem 1 in asymptotic form. Its non-asymptotic version will be given in Theorem 3 that only requires $\frac{np}{\sigma^{2}}$, $\frac{np}{\log n}$ to be greater than a certain constant. We note that in this paper, $p$ and $\sigma^{2}$ are not constants but functions of $n$. This dependence can be more explicitly represented as $p_{n}$ and $\sigma^{2}_{n}$. However, for simplicity of notation and readability, we choose to denote these as $p$ and $\sigma^{2}$ throughout the paper.

Theorem 1 shows that $\widehat{z}$ is not only rate-optimal but also achieves the exact minimax risk with the correct leading constant in front of the optimal rate. The conditions needed in Theorem 1 are necessary for consistent estimation of $z^{*}$ in the phase synchronization problem. The condition $\frac{np}{\sigma^{2}}\rightarrow\infty$ is needed so that $z^{*}$ can be estimated with a vanishing error according to the minimax lower bound (6). The condition $\frac{np}{\log n}\rightarrow\infty$ allows $p$ to decrease as $n$ grows and is close to the $\frac{np}{\log n}>(1+\epsilon)$ condition required for $A$ to be connected. If $A$ has disjoint subgraphs, it is impossible to estimate $z^{*}$ with a global phase. These two conditions are also needed in [18, 19] to establish the optimality of MLE (maximum likelihood estimation), GPM, and SDP (semidefinite programming). Under these two conditions, [18] used the spectral method as an initialization for the GPM and shows GPM achieves the minimax risk after $\log(1/\sigma^{2})$ iterations. On the contrary, Theorem 1 shows that the spectral method already achieves the minimax risk. This means that the spectral estimator is minimax optimal and achieves the correct leading constant whenever consistent estimation is possible, and in this parameter regime, it is as good as MLE, GPM, and SDP.

There are two key novel components toward establishing Theorem 1. The first is a new idea regarding the choice of the “population eigenvector” as $u$ can be viewed as its sample counterpart obtained from data. Due to (3) and the fact $pz^{*}z^{*{\mathrm{\scriptscriptstyle H}}}=np(z^{*}/\sqrt{n})(z^{*}/\sqrt{n})^{\mathrm{\scriptscriptstyle H}}$ is rank-one with the eigenvector $z^{*}/\sqrt{n}$, existing literature such as [18, 20, 27] treated $z^{*}/\sqrt{n}$ as the population eigenvector and study the perturbation of $u$ with respect to $z^{*}/\sqrt{n}$. This seems natural but turns out to be unappealing as it fails to explain why the spectral estimator is able to recover all phases exactly when $\sigma=0$, i.e., when there is no additive noise. Instead, we denote $\widecheck{u}\in\mathbb{R}^{n}$ as the leading eigenvector of $A$ and regard $u^{*}\in\mathbb{C}^{n}$, defined as

$\displaystyle u^{*}:=z^{*}\circ\widecheck{u},$

(8)

i.e., $u^{*}_{j}=z^{*}_{j}\widecheck{u}_{j}$ for each $j\in[n]$, as the population eigenvector. Note that $u^{*}$ is random as it depends on the graph $A$. A careful analysis of $u^{*}$ reveals that it is the leading eigenvector of $A\circ z^{*}z^{*{\mathrm{\scriptscriptstyle H}}}$ (see Lemma 1). In addition, Proposition 1 shows that with high probability, $u^{*}_{j}/|u^{*}_{j}|=z^{*}_{j}$ for each $j\in[n]$, up to a global phase. Since $u$ equals $u^{*}$ when $\sigma=0$, it successfully explains the exact recovery of the spectral method in the no-additive-noise case. Another advantage of viewing $u^{*}$ as the population eigenvector, instead of $z^{*}/\sqrt{n}$, is that intuitively $u^{*}$ is closer to $u$ than $z^{*}/\sqrt{n}$ is. This is because $A\circ z^{*}z^{*{\mathrm{\scriptscriptstyle H}}}$ is closer to the data matrix $X$ than $pz^{*}z^{*{\mathrm{\scriptscriptstyle H}}}$ is.

The second key component is a novel perturbation analysis for $u$. Classical matrix perturbation theory such as Davis-Kahan Theorem focuses on analyzing $\inf_{b\in\mathbb{C}_{1}}\|u-u^{*}b\|$. We go beyond it and show $u$ can be well-approximated by its first-order approximation $\widetilde{u}$ defined as

$\displaystyle\widetilde{u}:=\frac{Xu^{*}}{\left\|{Xu^{*}}\right\|},$

(9)

in the sense that the difference between these two vectors (up to a phase) has a small $\ell_{2}$ norm. This means that when $np\gtrsim\log n$ and $np\gtrsim\sigma^{2}$, we have

$\displaystyle\inf_{b\in\mathbb{C}_{1}}\left\|{u-\widetilde{u}b}\right\|\lesssim\frac{\sigma^{2}+\sigma}{np},$

(10)

with high probability (see Proposition 2). In fact, our perturbation analysis extends beyond the phase synchronization problem. What we establish is a general perturbation theory that can be applied to two arbitrary Hermitian matrices (see Lemma 2), which might be of independent interest.

With the help of these two key components, we then carry out an entrywise analysis for each $\widehat{z}_{j}=u_{j}/\left|u_{j}\right|$. Note that $u_{j}$ can be decomposed into $\widetilde{u}_{j}$ and the difference between $u_{j}$ and $\widetilde{u}_{j}$ (up to some global phase). We can decompose the error of $\widehat{z}_{j}$ into two parts, one is related to the estimation error of $\widetilde{u}_{j}/|\widetilde{u}_{j}|$, and the other is related to the magnitude of the difference between $u_{j}$ and $\widetilde{u}_{j}$ (up to some global phase). Summing over all coordinates, the first part eventually leads to the minimax risk $(1+o(1)){\sigma^{2}}/{2np}$ and the second part is essentially negligible due to (10), which leads to the exact minimax optimality of the spectral estimator.

The observations $\{\mathcal{X}_{jk}\}_{1\leq j<k\leq n}$ can be seen as submatrices of an $nd\times nd$ matrix $\mathcal{X}$ with $\mathcal{X}_{jj}:=0$ and $\mathcal{X}_{kj}:=\mathcal{X}_{jk}^{\mathrm{\scriptscriptstyle T}}$ for any $1\leq j<k\leq n$. Then $\mathcal{X}$ is symmetric and can be seen as a noisy version of

$\displaystyle\mathbb{E}\mathcal{X}=pZ^{*}Z^{*{\mathrm{\scriptscriptstyle T}}}-pI_{nd}$

(12)

whose leading eigenspace is $Z^{*}/\sqrt{n}$, where $Z^{*}\in\mathcal{O}(d)^{n}$ is an $nd\times d$ matrix such that its $j$th submatrix is $Z^{*}_{j}$. Similar to the phase synchronization, we have the following spectral method. Let $\lambda_{1}\geq\ldots\geq\lambda_{d}$ be the largest $d$ eigenvalues of $\mathcal{X}$ and $u_{1},\ldots,u_{d}\in\mathbb{R}^{nd}$ be their corresponding eigenvectors. Denote $U:=(u_{1},\ldots,u_{d})\in\mathbb{R}^{nd\times d}$ as the eigenspace that includes the top $d$ eigenvectors of $\mathcal{X}$. For each $j\in[n]$, denote $U_{j}\in\mathbb{R}^{d\times d}$ as its $j$th submatrix. Then the spectral estimator $\widehat{Z}_{j}\in\mathcal{O}(d)$ is defined as

$\displaystyle\widehat{Z}_{j}:=\begin{cases}\mathcal{P}(U_{j}),\text{ if }\det(U_{j})\neq 0,\\ I_{d},\quad\;\;\;\text{ if }\det(U_{j})=0,\end{cases}$

(13)

for each $j\in[n]$. Here the mapping $\mathcal{P}:\mathbb{R}^{d\times d}\rightarrow\mathcal{O}(d)$ is derived from the polar decomposition and serves as a normalization step for each $U_{j}$ such that $\widehat{Z}_{j}\in\mathcal{O}(d)$. Let $\widehat{Z}\in\mathcal{O}(d)^{n}$ be an $nd\times d$ matrix such that $\widehat{Z}_{j}$ is its $j$th submatrix for each $j\in[n]$. Then the performance of $\widehat{Z}$ can be quantified by a loss function $\ell^{\text{od}}(\widehat{Z},Z^{*})$ that is analogous to (5). The detailed definitions of $\mathcal{P}$ and $\ell^{\text{od}}$ are deferred to Section 3.

The spectral method $\widehat{Z}$ was used as an initialization in [20] for a variant of GPM to achieve the exact minimax risk $(1+o(1))\frac{d(d-1)\sigma^{2}}{2np}$ for $d=O(1)$. To conduct a sharp analysis of its statistical performance, we extend our novel perturbation analysis from analyzing the leading eigenvector to the leading eigenspace. Recall $\widecheck{u}$ is the leading eigenvector of $A$. Analogous to (8), we have a novel choice of the population eigenspace $U^{*}$, defined as

$\displaystyle U^{*}:=Z^{*}\circ(\widecheck{u}\otimes\mathds{1}_{d}),$

(14)

and view $U$ as its sample counterpart. This is different from existing literature [20, 40] which uses $Z^{*}/\sqrt{n}$ as the population eigenspace. Our choice of $U^{*}$ enables the establishment of the exact recovery of the spectral method when there is no additive noise (i.e., $\sigma=0$), as seen Proposition 3, and is closer to $U$ than $Z^{*}/\sqrt{n}$ is.

The first-order approximation of $U$ is a matrix determined by $\mathcal{X}U^{*}$ whose explicit expression will be given later in (22). We then show $U$ can be well-approximated by its first-order approximation, analogous to (9), with a remainder term of a small $\ell_{2}$ norm (see Proposition 4). This is a consequence of a more general eigenspace perturbation theory (see Lemma 4) for two arbitrary Hermitian matrices. Using the first-order approximation, we then carry out an entrywise analysis for $\widehat{Z}$. Our main result for the spectral method in the $\mathcal{O}(d)$ synchronization is summarized below in Theorem 2. The non-asymptotic version will be given in Theorem 4.

Theorem 2 shows that the spectral method $\widehat{Z}$ achieves exact minimax optimality as it matches the minimax lower bound $(1+o(1))\frac{d(d-1)\sigma^{2}}{2np}$ established in [20]. Similar to the phase synchronization, the two conditions needed in Theorem 2 so that consistent estimation of $Z^{*}$ is possible. They are also needed in [20] to achieve the minimax risk by a variant of GPM initialized by the spectral method. On the contrary, Theorem 2 shows that in this parameter regime, the spectral method is already minimax optimal with the correct leading constant.

Fine-grained perturbation analysis of eigenvectors has gained increasing attention in recent years for various low-rank matrix problems in machine learning and statistics. Existing literature has mostly focused on establishing $\ell_{\infty}$ bounds for eigenvectors [1, 12, 15] or $\ell_{2,\infty}$ bounds on eigenspaces [23, 11, 9, 3]. For instance, [1] developed $\ell_{\infty}$-type bounds for the difference between eigenvectors (or eigenspaces) and their first-order approximations. In this paper, we focus on developing sharp $\ell_{2}$-type perturbation bounds, where direct applications of existing $\ell_{\infty}$-type results will result in extra logarithm factors.

For the phase synchronization problem, [17, 21, 32] investigated variants of spectral methods based on Laplacian matrices. Instead of using the leading eigenvector of $X$ as in this paper, they utilize the eigenvector corresponding to the smallest eigenvalue of $D-X$ or $I_{n}-D^{-\frac{1}{2}}XD^{-\frac{1}{2}}$, where $D\in\mathbb{R}^{n\times n}$ is the degree matrix of $A$ with diagonal entries $D_{jj}:=\sum_{k\neq j}A_{jk}$ and off-diagonal entries set to zero. These studies have established upper bounds for the performance of their spectral methods applicable to general graphs $A$ and additive noise $W$. Our focus, however, is on Erdös-Rényi random graphs with Gaussian noise. Under our setting, their results imply an upper bound of $\frac{C\sigma^{2}}{np}$, where $C$ is a constant significantly greater than 1. In contrast, our work establishes a sharp upper bound with the correct leading constant $1/2$. Our analytical approach could potentially be extended to their methods to achieve the correct constant $1/2$.

The existing literature [18, 19, 20] explored the exact minimax risk in synchronization problems, focusing primarily on minimax lower bounds and analyzing MLE, GPM, and SDP. While our study shares a thematic resemblance with these prior efforts, it fundamentally diverges in both analysis and proof techniques. Previous studies hinge on contraction properties of the generalized power iteration (GPI), demonstrating the iterative reduction in GPM error until an optimal error is achieved. This approach further interprets MLE as a GPI fixed point and SDP as an extension of GPI in a higher-dimensional space, thereby establishing their optimality. In contrast, this paper employs a novel strategy specifically tailored for the spectral method. Instead of relying on the GPI framework, which proves inadequate for spectral analysis, we introduce a new perturbation toolkit designed for eigenvector analysis. This toolkit provides precise characterization of eigenvector perturbation and leads to the optimality of the spectral method. It opens new avenues for research and application beyond synchronization problems.

In this section, we will extend our analysis to matrix synchronizations where the quantities of interest are orthogonal matrices instead of phases. The orthogonal group synchronization problem has been briefly introduced in Section 1. Here we provide more details about the problem.

Let $d>0$ be an integer. Recall the definition of $\mathcal{O}(d)$ in (11) and that $Z^{*}_{1},\ldots,Z^{*}_{n}\in\mathcal{O}(d)$. For each $1\leq j<k\leq n$, the observation $\mathcal{X}_{jk}\in\mathbb{R}^{d\times d}$ is given by

$\displaystyle\mathcal{X}_{jk}:=\begin{cases}{Z_{j}^{*}Z_{k}^{*{\mathrm{\scriptscriptstyle T}}}+\sigma\mathcal{W}_{jk}},\text{ if }A_{jk}=1,\\ 0,\quad\quad\quad\quad\quad\quad\text{ if }A_{jk}=0,\end{cases}$

(20)

where $A_{jk}\sim\text{Bernoulli}(p)$ and $\mathcal{W}_{jk}\sim\mathcal{MN}(0,I_{d},I_{d})$, i.e., the standard matrix Gaussian distribution¹¹1A random matrix $X$ follows a matrix Gaussian distribution $\mathcal{MN}(M,\Sigma,\Omega)$ if its density function is proportional to $\exp\left(-\frac{1}{2}\text{Tr}\left(\Omega^{-1}(X-M)^{{\mathrm{\scriptscriptstyle T}}}\Sigma^{-1}(X-M)\right)\right)$.. We assume $\{A_{jk}\}_{1\leq j<k\leq n},\{\mathcal{W}_{j,k}\}_{1\leq j<k\leq n}$ are all independent of each other. Similar to the phase synchronization problem, the observations are missing at random with additive Gaussian noises. The goal is to recover $Z_{1}^{*},\ldots,Z_{n}^{*}$ from $\{\mathcal{X}_{jk}\}_{1\leq j<k\leq n}$ and $\{A_{j,k}\}_{1\leq j<k\leq n}$.

The data matrix $\mathcal{X}\in\mathbb{R}^{nd\times nd}$ can be written equivalently in a way that is analogous to (2). Define $A_{jj}:=0$ and $A_{kj}:=A_{jk}$ for all $1\leq j<k\leq n$. Define $\mathcal{W}\in\mathbb{C}^{nd\times nd}$ such that $\mathcal{W}_{jj}:=0_{d\times d}$ and $\mathcal{W}_{kj}:=\mathcal{W}_{jk}^{\mathrm{\scriptscriptstyle T}}$ for all $1\leq j<k\leq n$. Then we have the expression:

$\displaystyle\mathcal{X}=(A\otimes J_{d})\circ(Z^{*}Z^{*{\mathrm{\scriptscriptstyle T}}}+\sigma\mathcal{W})=(A\otimes J_{d})\circ Z^{*}Z^{*{\mathrm{\scriptscriptstyle T}}}+\sigma(A\otimes J_{d})\circ\mathcal{W}.$

(21)

From (12), the data matrix $\mathcal{X}$ can be seen as a noisy version of $pZ^{*}Z^{*{\mathrm{\scriptscriptstyle T}}}$. Since the columns of $Z^{*}$ are orthogonal to each other, we have the following eigendecomposition: $pZ^{*}Z^{*{\mathrm{\scriptscriptstyle T}}}=np(Z^{*}/\sqrt{n})(Z^{*}/\sqrt{n})^{\mathrm{\scriptscriptstyle T}}$ where $Z^{*}/\sqrt{n}\in\mathcal{O}(nd,d)$. That is, $np$ is the only non-zero eigenvalue of $pZ^{*}Z^{*{\mathrm{\scriptscriptstyle T}}}$ with multiplicity $d$.