# Projected Gradient Descent for Spectral Compressed Sensing via Symmetric
Hankel Factorization

Jinsheng Li  Wei Cui  and Xu Zhang  J. Li and W. Cui are with the School of
Information and Electronics, Beijing Institute of Technology, Beijing 100081,
China (e-mail<EMAIL_ADDRESS>cuiwei@bit.edu.cn). X. Zhang is with the
School of Artificial Intelligence, Xidian University, Xi’an 710126, China
(e-mail: zhang.xu@xidian.edu.cn). Corresponding author: Xu Zhang.

###### Abstract

Current spectral compressed sensing methods via Hankel matrix completion
employ symmetric factorization to demonstrate the low-rank property of the
Hankel matrix. However, previous non-convex gradient methods only utilize
asymmetric factorization to achieve spectral compressed sensing. In this
paper, we propose a novel nonconvex projected gradient descent method for
spectral compressed sensing via symmetric factorization named Symmetric Hankel
Projected Gradient Descent (SHGD), which updates only one matrix and avoids a
balancing regularization term. SHGD reduces about half of the computation and
storage costs compared to the prior gradient method based on asymmetric
factorization. Besides, the symmetric factorization employed in our work is
completely novel to the prior low-rank factorization model, introducing a new
factorization ambiguity under complex orthogonal transformation. Novel
distance metrics are designed for our factorization method and a linear
convergence guarantee to the desired signal is established with
$O(r^{2}\log(n))$ observations. Numerical simulations demonstrate the superior
performance of the proposed SHGD method in phase transitions and computation
efficiency compared to state-of-the-art methods.

###### Index Terms:

Spectral compressed sensing, Hankel matrix completion, symmetric matrix
factorization

##  I Introduction

In this paper, we study the spectrally compressed sensing problem, which aims
to recover a spectrally sparse signal from partial measurements. Spectrally
compressed sensing widely arises in applications such as target localization
in radar systems [1], analog-to-digital conversion [2], magnetic resonance
imaging [3], and channel estimation [4]. Denote the spectrally sparse signal
as $\bm{x}=[x_{0},x_{1},\cdots,x_{n-1}]^{T}\in\mathbb{C}^{n}$, whose elements
are a superposition of $r$ complex sinusoids

$\displaystyle x(t)=\sum_{k=1}^{r}d_{k}e^{(i2\pi f_{k}-\tau_{k})t},$ (1)

where $t=0,\cdots,n-1$, $f_{k}\in[0,1)$ is the $k$-th normalized frequency,
$d_{k}$ is the amplitude of the $k$-th sinusoids, and $\tau_{k}$ is the
damping factor.

Only a portion of the spectrally sparse signals can be observed in practice as
obtaining the complete set of sampling points is time-consuming and
technically challenging due to hardware limitations. Additionally, the
canonical full sampling on a uniform grid is inefficient due to the weakening
effects caused by the damping factors $\\{\tau_{k}\\}$. One natural question
is how to recover the entire elements of the desired signal $\bm{x}$ from
partial and nonuniform known entries, i.e.,

$\displaystyle\mbox{Find}\quad\bm{x}\quad\mbox{subject
to}\quad\mathcal{P}_{\Omega}(\bm{x})=\sum_{a\in{\Omega}}x_{a}\bm{e}_{a},$ (2)

where $\bm{e}_{a}$ is the canonical basis of $\mathbb{R}^{n}$,
${\Omega}\subseteq\\{0,\cdots,n-1\\}$ is the sampling index set with
cardinality as $m$, $\mathcal{P}_{\Omega}$ is an operator which projects the
signal $\bm{x}$ onto the sampling index set ${\Omega}$, and
$\mathcal{P}_{\Omega}(\bm{z})=\sum_{a\in{\Omega}}\left\langle\bm{z},\bm{e}_{a}\right\rangle\bm{e}_{a}$
for $\bm{z}\in\mathbb{C}^{n}$.

To reconstruct the spectrally sparse signals from partial observations, many
spectral compressed sensing approaches [5, 6, 7, 8] were proposed to exploit
the spectral sparsity via the low-rankness of the lifted Hankel matrix
$\mathcal{H}\bm{x}$, where
$\mathcal{H}:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n_{1}\times n_{2}}$ is the
Hankel lifting operator, and thus $\mathcal{H}\bm{x}$ is an $n_{1}\times
n_{2}~{}(n=n_{1}+n_{2}-1)$ Hankel matrix. Especially, the low-rankness of the
lifted Hankel matrix $\mathcal{H}\bm{x}$ arises from the following Vandermonde
decomposition

$\mathcal{H}\bm{x}=\bm{E}_{L}\bm{D}\bm{E}_{R}^{T},$

where $\bm{E}_{L}$ is an $n_{1}\times r$ Vandermonde matrix whose $k$-th
column is $[1,w_{k},\cdots,w_{k}^{n_{1}-1}]^{T}$, $\bm{E}_{R}$ is an
$n_{2}\times r$ Vandermonde matrix whose $k$-th column is
$[1,w_{k},\cdots,w_{k}^{n_{2}-1}]^{T}$, $w_{k}=e^{(i2\pi f_{k}-\tau_{k})}$ and
$\bm{D}=\operatorname{diag}(d_{1},\cdots,d_{r})$. When the frequencies
$\\{f_{k}\\}$ are distinct and each diagonal entry $d_{k}$ is nonzero, we have
$\operatorname{rank}(\mathcal{H}\bm{x})=r$.

Using the low rankness of the lifted Hankel matrix, a convex relaxation
approach named Enhanced Matrix Completion (EMaC) was proposed in [5] to
reconstruct spectrally sparse signals. However, the convex approach faces a
computational challenge when dealing with large-scale problems. To address the
computational challenge, several non-convex approaches were proposed,
including PGD [6], FIHT [7], and ADMM with prior information [8], which
demonstrate a lower computation cost than convex approaches [5, 9], especially
in high-dimensional situations. In particular, PGD is a non-convex projected
gradient descent method that applies asymmetric low-rank factorization to
parameterize the lifted low-rank Hankel matrix, i.e.,
$\mathcal{H}\bm{x}=\bm{Z}_{U}\bm{Z}_{V}^{H}$ and introduces a regularization
term to reduce the solution space, where $\bm{Z}_{U}\in\mathbb{C}^{n_{1}\times
r}$ and $\bm{Z}_{V}\in\mathbb{C}^{n_{2}\times r}$.

Figure 1: Left: Asymmetric factorization and update rule of PGD. Right:
Symmetric factorization and update rule of the proposed SHGD. The loss
function, $f(\cdot)$, is defined separately in the PGD and SHGD approaches,
while the regularization function $g(\cdot)$ is specific to PGD. SHGD does not
need a balancing regularization term and updates one single factor at a time
compared to PGD, reducing at least half of operations and storage costs.

###  I-A Motivation and Contributions

Current spectral compressed sensing methods via Hankel matrix completion
employ symmetric factorization to demonstrate the low-rankness of the lifted
Hankel matrix111Without loss of generality, we suppose that the dimension of
the signal is odd; otherwise, we can make it odd by either zero padding after
the end of the signal or deleting the last entry. . However, these methods
employ a non-convex gradient descent algorithm using asymmetric factorization
to solve the spectral compressed sensing problem. This raises the question of
whether a non-convex gradient method based on symmetric factorization can be
used instead.

Fortunately, when the dimension of the signal is odd, the lifted matrix can be
a square Hankel matrix $\mathcal{H}\bm{x}\in\mathbb{C}^{n_{s}\times
n_{s}}~{}(n=2n_{s}-1)$ that enjoys the following complex symmetric
factorization via Vandermonde decomposition

$\displaystyle\mathcal{H}\bm{x}=\bm{E}_{L}\bm{D}\bm{E}_{R}^{T}=\bm{E}\bm{D}\bm{E}^{T}=\bm{Z}_{\natural}\bm{Z}_{\natural}^{T},$

where we can set the $n_{s}\times r$ Vandermonde matrix
$\bm{E}_{L}=\bm{E}_{R}\triangleq\bm{E}$ in square case and the $n_{s}\times r$
matrix $\bm{Z}_{\natural}\triangleq\bm{E}\bm{D}^{\frac{1}{2}}$. This inspires
us to use a low-rank complex symmetric factorization
$\bm{M}=\mathcal{H}\bm{z}=\bm{Z}\bm{Z}^{T}$ with the Hankel constraint to
estimate $\mathcal{H}\bm{x}$ where $\bm{Z}\in\mathbb{C}^{n_{s}\times r}$, and
further recover the desired signal $\bm{x}$. Using symmetric factorization
offers two advantages. First, we can avoid the need for a regularization term
to balance asymmetric matrices. Second, we can reduce computation time and
storage costs by optimizing a single matrix instead of two.

In this paper, we propose a projected gradient descent method under symmetric
Hankel factorization to solve the spectrally compressed sensing problem, named
Symmetric Hankel Projected Gradient Descent (SHGD). The differences between
SHGD and PGD are shown in Fig. 1, where SHGD reduces at least half of
operations and storage costs compared to PGD. Our main contributions are
twofold:

1) A new nonconvex symmetrically factorized gradient descent method named SHGD
is proposed. It is demonstrated that SHGD enjoys a linear convergence rate
towards the desired spectrally sparse signal with high probability, provided
the number of observations is $O(r^{2}\log(n))$. Extensive numerical
simulations are conducted to validate the performance of the proposed SHGD.
Compared to PGD, SHGD reduces about half of the computational cost and storage
requirement. Additionally, SHGD demonstrates superior phase transition
behavior compared to FIHT, while maintaining comparable computational
efficiency.

2) The complex symmetric factorization $\bm{M}=\bm{Z}\bm{Z}^{T}$ employed in
our work is novel to the prior low-rank factorization model in the literature
222There are two types of factorization for low-rank real matrices:
$\bm{M}=\bm{L}\bm{R}^{T}$ for asymmetric matrices and
$\bm{M}=\bm{Z}\bm{Z}^{T}$ for positive semi-definite matrices. However, there
are three types of factorization for complex matrices:
$\bm{M}=\bm{L}\bm{R}^{H}$ for asymmetric matrices, $\bm{M}=\bm{Z}\bm{Z}^{H}$
for Hermitian positive semi-definite matrices, and $\bm{M}=\bm{Z}\bm{Z}^{T}$
for complex symmetric matrices, which is first pointed out in our work..
Complex symmetric factorization introduces a new factorization ambiguity
$\bm{M}=\bm{Z}\bm{Z}^{T}=\bm{Z}\bm{Q}(\bm{Z}\bm{Q})^{T}$ where $\bm{Q}$ is
called complex orthogonal matrix (not unitary) [10, 11, 12, 13] such that
$\bm{Q}\bm{Q}^{T}=\bm{I}$ for $\bm{Q}\in\mathbb{C}^{r\times r}$. Such
factorization ambiguity under complex orthogonal transformation is first
pointed out in our work. Consequently, the analyses are non-trivial and highly
technical to deal with a new factorization model and algorithm. Novel distance
metrics are designed to analyze such factorization model.

###  I-B Related work

When a spectral sparse signal is undamped and its frequencies are discretized
on a uniform grid, one can use conventional compressed sensing [14, 15] to
estimate its spectrum. However, the true frequencies are usually continuous
and off-the-grid in practice, thus mismatch errors arise when applying
conventional compressed sensing approaches. A grid-free approach named atomic
norm minimization (ANM) was proposed by exploiting the sparsity of the
frequency domain in a continuous way [9]. When frequencies are well-separated,
ANM can achieve exact recovery with high probability from $O(r\log(r)\log(n))$
random samples. Another grid-free approach was proposed to exploit the
sparsity of frequencies via the low-rank Hankel completion model [5]. By
formulating it as a nuclear norm minimization problem, the authors showed
$O(r\log^{4}(n))$ random observations are sufficient to recover the true
signals exactly with high probability, equipped with the incoherence
conditions. Both of them are convex approaches and suffer from computation
issues for large-scale problems.

Several provable non-convex approaches were proposed to solve the previous
low-rank Hankel matrix completion problem towards spectrally compressed
sensing, such as PGD [6], FIHT [7] and pMAP [16]. These methods offer lower
computational costs compared to convex approaches, particularly in high-
dimensional scenarios. Inspired by a factorized projected gradient method
towards matrix completion [17], a similar gradient method PGD [6] that employs
low-rank factorization was proposed for Hankel matrix completion. Equipped
with $O(r^{2}\log(n))$ random samples, PGD enjoys an exact recovery with high
probability [6]. However, the inherent symmetric structure of Hankel matrices
wasn’t exploited thoroughly in [6]. The complex symmetry of the square Hankel
matrix was exploited in [18] for fast symmetric SVD of Hankel matrices and in
[19] for a fast alternating projection method towards complex frequency
estimation, but there are no theoretical guarantees in [18] and [19]. FIHT [7]
is an algorithm that employs the Riemannian optimization technique towards
low-rank matrices [20, 21] to solve the Hankel matrix completion problem, and
the complex symmetry could be exploited for further reduction of storage and
computation. FIHT achieves an exact recovery with high probability equipped
with $O(r^{2}\log^{2}(n))$ random samples, which is worse with an order of
$\log(n)$ compared to PGD and SHGD. A penalized version of the Method of
Alternating Projections (MAP) was proposed in [16] towards weighted low-rank
Hankel optimization, of which the low-rank Hankel matrix completion was a
special case. pMAP enjoys a linear convergence rate to the truth with a
reasonable computation cost. Considering the prior information in the spectral
compressed sensing problem, the authors in [8] employed low-rank
factorization, and proposed an ADMM framework with prior information, but
there is no convergence rate analysis for their nonconvex approach. A
concurrent work called HT-GD [22] was proposed to employ the symmetric
property of low-rank Hankel matrices and Hermitian property of low-rank
Toeplitz matrices, but it focuses on undamped signals and lacks theoretical
guarantees [22]. Inspired by the previous PGD method [6], the authors in [23]
proposed a projected gradient descent method via vectorized Hankel lift (PGD-
VHL) towards blind super-resolution. They also apply the approach which is
based on low-rank asymmetric factorization with a balancing regularization
term.

Symmetric factorization arises in a variety of fields, such as non-negative
matrix factorization for clustering [24, 25], matrix completion [26], mobile
communications, and factor analysis where symmetric tensor factorization is
employed [27, 28, 29], further reducing computation and storage costs.
However, the complex symmetric factorization in this work is a new type of
factorization, which hasn’t appeared in the literature to our knowledge,
enriching the range of low-rank matrix factorization approaches.

Notations. We denote vectors with bold lowercase letters, matrices with bold
uppercase letters, and operators with calligraphic letters. For matrix
$\bm{Z}$, we use $\bm{Z}^{T}$, $\bm{Z}^{H}$, $\bar{\bm{Z}}$,
$\left\lVert\bm{Z}\right\rVert$, and $\left\lVert\bm{Z}\right\rVert_{F}$ to
denote its transpose, conjugate transpose, complex conjugate, spectral norm,
and Frobenius norm, respectively. Besides, we define
$\left\lVert\bm{Z}\right\rVert_{2,\infty}$ as the largest $\ell_{2}$-norm of
its rows. For vector $\bm{z}$, $\left\lVert\bm{z}\right\rVert_{1}$ and
$\left\lVert\bm{z}\right\rVert_{2}$ denotes its $\ell_{1}$ and $\ell_{2}$
norms. Define the inner product of two matrices $\bm{Z}_{1}$ and $\bm{Z}_{2}$
as
$\left\langle\bm{Z}_{1},\bm{Z}_{2}\right\rangle=\mathrm{trace}(\bm{Z}_{1}^{H}\bm{Z}_{2})$.
$[n]$ denotes the set $\\{0,1,\cdots,n-1\\}$ where $n$ is a natural number. We
denote the identity matrix and operator as $\bm{I}$ and $\mathcal{I}$,
respectively. The adjoint of the operator $\mathcal{A}$ is denoted as
$\mathcal{A}^{*}$. $\operatorname{{Re}}(\cdot)$ denotes the real part of a
complex number.

##  II Problem formulation

Without loss of generality, we suppose the length $n$ of the desired signal
$\bm{x}$ is odd in the following. We aim to recover the desired signal via
symmetric Hankel matrix completion. Firstly, we form a $n_{s}\times n_{s}$
complex symmetric Hankel matrix

$\displaystyle\mathcal{H}\bm{x}=\begin{bmatrix}x_{0}&x_{1}&\cdots&x_{n_{s}-1}\\\
x_{1}&x_{2}&\cdots&x_{n_{s}}\\\ \vdots&\vdots&\vdots&\vdots\\\
x_{n_{s}-1}&x_{n_{s}}&\cdots&x_{n-1}\end{bmatrix},$

where $\mathcal{H}:\mathbb{C}^{n}\rightarrow\mathbb{C}^{n_{s}\times
n_{s}}~{}(n=2n_{s}-1)$ is the Hankel lifting linear operator. The lifted
symmetric Hankel matrix $\mathcal{H}\bm{x}$ has a Vandermonde decomposition as
pointed out in the introduction:

$\mathcal{H}\bm{x}=\bm{E}\bm{D}\bm{E}^{T},$

where $\bm{E}$ is an $n_{s}\times r$ Vandermonde matrix with its $k$-th column
as $[1,w_{k},\cdots,w_{k}^{n_{s}-1}]^{T}$ and $w_{k}=e^{(i2\pi
f_{k}-\tau_{k})}$. When frequencies are distinct and each diagonal entry
$d_{k}$ of $\bm{D}$ is non-zero, $\operatorname{rank}(\mathcal{H}\bm{x})=r$.

We denote $\mathcal{H}^{*}$ as the adjoint of $\mathcal{H}$, which is to map a
matrix $\bm{M}\in\mathbb{C}^{n_{s}\times n_{s}}$ to a vector
$\mathcal{H}^{*}\bm{M}=\\{\sum_{i+j=a}\bm{M}_{(i,j)}\\}_{a=0}^{n-1}$. Let
$\mathcal{D}^{2}=\mathcal{H}^{*}\mathcal{H}$ and $\mathcal{D}$ is a linear
operator that maps a vector $\bm{x}\in\mathbb{C}^{n}$ to
$\mathcal{D}\bm{x}\in\mathbb{C}^{n}$, where
$[\mathcal{D}\bm{x}]_{a}=\sqrt{w_{a}}x_{a}$ and $w_{a}$ is the length of the
$a$-th skew-diagonal of an $n_{s}\times n_{s}$ matrix with $a\in[n-1]$.
Defining $\mathcal{G}=\mathcal{H}\mathcal{D}^{-1}$ as $\mathcal{D}$ is
invertible, then we have $\mathcal{G}^{*}\mathcal{G}=\mathcal{I}$.

As stated in [6], a rank constraint weighted least square problem is
constructed to recover $\bm{x}$, utilizing the exact low-rankness property of
$\mathcal{H}\bm{x}$:

$\displaystyle\min_{\bm{z}\in\mathbb{C}^{n}}\left\langle\mathcal{P}_{\Omega}\left(\mathcal{D}(\bm{z}-\bm{x})\right),\mathcal{D}(\bm{z}-\bm{x})\right\rangle~{}~{}\mbox{s.t.}~{}\operatorname{rank}(\mathcal{H}\bm{z})=r.$
(3)

The problem (3) can be reformulated as the following optimization problem
after making the substitutions that $\bm{y}\leftarrow\mathcal{D}\bm{x}$ and
$\bm{z}\leftarrow\mathcal{D}\bm{z}$

$\displaystyle\min_{\bm{z}\in\mathbb{C}^{n}}\left\langle\mathcal{P}_{\Omega}\left(\bm{z}-\bm{y}\right),\bm{z}-\bm{y}\right\rangle~{}\mbox{s.t.}~{}\operatorname{rank}(\mathcal{G}\bm{z})=r.$
(4)

Note that $\mathcal{G}\bm{z}\in\mathbb{C}^{n_{s}\times n_{s}}$ is a rank-$r$
square Hankel matrix and consequently complex symmetric.

For any complex symmetric matrix $\bm{M}=\bm{M}^{T}\in\mathbb{C}^{n_{s}\times
n_{s}}$, there is a special form of SVD called Takagi factorization [10],
which is defined as $\bm{M}=\bm{U}\bm{\Sigma}\bm{U}^{T},$ where $\bm{U}$ is a
unitary matrix and $\bm{\Sigma}$ is a diagonal singular value matrix. If we
set $\bm{Z}=\bm{U}\bm{\Sigma}^{\frac{1}{2}}$, then the complex symmetric
matrix $\bm{M}$ enjoys symmetric factorization $\bm{M}=\bm{Z}\bm{Z}^{T}$.

Therefore, we can employ low-rank symmetric factorization
$\mathcal{G}\bm{z}=\bm{Z}\bm{Z}^{T}$ to eliminate the rank constraint, where
$\bm{Z}\in\mathbb{C}^{n_{s}\times r}$, and enforce Hankel structure by the
following Hankel constraint

$\displaystyle(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}\bm{Z}^{T})=\bm{0},$

where $\mathcal{G}\mathcal{G}^{*}$ is a projector that maps a matrix to a
Hankel matrix. Thus (4) can be rewritten as

$\displaystyle\min_{\bm{Z}\in\mathbb{C}^{n_{s}\times
r}}\left\langle\mathcal{P}_{\Omega}\left(\mathcal{G}^{*}\left(\bm{Z}\bm{Z}^{T}\right)-\bm{y}\right),\mathcal{G}^{*}\left(\bm{Z}\bm{Z}^{T}\right)-\bm{y}\right\rangle$
$\displaystyle\quad\mbox{s.t.}\quad(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}\bm{Z}^{T})=\bm{0},$
(5)

where we use the fact that
$\bm{z}=\mathcal{G}^{*}\left(\mathcal{G}\bm{z}\right)=\mathcal{G}^{*}\left(\bm{Z}\bm{Z}^{T}\right)$,
and substitute this in (4).

Furthermore, we consider a penalized version of (5) to recover the weighted
desired signal $\bm{y}$,

$\displaystyle f(\bm{Z})=$
$\displaystyle\frac{1}{4p}\left\langle\mathcal{P}_{\Omega}(\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y}),\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y}\right\rangle$
$\displaystyle+\frac{1}{4}\left\|(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}\bm{Z}^{T})\right\|_{F}^{2},$
(6)

where $p=m/n$ is the sampling ratio. We can interpret (6) as that one uses a
low-rank symmetrically factorized matrix $\bm{M}=\bm{Z}\bm{Z}^{T}$ with Hankel
structure (Hankel structure enforced in (6) by a penalty term) to estimate the
lifted matrix $\mathcal{G}\bm{y}$ via minimizing the mismatch in the
measurement domain.

##  III Definitions and algorithms

###  III-A Definitions

Denoting the lifted truth matrix as $\bm{M}_{\star}=\mathcal{G}\bm{y}$, which
is a complex symmetric matrix as pointed out in the previous section. Let the
Takagi factorization of $\bm{M}_{\star}$ being
$\bm{M}_{\star}=\bm{U}_{\star}\bm{\Sigma}_{\star}\bm{U}_{\star}^{T}$ and set

$\displaystyle\bm{Z}_{\star}=\bm{U}_{\star}\bm{\Sigma}_{\star}^{\frac{1}{2}},$
(7)

then $\bm{M}_{\star}=\bm{Z}_{\star}\bm{Z}_{\star}^{T}$.

If the singular vectors of $\bm{M}_{\star}$ are aligned with the sampling
basis, we can’t recover $\bm{M}_{\star}$ from element-wise observations.
Consequently, we suppose that $\bm{M}_{\star}=\mathcal{G}\bm{y}$ is
$\mu_{0}$-incoherent as in [6, 30],

######  Definition 1.

The $n_{s}\times n_{s}$ square Hankel matrix
$\bm{M}_{\star}=\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent, if there exists an
absolute numerical constant $\mu_{0}>0$ such that

$\displaystyle\left\|\bm{U}_{\star}\right\|_{2,\infty}\leq\sqrt{\frac{2\mu_{0}r}{n}},$

where $n=2n_{s}-1$, $\bm{U}_{\star}\bm{\Sigma}_{\star}\bm{U}_{\star}^{T}$ is
the Takagi factorization of $\bm{M}_{\star}$.

######  Remark 1.

We remove the dependence on $c_{s}$ for the incoherence condition defined in
[5, 6, 7], as $c_{s}=n/n_{s}\leq 2$ for the square case and only suppose the
incoherence condition for one single matrix $\bm{U}_{\star}$.

It has been demonstrated in [6] and [31, Thm. 2] that
$\bm{M}_{\star}=\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent as long as the
minimum wrap-around distance between the frequencies is greater than about
$2/n$, and the spectrally sparse signals are undamped.

Let $\mu$ and $\sigma$ be numerical constants such that $\mu\geq\mu_{0}$ and
$\sigma\geq\sigma_{1}(\bm{M}_{\star})$. When $\bm{M}_{\star}$ is
$\mu_{0}$-incoherent, the matrix $\bm{Z}_{\star}$ satisfies
$\left\|\bm{Z}_{\star}\right\|_{2,\infty}\leq\sqrt{2\mu r\sigma/n}$. Moreover,
let $\mathcal{C}$ be a convex set defined as

$\displaystyle\mathcal{C}=\left\\{\bm{Z}\in\mathbb{C}^{n_{s}\times
r}~{}|~{}\left\|\bm{Z}\right\|_{2,\infty}\leq 2\sqrt{\frac{\mu
r\sigma}{n}}\right\\}.$ (8)

We immediately have $\bm{Z}_{\star}\in\mathcal{C}$ and it is reasonable to
project the factor $\bm{Z}$ onto this set. We define the following projection
operator $\mathcal{P}_{\mathcal{C}}(\cdot)$, for
$\bm{Z}\in\mathbb{C}^{n_{s}\times r}$

$\displaystyle[\mathcal{P}_{\mathcal{C}}(\bm{Z})]^{(i,:)}=\begin{cases}\bm{Z}^{(i,:)}&\mbox{if
}\|\bm{Z}^{(i,:)}\|_{2}\leq 2\sqrt{\frac{\mu r\sigma}{n}},\\\
\frac{\bm{Z}^{(i,:)}}{\|\bm{Z}^{(i,:)}\|_{2}}2\sqrt{\frac{\mu
r\sigma}{n}}&\mbox{otherwise}.\end{cases}$

###  III-B Algorithms

As discussed above, we construct the following constrained optimization
problem under symmetric factorization

$\displaystyle\min_{\bm{Z}\in\mathcal{C}}~{}f(\bm{Z})=$
$\displaystyle\frac{1}{4p}\left\langle\mathcal{P}_{\Omega}(\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y}),\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y}\right\rangle$
$\displaystyle+\frac{1}{4}\left\|(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}\bm{Z}^{T})\right\|_{F}^{2}.$
(9)

Under Wirtinger calculus, the gradient of the loss function $f(\bm{Z})$ is

$\displaystyle\nabla
f(\bm{Z})=p^{-1}\left(\mathcal{G}\mathcal{P}_{\Omega}(\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y})\right)\bar{\bm{Z}}+(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}\bm{Z}^{T})\bar{\bm{Z}}.$
(10)

Note that the gradient of our loss function avoids the computation and storage
for two factors as well as a balancing regularization term, compared to PGD
[6].

Algorithm 1 Symmetric Hankel Projected Gradient Descent

Preprocessing: 1\. Make each dimension of the signal odd by zero-padding.

2\. Partition ${\Omega}$ into disjoint sets ${\Omega}_{0},\cdots,{\Omega}_{K}$
of equal size $\hat{m}$, let $\hat{p}=\frac{\hat{m}}{n}$.

Initialization:
$\bm{M}^{0}=\mathcal{T}_{r}\left(\hat{p}^{-1}\mathcal{G}\mathcal{P}_{{\Omega}_{0}}(\bm{y})\right)=\bm{U}^{0}\bm{\Sigma}^{0}(\bm{U}^{0})^{T}$,
$\tilde{\bm{Z}}^{0}=\bm{U}^{0}(\bm{\Sigma}^{0})^{1/2}$ and
$\bm{Z}^{0}=\mathcal{P}_{\mathcal{C}}(\tilde{\bm{Z}}^{0})$.

for $k=0,1,\cdots,K$ do

1\. ${\tilde{\bm{Z}}^{k+1}}=\bm{Z}^{k}-\eta\nabla f^{(k)}(\bm{Z}^{k})$

2\. $\bm{Z}^{k+1}=\mathcal{P}_{\mathcal{C}}(\tilde{\bm{Z}}^{k+1})$

end for

Output: $\bm{Z}^{K}$ in the last iteration,
$\bm{y}^{K}=\mathcal{G}^{*}(\bm{Z}^{K}(\bm{Z}^{K})^{T})$ and
$\bm{x}^{K}=\mathcal{D}^{-1}\bm{y}^{K}$.

We design a projected gradient descent algorithm with sample-splitting for
solving (9), named as Symmetric Hankel Projected Gradient Descent, seeing
Algorithm 1. In preprocessing, we first transform the dimension of the
observed signal odd by zero-padding if not. Secondly, the sampling set
${\Omega}$ is partitioned into $K+1$ disjoint subsets of equal size $\hat{m}$.
Such partitions of the sample set are commonly used in analyzing matrix
completion problems [32, 33] and Hankel matrix completion problems [7, 34,
35]. The sample-splitting technique keeps the independence between the current
sampling set and the previous iterates, which simplifies the theoretical
analyses.

The following is the initialization, which is one-step hard thresholding after
Takagi factorization and then projecting the initial factor to the convex set
$\mathcal{C}$. There are several algorithms to implement Takagi factorization
[18, 36, 37] and one may use the algorithm in [18] directly.

Finally, we enter the stage that iteratively updates on single factor
$\bm{Z}$. Considering the sample splitting setting, the gradient term $\nabla
f^{(k)}(\bm{Z}^{k})$ in the $k$-th iteration is

$\displaystyle\nabla f^{(k)}(\bm{Z}^{k})$
$\displaystyle=\hat{p}^{-1}\left(\mathcal{G}\mathcal{P}_{{\Omega}_{k}}(\mathcal{G}^{*}(\bm{Z}^{k}({\bm{Z}^{k}})^{T})-\bm{y})\right)\bar{\bm{Z}}^{k}$
$\displaystyle\quad+(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}^{k}({\bm{Z}^{k}})^{T})\bar{\bm{Z}}^{k}.$

Our model and algorithm can be generalized to multi-dimensional spectral
sparse signals. We omit the details due to space limitations.

###  III-C Computational complexity for SHGD and PGD

The empirical computational efficiency improvement of SHGD relies on specific
implementations and characteristics of the problem being considered, such as
the problem scale $n$ and the number of sinusoids $r$. A detailed comparison
of the empirical computational time between SHGD and PGD is shown in Section
V-B. In this part, we provide an analysis of computational complexity for SHGD
and PGD as follows.

The gradient of SHGD (10) can be reformulated as $\nabla
f(\bm{Z})=\mathcal{G}\left(p^{-1}\mathcal{P}_{\Omega}(\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y})-\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})\right)\bar{\bm{Z}}+\bm{Z}(\bm{Z}^{T}\bar{\bm{Z}})$.
We compute $\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})$ by $r$ fast convolutions.
Setting
$\bm{w}=p^{-1}\mathcal{P}_{\Omega}(\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})-\bm{y})-\mathcal{G}^{*}(\bm{Z}\bm{Z}^{T})$,
$\left(\mathcal{G}\bm{w}\right)\bar{\bm{Z}}$ can be computed by $r$ fast
Hankel matrix-vector multiplications. Each of the aforementioned two steps
needs $Cnr\log(n)$ flops, where $C>0$ is a constant333The constant factors may
differ for the two steps, but we can choose the larger one as $C$.. Besides,
the computation of $\bm{Z}(\bm{Z}^{T}\bar{\bm{Z}})$ requires
$2n_{s}r^{2}(\approx nr^{2})$ flops. For PGD [6], there are three steps of
convolution type, each requiring $Cnr\log(n)$ flops, and four steps associated
with balancing regularization, each using about $nr^{2}$ flops. Thus the
computation time ratio between SHGD and PGD per-iteration is
$\frac{2Cnr\log(n)+nr^{2}}{3Cnr\log(n)+4nr^{2}}=\frac{2C\log(n)+r}{3C\log(n)+4r}$.
Depending on the relative scale between $r$ and $C\log(n)$, this ratio ranges
from $1/4$ to $2/3$.

##  IV Theoretical results and Analysis

In this section, we first discuss the challenges associated with analyzing the
proposed SHGD algorithm. We then present the linear convergence result of
SHGD. Lastly, we introduce our analysis framework and provide the proof of
Theorem 1.

###  IV-A Theoretical challenges

The previous analysis for PGD primarily relies on asymmetric factorization and
balancing regularization [6], which makes it feasible to design a distance
metric to effectively address unitary ambiguity. However, extending the
analysis from PGD to SHGD is challenging. The complex symmetric factorization
is novel to the prior low-rank factorization model, introducing a new
factorization ambiguity
$\bm{M}_{\star}=\bm{Z}_{\star}\bm{Z}_{\star}^{T}=\bm{Z}_{\star}\bm{Q}(\bm{Z}_{\star}\bm{Q})^{T}$
for $\bm{M}_{\star}\in\mathbb{C}^{n_{s}\times n_{s}}$, where $\bm{Q}$ is
called complex orthogonal matrix (not unitary unless real) [10, 11, 12] such
that $\bm{Q}\in\mathcal{Q}\triangleq\\{\bm{S}\in\mathbb{C}^{r\times
r}|\bm{S}\bm{S}^{T}=\bm{S}^{T}\bm{S}=\bm{I}\\}$. Such factorization ambiguity
under complex orthogonal transformation is first pointed out in our work. The
analysis of complex orthogonal ambiguity is more challenging than that of the
unitary ambiguity, since $\|\bm{Z}\bm{R}\|=\|\bm{Z}\|$ holds for all unitary
matrix $\bm{R}$ while $\|\bm{Z}\bm{Q}\|\neq\|\bm{Z}\|$ can occur for some
complex orthogonal matrix $\bm{Q}$ as shown in Example 1. To characterize the
distance from any matrix $\bm{Z}$ to the set $\bm{Z}_{\star}\bm{Q}$ as
indicated in [6, 17, 38, 39], it is essential to define a distance metric

$\displaystyle{\mbox{dist}_{Q}}(\bm{Z},\bm{Z}_{\star})=\inf_{\bm{Q}\in\mathcal{Q}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\right\|_{F}.$
(11)

However, it is hard to find the closed-form solution or a concise first-order
optimality condition for (11). Besides, the set of the complex orthogonal
matrix $\bm{Q}$ may be unbounded and therefore non-compact, bringing existence
issues to the solution of the previous definition. See examples as follows.

######  Example 1.

Define hyperbolic functions $\cosh{t}=(e^{t}+e^{-t})/2$ and
$\sinh{t}=(e^{t}-e^{-t})/2$ . Let $\bm{I}_{2}$ be a $2\times 2$ identity
matrix and $\bm{S}=\left[\begin{smallmatrix}0&1\\\
-1&0\end{smallmatrix}\right]$, then
$\bm{Q}(t)=(\cosh{t})\bm{I}_{2}+(j\sinh{t})\bm{S}$ is a complex orthogonal
matrix for all $t\in\mathbb{R}$ and the set
$\bm{Q}=\\{\bm{Q}(t)|t\in\mathbb{R}\\}$ is non-compact.

###### Proof.

$\displaystyle\bm{Q}(t)\bm{Q}(t)^{T}=\bm{Q}(t)^{T}\bm{Q}(t)=\left((\cosh{t})^{2}-(\sinh{t})^{2}\right)\bm{I}_{2}=\bm{I}_{2},$

thus $\bm{Q}(t)$ is a complex orthogonal matrix. Besides, it is obvious that
$\left\|\bm{Q}(t)\right\|$ and $\left\|\bm{Q}(t)\right\|_{F}$
$\rightarrow\infty$ as $t\rightarrow\infty$. So the set $\bm{Q}$ is not
bounded and therefore non-compact. ∎

Notice that $\mbox{dist}^{2}_{Q}(\bm{Z},\bm{Z}_{\star})$ can be reformulated
as

$\displaystyle 2\mbox{dist}^{2}_{Q}(\bm{Z},\bm{Z}_{\star})$
$\displaystyle=\inf_{\begin{subarray}{c}\bm{Q}^{T}\bm{Q}=\bm{Q}\bm{Q}^{T}=\bm{I},\\\
\bm{Q}\in\mathbb{C}^{r\times
r}\end{subarray}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\right\|_{F}^{2}+\|\bm{Z}-\bm{Z}_{\star}\bm{Q}^{-T}\|_{F}^{2},$

$\bm{Q}$ must be invertible as $\bm{Q}^{T}\bm{Q}=\bm{I}$ for
$\bm{Q}\in{\mathbb{C}^{r\times r}}$. We drop the complex orthogonal constraint
$\bm{Q}^{T}\bm{Q}=\bm{I}$ but keep $\bm{Q}$ invertible, and thus define a
novel distance metric444Rigorously, we should use $\inf$ instead of $\min$ for
the distance definition. However, we will prove the existence of a minimizer
in Lemma 4 under some conditions that the iterates of our algorithm satisfy.

$\displaystyle{\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star}){=}\min_{\begin{subarray}{c}\bm{P}\in\mathbb{C}^{r\times
r}\\\ \text{ invertible
}\end{subarray}}\sqrt{\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}\right\|_{F}^{2}+\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}^{-T}\right\|_{F}^{2}}.$
(12)

Additionally, we reveal the relationship between the previous distance metrics
through the following Lemma.

######  Lemma 1.

The distance metrics ${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ and
${\mbox{dist}_{Q}}(\bm{Z},\bm{Z}_{\star})$ satisfy the following
relationships:

1)
${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})\leq\sqrt{2}{\mbox{dist}_{Q}}(\bm{Z},\bm{Z}_{\star})$.

2) Suppose
${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})\leq\varepsilon\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}$,
then ${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ asymptotically reduces to
$\sqrt{2}{\mbox{dist}_{Q}}(\bm{Z},\bm{Z}_{\star})$ when
$\varepsilon\rightarrow 0$.

###### Proof.

See Appendix D-C. ∎

Lemma 1 tells us that ${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ approximates
${\mbox{dist}_{Q}}(\bm{Z},\bm{Z}_{\star})$ well under some conditions. This
inspires us to study ${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ and build the
corresponding convergence. The distance metric
${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ has a nice first-order optimality
condition:

$\displaystyle(\bm{Z}_{\star}\bm{P}_{\bm{Z}})^{H}({\bm{Z}-\bm{Z}_{\star}\bm{P}_{\bm{Z}}})=(\bm{Z}-\bm{Z}_{\star}\bm{P}_{\bm{Z}}^{-T})^{T}\big{(}\overline{\bm{Z}_{\star}\bm{P}_{\bm{Z}}^{-T}}\big{)},$
(13)

where we define the optimal solution $\bm{P}_{\bm{Z}}$ for (12), if exists,
as:

$\displaystyle\bm{P}_{\bm{Z}}:=\underset{\begin{subarray}{c}\bm{P}\in\mathbb{C}^{r\times
r}\\\ \text{ invertible
}\end{subarray}}{\operatorname{argmin}}\sqrt{\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}\right\|_{F}^{2}+\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}^{-T}\right\|_{F}^{2}}.$
(14)

Establishing the convergence mechanism for SHGD in terms of the distance
metric ${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ reveals a distinct deviation
from that for PGD [6]. Unlike PGD, SHGD requires an analysis towards a
different distance metric ${\mbox{dist}_{P}}(\bm{Z},\bm{Z}_{\star})$ and a
different gradient direction under symmetric factorization. Also, the
conditions about the local basin of attraction towards SHGD are characterized
differently from PGD, which are small deviations from the true solution in
terms of our distance metric and well-conditioness of
$\bm{Z}_{\star}\bm{P}_{\bm{Z}}$, $\bm{Z}_{\star}\bm{P}_{\bm{Z}}^{-T}$. At
last, we establish an inductive convergence analysis framework to verify that
the iterates in Algorithm 1 satisfy such conditions and thus enjoy linear
convergence to the ground truth.

###  IV-B Main results

We consider the sampling with replacement model as in [6, 7, 34, 35]. To
distinguish from the sampling index set $\Omega$ and $\Omega_{k}$, we denote
${\hat{\Omega}}=\\{a_{j}|j=1,\cdots,m\\}$ and
${\hat{\Omega}}_{k}=\\{a_{j}|j=1,\cdots,\hat{m}\\}$ as the corresponding
sampling set with replacement, where the indices $a_{j}$ are drawn
independently and uniformly from $\\{0,\cdots,n-1\\}$.

######  Theorem 1.

Assume $\bm{M}_{\star}=\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent. Let
$\eta=\frac{\eta^{\prime}}{\sigma_{1}(\bm{M}_{\star})}$, where
$0<\eta^{\prime}\leq\frac{1}{54}$. Set
$\sigma=\sigma_{1}(\bm{M}_{0})/(1-\varepsilon_{0})$ and $\mu\geq\mu_{0}$,
where $\varepsilon_{0}$ is a small enough constant. If the number of
observations satisfies $m\gtrsim
O\left(\varepsilon_{0}^{-2}\mu^{2}\kappa^{4}r^{2}\log(n)\right)$, the iterates
$\bm{Z}^{k}$ of Algorithm 1 satisfy

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})\leq\left(1-\frac{11\eta^{\prime}}{100\kappa}\right)^{k}c_{1}\hskip
1.42271pt\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}},$

with probability at least $1-O(n^{-2})$, where $c_{1}>0$ is a universal
constant, and $\kappa=\sigma_{1}(\bm{M}_{\star})/\sigma_{r}(\bm{M}_{\star})$.

######  Remark 2.

By establishing the relationship between $\|\bm{x}^{k}-\bm{x}\|_{2}$ and
${\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})\leq\varepsilon$, we can obtain
the $\varepsilon$ recovery accuracy for $\|\bm{x}^{k}-\bm{x}\|_{2}$, i.e.,
$\|\bm{x}^{k}-\bm{x}\|_{2}\leq O(\varepsilon)$, when
${\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})\leq\varepsilon$. See the proof
in Appendix D-B.

######  Remark 3.

Our step size is on the order of $O(\frac{1}{\sigma_{1}(\bm{M}_{\star})})$
independent of $r$ and $\mu$, compared to the conservative stepsize
$O(\frac{\sigma_{r}(\bm{M}_{\star})}{\mu^{2}c_{s}^{2}r^{2}\sigma_{1}^{2}(\bm{M}_{\star})})$
of PGD [6]. To attach $\epsilon$ recovery accuracy, the iteration complexity
is $O(\kappa\log(\frac{1}{\epsilon}))$, independent of $r$ and $\mu$ also,
compared to PGD [6].

######  Remark 4.

With a more careful sample-splitting strategy555One can use a non-uniform
sample-splitting strategy, which is to set the size of the sample set for the
initialization satisfy $\tilde{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu\kappa^{4}r^{2}\log(n)\right)$ from Lemma 2, and
$\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu^{2}\kappa^{2}r^{2}\log(n)\right)$ for the
following $\kappa\log(1/\varepsilon)$ steps from Lemma 3., the sample
complexity of SHGD is
$O(\varepsilon_{0}^{-2}\mu^{2}\kappa^{4}r^{2}\log(n)\log(1/\varepsilon))$ to
attach $\varepsilon$ accuracy in terms of our distance metric. We believe that
the requirement for sample-splitting is just an artifact of our proof and can
be removed safely, which we leave as a potential direction for future work.

######  Remark 5.

By combining the techniques from Theorem 1 and [26], we can generalize the
analysis to noisy measurements
$\bm{y}_{e}=\mathcal{P}_{\Omega}(\bm{x}+\bm{e})$ and give an order-wise
recovery bound for simplicity. Suppose the elements of the noise vector
$\bm{e}$ are independent zero-mean sub-Guassian random variables with
parameter $\sigma$ [40]. With high probability, one has the following robust
recovery guarantee:

$\displaystyle\|\bm{x}-\bm{x}^{k}\|_{2}\lesssim\|\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}\|_{F}\lesssim\sigma\sqrt{\frac{n^{2}}{m}},$

as long as $\sigma\ll\sigma_{r}(\bm{M}_{\star})$ and $m\geq O(r^{2}\log(n))$.

###  IV-C Analytical framework

First of all, we introduce the existence of the minimizer of the problem (12)
under some conditions; see Lemma 4 in Appendix A for the details. One can
easily verify the existence of the minimizer of the distance metric for the
iterates in our Algorithm 1 inductively by invoking Lemma 4. Consequently, we
omit the debation about such existence issues for simplicity of presentation.

We apply a two-stage analysis, combined with an inductive framework. Firstly
we show some good properties of the initialization, which are about the local
basin of attraction specified for SHGD. Then we show that the iterates
returned by Algorithm 1 can converge linearly to the true solution provided
they are in such a local basin of attraction.

######  Lemma 2 (Initialization conditions).

Suppose $\bm{M}_{\star}=\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent,
$\sigma=\sigma_{1}(\bm{M}_{0})/(1-\varepsilon_{0})$, $\mu\geq\mu_{0}$ and
$\varepsilon_{0}$ is a small enough constant. When $\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu\kappa^{4}r^{2}\log(n)\right)$, one has the
following results with probability at least $1-n^{-2}$,

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{0},\bm{Z}_{\star})$ $\displaystyle\leq
1.6\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}};$ (15)
$\displaystyle\sigma_{r}\left(\bm{Z}_{\star{i}}^{0}\right)$
$\displaystyle\geq\left(1-6.4\kappa^{-1}\varepsilon_{0}\right)\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},i=1,2,$
(16) $\displaystyle\sigma_{1}\left(\bm{Z}_{\star{i}}^{0}\right)$
$\displaystyle\leq\left(1+6.4\kappa^{-1}\varepsilon_{0}\right)\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},i=1,2;$
(17) $\displaystyle\sigma$
$\displaystyle=\frac{\sigma_{1}(\bm{M}^{0})}{1-\varepsilon_{0}}\geq\sigma_{1}(\bm{M}_{\star}).$
(18)

Besides, the optimal solution $\bm{P}^{0}$ to
${\mbox{dist}_{P}}(\bm{Z}^{0},\bm{Z}_{\star})$ exists and we set
$\bm{Z}^{0}_{\star{1}}=\bm{Z}_{\star}\bm{P}^{0},\bm{Z}^{0}_{\star{2}}=\bm{Z}_{\star}(\bm{P}^{0})^{-T}$.

###### Proof.

See Appendix B. ∎

Next, we suppose the iterate in the $k$-th step satisfies the following
conditions:

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})\leq
1.6\kappa^{-1}\varepsilon_{0}\left(1-\frac{11\eta^{\prime}}{100\kappa}\right)^{k}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}};$
(19)
$\displaystyle\sigma_{r}(\bm{Z}_{\star{i}}^{k})\geq\Big{[}1-6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{k}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\Big{]}\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},$
(20)
$\displaystyle\sigma_{1}(\bm{Z}_{\star{i}}^{k})\leq\Big{[}(1+6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{k}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\Big{]}\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},$
(21)

where $i=1,2$. Besides, the optimal solution $\bm{P}^{k}$ to
${\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})$ exists and
$\bm{Z}^{k}_{\star{1}}=\bm{Z}_{\star}\bm{P}^{k},\bm{Z}^{k}_{\star{2}}=\bm{Z}_{\star}(\bm{P}^{k})^{-T}$.

Next, we point out the following distance contraction result for the gradient
update that $\tilde{\bm{Z}}^{k+1}=\bm{Z}^{k}-\eta\nabla f^{(k)}(\bm{Z}^{k})$
under conditions (19), (20), and (21).

Figure 2: The comparisons of different algorithms in terms of phase
transitions under frequencies without separation. The red curve is the 90%
success rate curve.

Figure 3: The comparisons of different algorithms in terms of phase
transitions under frequencies with separation. The red curve is the 90%
success rate curve.

######  Lemma 3 (Distance contraction).

Suppose $\bm{M}_{\star}=\mathcal{G}\bm{y}$ is $\mu_{0}$-incoherent, and the
conditions of (19), (20), and (21) hold. Let
$\eta=\frac{\eta^{\prime}}{\sigma_{1}(\bm{M}_{\star})}$ where
$0<\eta^{\prime}\leq\frac{1}{54}$. When $\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu^{2}\kappa^{2}r^{2}\log(n)\right)$, the iterate
$\tilde{\bm{Z}}^{k+1}$ of Algorithm 1 satisfies

$\displaystyle{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})\leq\left(1-\frac{11\eta^{\prime}}{100\kappa}\right){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}),$
(22)

with probability at least $1-c_{2}n^{-2}$, and $\varepsilon_{0}$ is a small
enough constant.

###### Proof.

See Appendix C. ∎

Lemma 3 reveals that when the iterates satisfy conditions about the local
basin of attraction, they enjoy a linear convergence rate to the truth after
the gradient update step. We verify the iterates returned by Algorithm 1
satisfy such conditions via an inductive manner, the details of which are
deferred to subsection IV-D.

###  IV-D Proof of Theorem 1

###### Proof.

We establish Theorem 1 via an inductive approach with the conditions (19),
(20), and (21) as the inductive hypotheses. If we can inductively establish
(19), we finish the proof of Theorem 1. The hypotheses (20) and (21) are
necessary to establish (19) for the next iteration (eg: from the $k$-th to the
$(k+1)$-th step). Therefore we need to establish (19), (20), and (21) together
in an inductive manner.

First, it is evident that the hypotheses are true for $k=0$ from Lemma 2, with
probability at least $1-n^{-2}$ when $\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu\kappa^{4}r^{2}\log(n)\right)$,
$\sigma=\sigma_{1}(\bm{M}_{0})/(1-\varepsilon_{0})$, $\mu\geq\mu_{0}$ and
$\varepsilon_{0}$ is a small enough constant.

Next, we provide the proof sketch from the $k$-step to the $(k+1)$-th step;
see the complete proof in Appendix D-A. The key idea is to invoke Lemma 3 to
establish a linear convergence mechanism and that the projection step makes
the iterates contract in terms of our distance metric, as well as keeping
well-conditioness of $\bm{Z}_{\star{{1}}}$ and $\bm{Z}_{\star{2}}$. Each
induction step holds with probability at least $1-c_{2}n^{-2}$ when
$\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu^{2}\kappa^{2}r^{2}\log(n)\right)$ and
$\eta=\frac{\eta^{\prime}}{\sigma_{1}(\bm{M}_{\star})}$, where
$0<\eta^{\prime}\leq\frac{1}{54}$.

combining the above results, we can inductively establish the following linear
convergence result of Theorem 1

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})\leq\left(1-\frac{11\eta^{\prime}}{100\kappa}\right)^{k}c_{1}\hskip
1.42271pt\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}},$

with probability at least $1-O(n^{-2})$ and $m=(k+1)\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu^{2}\kappa^{4}r^{2}\log(n)\right)$. ∎

##  V Numerical Simulations

In this section, we demonstrate the performance of SHGD via extensive
numerical simulations. Specifically, the signal can not be lifted as a
symmetric Hankel matrix at first glance as the length of the signal is set
even on purpose. The simulations are run in MATLAB R2019b on a 64-bit Windows
machine with multi-core Intel CPU i9-10850K at 3.60 GHz and 16GB RAM. For each
gradient updating of SHGD, we use the whole observation set rather than the
disjoint subsets, as done in [35, 7, 32]. We first compare the SHGD’s
performance in phase transition to nonconvex methods PGD and FIHT, and convex
approaches EMaC and ANM. Then the computational efficiency of SHGD is shown in
Section V-B. We choose the square Hankel version of FIHT to compare, which
further reduces the computation and storage costs. Also, the robustness of
SHGD to additive noise is presented in Section V-C. Lastly, we present the
application of SHGD in delay-Doppler estimation in Section V-D.

###  V-A Phase transition

We compare the performance of SHGD with PGD [6], FIHT [7], EMaC [5], and ANM
[9] in terms of phase transition. The frequency $f_{k}$ is random generated
from $[0,1)$, and the amplitudes $d_{k}$ are selected as
$d_{k}=(1+10^{0.5c_{k}})e^{-i\phi_{k}}$, where $c_{k}$ is uniformly
distributed on $[0,1)$ and $\phi_{k}$ is uniformly sampled from $[0,2\pi)$.

We consider 1-D signals, of which the length is $n=126$. The lifted Hankel
matrix is a $63\times 64$ rectangular matrix for PGD and EMaC. The observed
signal is zero-padded to be a length of the $n=127$ signal and then lifted to
a $64\times 64$ square Hankel matrix for SHGD and FIHT. ANM and EMaC are
implemented using CVX. We use the backtracking line search to choose the
stepsize of SHGD and PGD. SHGD, PGD and FIHT terminate when
$\left\lVert\bm{x}^{k+1}-\bm{x}^{k}\right\rVert_{2}/\left\lVert\bm{x}^{k}\right\rVert_{2}\leq
10^{-7}$ or the maximum number of iterations is reached. We run 50 random
simulations for a prescribed $(r,m)$, where $m=\lfloor pn\rfloor$ is the
number of observations, sample ratio $p$ takes $19$ values from $0.05$ to
$0.95$, and $r$ starts from $r=1$ and increases by $1$ until $r=35$. We run
simulations for two different settings of frequencies, one is that there is no
separation between $\\{f_{k}\\}_{k=1}^{r}$ and the other is that the
separation condition $\Delta=\min_{j\neq k}\|f_{j}-f_{k}\|\geq 1.5/n$ is
satisfied. A test is supposed to be successful when
$\left\lVert\bm{x}^{k}-\bm{x}\right\rVert_{2}/\left\lVert\bm{x}\right\rVert_{2}\leq
10^{-3}$.

From Fig. 2 and Fig. 3, SHGD shows almost the same performance as PGD and
better performance than EMaC in both frequencies without separation and with
separation settings. Also, SHGD outperforms ANM when frequencies are generated
without separation. Compared to FIHT, SHGD shows a better phase transition
performance when the sampling ratio $p<0.5$, which is more typical in reality.

Figure 4: The comparisons of computation time versus relative error. Left:
Frequencies without separation. Right: Frequencies with separation.

Figure 5: The comparisons of computation time versus different problem scales.
Left: The average computation time of three algorithms SHGD, PGD, and FIHT.
Right: The average computation time ratio between SHGD and PGD (SHGD/PGD).

###  V-B Computation time

This subsection conducts simulations to evaluate the computational efficiency
of SHGD. In the first experiment, we examine the computation time versus the
different prescribed recovery error
$\left\lVert\bm{x}^{k}-\bm{x}\right\rVert_{2}/\left\lVert\bm{x}\right\rVert_{2}$
and run 30 random simulations for different nonconvex fast algorithms, namely
SHGD, PGD, and FIHT. We choose the length of the signal as $n=2046$, and the
number of sinusoids as $r=150$, which is a high-dimensional situation. The
number of measurements is set as $m=876$. After zero-padding preprocessing,
the signal is lifted to a $1024\times 1024$ square Hankel matrix for SHGD and
FIHT. We apply a fixed stepsize strategy as
$\eta=\frac{0.75}{\sigma_{1}(\bm{M}_{0})}$ in SHGD and PGD for faster
computing. The simulations are conducted both for frequencies with separation
and without separation. In Fig. 4, the average computation time of SHGD is
approximately half that of PGD in both settings. Additionally, it can be
observed that FIHT exhibits better performance in the initial stage, but
becomes slower than SHGD when higher accuracy is required in both frequency
settings.

In the second experiment, we examine the computation time for SHGD, PGD, and
FIHT under different problem scales. Each problem scale runs 20 random
simulations. The target recovery accuracy is set as
$\left\lVert\bm{x}^{k}-\bm{x}\right\rVert_{2}/\left\lVert\bm{x}\right\rVert_{2}=10^{-7}$.
We set the number of sinusoids as $r=30$, and the frequencies
$\\{f_{k}\\}_{k=1}^{r}$ are generated without separations. The number of
measurements is set as $m=512$. Firstly, we set the problem dimension as
$n=2^{j}-2$ with $j\in\\{11,...,16\\}$. The lifted matrix is a square
$2^{j-1}\times 2^{j-1}$ Hankel matrix for SHGD and FIHT after zero-padding
preprocessing and a $(2^{j-1}-1)\times 2^{j-1}$ rectangular Hankel matrix for
PGD. We examine the average computation time for SHGD, PGD, and FIHT to attain
the fixed accuracy $10^{-7}$. Secondly, we evaluate the time ratio between
SHGD and PGD. And we set the problem dimension as $n=2\lfloor
2^{j-1}\rfloor-2$, where $j$ takes 30 equally spaced values from 11 to 16.
From Fig. 5, we observe that SHGD behaves faster than FIHT for most problem
scenarios. Besides, the computation time of SHGD is about 0.55-0.65 times that
of PGD across different problem dimensions.

###  V-C Robust recovery from noisy observations

In this subsection, we demonstrate the robustness of SHGD to additive noise.
The experiments are conducted both in frequencies with separation and
frequencies without separation settings. We set the length of signals as
$n=127$ and the number of sinusoids as $r=12$. The observations are perturbed
by the noise vector
$\bm{e}=\sigma_{e}\cdot\|\mathcal{P}_{\Omega}(\bm{x})\|_{2}\cdot\bm{w}/\|\bm{w}\|_{2}$
where $\bm{x}$ is referred to the desired signal, $\sigma_{e}$ be the noise
level, and the elements of $\bm{w}$ are i.i.d. standard complex Gaussian
random variables. SHGD is terminated when
$\left\lVert\bm{x}^{k+1}-\bm{x}^{k}\right\rVert_{2}/\left\lVert\bm{x}^{k}\right\rVert_{2}\leq
10^{-7}$. The number of measurements is set as $m=60$ and $m=120$
respectively. The noise level $\sigma_{e}$ varies from $10^{-3}$ to 1,
corresponding to signal-to-noise ratios (SNR) from 60 to 0 dB. For each noise
level, we conduct 20 random tests and record the average RMSE of the
reconstructed signals versus SNR. Fig. 6 shows that the relative
reconstruction error decreases linearly as SNR increases (the noise level
decreases). Besides, the relative reconstruction error reduces when the number
of observations increases.

Figure 6: The stable recovery of SHGD under different noise levels. Left:
Frequencies without separation. Right: Frequencies with separation.

###  V-D SHGD for delay-Doppler estimation based on OFDM signals

We consider the problem of delay-Doppler estimation from OFDM signals [41].
The channel matrix $\bm{D}_{R}$ can be obtained when the reference symbols are
known, which is defined as

$\displaystyle\bm{D}_{R}(l_{1},l_{2})=\sum_{k=1}^{r}d_{k}e^{i2\pi(l_{2}Tf_{k}-l_{1}\Delta
f\tau_{k})},(l_{1},l_{2})\in[N_{1}]\times[N_{2}],$

where $N_{1}$ denotes the number of orthogonal subcarriers, $N_{2}$ denotes
the number of OFDM symbols, $\\{\tau_{k},f_{k}\\}_{k=1}^{r}$ are the delays
and Doppler frequencies, $\\{d_{k}\\}_{k=1}^{r}$ are channel coefficients, and
$T$ is the symbol duration. For simplicity of notation, we define
$\psi_{k}=Tf_{k}\in[0,1)$, $\phi_{k}=\Delta f\tau_{k}\in[0,1)$.

Figure 7:  SHGD for delay-Doppler estimation from OFDM signals. Left: SHGD as
a denoiser when subcarriers of OFDM are equally spaced. Right: SHGD as a
recovery algorithm when the subcarriers of OFDM are unequally spaced.

SHGD can be viewed as a denoiser when subcarriers of OFDM are equally spaced
as this corresponds to the full sampling of the channel matrix $\bm{D}_{R}$.
When subcarriers of OFDM are unequally spaced as [42], it can be interpreted
as the partially non-uniform row sampling of channel matrix $\bm{D}_{R}$, and
we use SHGD to recover the channel matrix. The strengths of non-uniform
subcarriers are lower complexity of the system design and less power of
intercarrier interference (ICI) [42, 43]. We conduct experiments on delay-
Doppler estimations from standard OFDM signals and OFDM signals of non-uniform
subcarriers. After executing SHGD, we use 2-D MUSIC [44] to estimate delay-
Doppler parameters. In simulations, we set $N_{1}=30$, $N_{2}=14$, and the
number of targets $r=4$. We suppose the channel matrix is corrupted by the
noise $\bm{E}=\sigma_{e}\cdot\|\bm{D}_{R}\|_{F}\cdot\bm{W}/\|\bm{W}\|_{F}$
where the entries of $\bm{W}$ are i.i.d standard complex Gaussian random
variables and $\sigma_{e}=0.15$ is the noise level. The number of non-uniform
subcarriers is set as $m=8$ and thus the number of observations of
$\bm{D}_{R}$ is $M=mN_{2}=112$. The 2-D channel matrix $\bm{D}_{R}$ is lifted
to a $128\times 128$ symmetric block-Hankel matrix after zero-padding
preprocessing. Fig. 7 shows that our algorithm SHGD estimates delays, Doppler
frequencies almost exactly, and channel coefficients with minor deviations.

##  VI Conclusions

In this paper, a new non-convex gradient method called SHGD based on symmetric
factorization is proposed for the spectrally sparse signal recovery problem.
Theoretical analysis guarantees that if the number of random samples is larger
than $O(r^{2}\log(n))$, a linear convergence guarantee towards the true
spectrally sparse signal can be established with high probability. The
proposed SHGD algorithm reduces almost half of the cost of computation and
storage compared to PGD and shows comparable computation efficiency as FIHT in
numerical simulations. It is interesting to extend our symmetric factorization
approach to robust spectrally sparse signal recovery and explore vanilla
gradient descent via symmetric factorization when the number of sinusoids is
unknown.

##  Appendix A Supporting lemmas

######  Lemma 4 (Existence of minimizer of the distance metric).

Suppose there exists an invertible matrix $\hat{\bm{P}}$, and set
$\hat{\bm{Z}}_{\star{1}}=\bm{Z}_{\star}\hat{\bm{P}}$,
$\hat{\bm{Z}}_{\star{2}}=\bm{Z}_{\star}\hat{\bm{P}}^{-T}$. When

$\displaystyle\sqrt{\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\|_{F}^{2}+\|\bm{Z}-\hat{\bm{Z}}_{\star{2}}\|_{F}^{2}}\leq\varepsilon\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}},$

where $\varepsilon\leq\sqrt{2}/8$ and
$\sigma_{r}(\hat{\bm{Z}}_{\star{1}})~{}\mathrm{or}~{}\sigma_{r}(\hat{\bm{Z}}_{\star{2}})\geq\frac{\sqrt{2}}{2}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}$,
the minimizer of the problem (12) exists.

###### Proof.

See Appendix E-A. ∎

######  Lemma 5.

For symmetric $\bm{M}_{\star},\bm{M}\in\mathbb{C}^{n_{s}\times n_{s}}$, their
Takagi factorizations are given as
$\bm{M}_{\star}=\bm{U}_{\star}\bm{\Sigma}_{\star}\bm{U}_{\star}^{T}$ and
$\bm{M}=\bm{U}\bm{\Sigma}\bm{U}^{T}$. Let
$\bm{Z}=\bm{U}\bm{\Sigma}^{\frac{1}{2}}$,
$\bm{Z}_{\star}=\bm{U}_{\star}\bm{\Sigma}_{\star}^{\frac{1}{2}}$, then

$\displaystyle{\mbox{dist}^{2}_{P}}(\bm{Z},\bm{Z}_{\star})\leq
2\underset{\bm{Q}\in\mathcal{O}}{\operatorname{min}}\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\|_{F}^{2}\leq\frac{\sqrt{2}+1}{\sigma_{r}(\bm{M}_{\star})}\|\bm{M}-\bm{M}_{\star}\|_{F}^{2},$
(23)

where $\mathcal{O}\triangleq\\{\bm{S}\in\mathbb{R}^{r\times
r}|\bm{S}\bm{S}^{T}=\bm{S}^{T}\bm{S}=\bm{I}\\}$.

###### Proof.

See Appendix E-B. ∎

We remove the dependence on $c_{s}$ in Lemma 6 and Lemma 7 as $c_{s}\leq 2$ in
the square case.

######  Lemma 6 ([7, Lemma 2]).

Assume $\mathcal{G}{\bm{y}}$ is $\mu_{0}$-incoherent. Set
$\bm{L}=\mathcal{T}_{r}\mathcal{G}(p^{-1}\mathcal{P}_{{\hat{\Omega}}}{(\bm{y})})$,
and there exists a universal constant $c_{3}>0$ such that

$\|\bm{L}-\mathcal{G}\bm{y}\|_{2}\leq
c_{3}\sqrt{\frac{\mu_{0}r\log(n)}{m}}\|\mathcal{G}{\bm{y}}\|_{2}$

holds with probability at least $1-n^{-2}$.

######  Lemma 7 ([5, Lemma 3]).

Assume $\mathcal{G}{\bm{y}}$ is $\mu_{0}$-incoherent, and let $T$ be the
tangent space of the rank $r$ matrix manifold at $\mathcal{G}\bm{y}$. Then

$\|\mathcal{P}_{T}\mathcal{G}(\mathcal{I}-p^{-1}\mathcal{P}_{{\hat{\Omega}}})\mathcal{G}^{*}\mathcal{P}_{T}\|_{2}\leq\sqrt{\frac{64\mu_{0}r\log(n)}{m}}$

holds with probability at least $1-n^{-2}$.

######  Remark 6.

Let $\mathcal{G}\bm{y}=\bm{U}\bm{\Sigma}\bm{V}^{H}$ being the SVD of an
$n_{s}\times n_{s}$ matrix $\mathcal{G}\bm{y}$. The subspace $T$ is denoted as

$\displaystyle
T=\\{\bm{U}\bm{C}^{H}+\bm{D}\bm{V}^{H}~{}|~{}\bm{C}\in\mathbb{C}^{n_{s}\times
r},~{}\bm{D}\in\mathbb{C}^{n_{s}\times r}\\}.$

The Takagi factorization of $\mathcal{G}\bm{y}$ is
$\mathcal{G}\bm{y}=\bm{U}_{\star}\bm{\Sigma}_{\star}\bm{U}_{\star}^{T}$ where
$\bm{\Sigma}=\bm{\Sigma}_{\star}$, which is a special form of SVD. Therefore
the subspace $T$ can be reformulated as
$T=\\{\bm{U}_{\star}\bm{C}^{T}+\bm{D}\bm{U}_{\star}^{T}~{}|~{}\bm{C}\in\mathbb{C}^{n_{s}\times
r},~{}\bm{D}\in\mathbb{C}^{n_{s}\times r}\\}$.

######  Lemma 8.

For any fixed matrices $\bm{A}\in\mathbb{C}^{n_{s}\times r}$,
$\bm{B}\in\mathbb{C}^{n_{s}\times r}$, $\bm{C}\in\mathbb{C}^{n_{s}\times r}$,
$\bm{D}\in\mathbb{C}^{n_{s}\times r}$

$\displaystyle\left|\left\langle\mathcal{G}(p^{-1}\mathcal{P}_{\hat{\Omega}}-\mathcal{I})\mathcal{G}^{*}(\bm{A}\bm{B}^{T}),\bm{C}\bm{D}^{T}\right\rangle\right|$
$\displaystyle\leq
c_{4}\sqrt{\frac{n^{2}\log(n)}{m}}\min\\{\left\|\bm{A}\right\|_{F}\left\|\bm{B}\right\|_{2,\infty},\left\|\bm{A}\right\|_{2,\infty}\left\|\bm{B}\right\|_{F}\\}$
$\displaystyle\quad\cdot\min\\{\left\|\bm{C}\right\|_{F}\left\|\bm{D}\right\|_{2,\infty},\left\|\bm{C}\right\|_{2,\infty}\left\|\bm{D}\right\|_{F}\\},$

holds with probability at least $1-n^{-2}$ when $m\gtrsim O(\log(n))$.

###### Proof.

See Appendix E-C. ∎

##  Appendix B Proof of Lemma 2

The proof is partly inspired by [45], and we make adaptations to our (complex)
symmetric Hankel matrix completion problem.

By setting $\bm{L}=\bm{M}^{0}$ in Lemma 6, one has

$\displaystyle\left\|\bm{M}^{0}-\mathcal{G}\bm{y}\right\|_{F}\leq\sqrt{2r}\left\|\bm{M}^{0}-\mathcal{G}\bm{y}\right\|_{2}$
$\displaystyle\leq
c_{3}\sqrt{\frac{2\mu_{0}r^{2}\log(n)}{\hat{m}}}\sigma_{1}(\bm{M}_{\star})$
$\displaystyle\leq\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star}),$ (24)

with probability at least $1-n^{-2}$. The last inequality follows from the
assumption on $\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu\kappa^{4}r^{2}\log(n)\right)$.

As $\bm{M}_{\star}=\mathcal{G}\bm{y}$, a simple corollary of (24) by Weyl’s
Theorem is that

$\displaystyle\sigma_{1}(\bm{M}^{0})\geq(1-\kappa^{-1}\varepsilon_{0})\sigma_{1}(\bm{M}_{\star})\geq(1-\varepsilon_{0})\sigma_{1}(\bm{M}_{\star}).$

Consequently
$\sigma=\frac{\sigma_{1}(\bm{M}^{0})}{1-\varepsilon_{0}}\geq\sigma_{1}(\bm{M}_{\star})$
with high probability.

In Lemma 5, we set
$\bm{M}=\bm{M}^{0}=\bm{U}^{0}\bm{\Sigma}^{0}(\bm{U}^{0})^{T}$, then
$\bm{Z}=\tilde{\bm{Z}}^{0}=\bm{U}^{0}(\bm{\Sigma}^{0})^{\frac{1}{2}}$.
combining Lemma 5 and (24) to obtain

$\displaystyle{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{0},\bm{Z}_{\star})$
$\displaystyle\leq\sqrt{2}\min_{\begin{subarray}{c}\bm{Q}\in\mathcal{O}\end{subarray}}\|\tilde{\bm{Z}}^{0}-\bm{Z}_{\star}\bm{Q}\|_{F}$
$\displaystyle\leq\frac{(\sqrt{2}+1)^{\frac{1}{2}}}{\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}}\|\bm{M}^{0}-\mathcal{G}\bm{y}\|_{F}\leq\frac{8\varepsilon_{0}}{5\kappa}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}.$
(25)

Denote the optimal solution for problem
$\min_{\bm{Q}\in\mathcal{O}}\|\tilde{\bm{Z}}^{0}-\bm{Z}_{\star}\bm{Q}\|_{F}$
as $\bm{Q}_{z}$ and we have the following result from (25)

$\displaystyle\|\tilde{\bm{Z}}^{0}-\bm{Z}_{\star}\bm{Q}_{z}\|\leq
1.6\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}.$

Besides, denote
$\tilde{\bm{Z}}^{0}_{\star{1}}=\bm{Z}_{\star}\tilde{\bm{P}}^{0},\tilde{\bm{Z}}^{0}_{\star{2}}=\bm{Z}_{\star}(\tilde{\bm{P}}^{0})^{-T}$
where $\tilde{\bm{P}}^{0}$ is the optimal solution to
${\mbox{dist}_{P}}(\tilde{\bm{Z}}^{0},\bm{Z}_{\star})$, as defined in (14).
The existence of $\tilde{\bm{P}}^{0}$ can be verified from Lemma 4 by setting
$\hat{\bm{P}}=\bm{Q}_{z}$. Thus

$\displaystyle{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{0},\bm{Z}_{\star})=\sqrt{\|\tilde{\bm{Z}}^{0}-\tilde{\bm{Z}}^{0}_{\star{1}}\|_{F}^{2}+\|\tilde{\bm{Z}}^{0}-\tilde{\bm{Z}}^{0}_{\star{2}}\|_{F}^{2}},$

and

$\displaystyle\|\tilde{\bm{Z}}^{0}-\tilde{\bm{Z}}^{0}_{\star{1}}\|\leq{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{0},\bm{Z}_{\star})\leq
1.6\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}.$

Invoking triangle inequality, we have

$\displaystyle\sigma_{1}\left(\tilde{\bm{Z}}^{0}_{\star{1}}\right)$
$\displaystyle\leq\|\bm{Z}_{\star}\bm{Q}_{z}\|+\|\tilde{\bm{Z}}^{0}-\bm{Z}_{\star}\bm{Q}_{z}\|+\|\tilde{\bm{Z}}^{0}_{\star{1}}-\tilde{\bm{Z}}^{0}\|$
$\displaystyle\leq(1+3.2\kappa^{-1}\varepsilon_{0})\sigma_{1}(\bm{M}_{\star})^{\frac{1}{2}}.$
(26)

Additionally, from Weyl’s theorem

$\displaystyle\sigma_{r}\left(\tilde{\bm{Z}}_{\star{1}}^{0}\right)$
$\displaystyle\geq\sigma_{r}\left(\bm{Z}_{\star}\bm{Q}_{z}\right)-\|\tilde{\bm{Z}}_{\star{1}}^{0}-\bm{Z}_{\star}\bm{Q}_{z}\|$
$\displaystyle\geq\sigma_{r}\left(\bm{Z}_{\star}\right)-\Big{(}\|\tilde{\bm{Z}}^{0}_{\star{1}}-\tilde{\bm{Z}}^{0}\|+\|\tilde{\bm{Z}}^{0}-\bm{Z}_{\star}\bm{Q}_{z}\|\Big{)}$
$\displaystyle\geq\left(1-3.2\kappa^{-1}\varepsilon_{0}\right)\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}.$
(27)

Following a similar route, we can establish the same bounds for
$\sigma_{1}\left(\tilde{\bm{Z}}^{0}_{\star{2}}\right)$ and
$\sigma_{r}\left(\tilde{\bm{Z}}_{\star{1}}^{0}\right)$. Then we have

$\displaystyle\|\tilde{\bm{Z}}^{0}_{\star{1}}\|_{2,\infty}$
$\displaystyle\leq\left\|\bm{U}_{\star}\right\|_{2,\infty}\sigma_{1}\left(\tilde{\bm{Z}}^{0}_{\star{1}}\right)\leq
2\sqrt{\frac{\mu r\sigma}{n}},$

where we invoke (26) and the fact that $\varepsilon_{0}$ is a small constant
such that $(1+3.2\kappa^{-1}\varepsilon_{0})\leq\sqrt{2}$. The same holds for
$\|\tilde{\bm{Z}}^{0}_{\star{2}}\|_{2,\infty}$.

Thus it is evidently that
$\tilde{\bm{Z}}^{0}_{\star{1}},\tilde{\bm{Z}}^{0}_{\star{2}}\in\mathcal{C}$
where $\mathcal{C}$ is a convex set defined in (8). Considering
$\bm{Z}^{0}=\mathcal{P}_{\mathcal{C}}(\tilde{\bm{Z}}^{0})$ and the definition
of ${\mbox{dist}_{P}}(\bm{Z}^{0},\bm{Z}_{\star})$, we immediately have

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{0},\bm{Z}_{\star})$
$\displaystyle\leq\sqrt{\|\bm{Z}^{0}-\tilde{\bm{Z}}^{0}_{\star{1}}\|_{F}^{2}+\|\bm{Z}^{0}-\tilde{\bm{Z}}^{0}_{\star{2}}\|_{F}^{2}}$
$\displaystyle\leq\sqrt{\|\tilde{\bm{Z}}^{0}-\tilde{\bm{Z}}^{0}_{\star{1}}\|_{F}^{2}+\|\tilde{\bm{Z}}^{0}-\tilde{\bm{Z}}^{0}_{\star{2}}\|_{F}^{2}}$
$\displaystyle\leq
1.6\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}.$ (28)

The existence of ${\bm{P}}^{0}$ can be verified from Lemma 4 by setting
$\hat{\bm{P}}=\tilde{\bm{P}}^{0}$. Two immediate results of (28) are

$\displaystyle\|\bm{Z}_{\star{1}}^{0}-\bm{Z}^{0}\|\leq
1.6\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}},$ (29)
$\displaystyle\|\bm{Z}^{0}-\tilde{\bm{Z}}_{\star{1}}^{0}\|\leq
1.6\kappa^{-1}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}.$ (30)

We are still left to prove (16) and (17), which are the bounds for
$\bm{Z}_{\star{1}}^{0}$ and $\bm{Z}_{\star{2}}^{0}$. We only give the proof of
bounding $\bm{Z}_{\star{1}}^{0}$, and the same holds for bounding
$\bm{Z}_{\star{2}}^{0}$.

Invoking triangle inequality and combining (26), (27), (29), (30), we have

$\displaystyle\sigma_{1}\left(\bm{Z}_{\star{1}}^{0}\right)$
$\displaystyle\leq\|\bm{Z}_{\star{1}}^{0}-\bm{Z}^{0}\|+\|\bm{Z}^{0}-\tilde{\bm{Z}}_{\star{1}}^{0}\|+\|\tilde{\bm{Z}}^{0}_{\star{1}}\|$
$\displaystyle\leq\left(1+6.4\kappa^{-1}\varepsilon_{0}\right)\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},$

and additionally invoking Weyl’s theorem

$\displaystyle\sigma_{r}\left(\bm{Z}_{\star{1}}^{0}\right)$
$\displaystyle\geq\sigma_{r}(\tilde{\bm{Z}}_{\star{1}}^{0})-\|\bm{Z}_{\star{1}}^{0}-\tilde{\bm{Z}}_{\star{1}}^{0}\|$
$\displaystyle\geq\sigma_{r}(\tilde{\bm{Z}}_{\star{1}}^{0})-\Big{(}\|\bm{Z}_{\star{1}}^{0}-\bm{Z}^{0}\|+\|\bm{Z}^{0}-\tilde{\bm{Z}}_{\star{1}}^{0}\|\Big{)}$
$\displaystyle\geq\left(1-6.4\kappa^{-1}\varepsilon_{0}\right)\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}.$

##  Appendix C Proof of Lemma 3

As shown in [34, 35, 46], it is sufficient to study the sampling model with
replacement. In the following derivations, ${\hat{\Omega}}_{k}$ is a sampling
set with replacement corresponding to the index set $\Omega_{k}$ in the $k$-th
iteration. The following facts are immediate results of (20) and (21), which
are helpful for establishing (22).

######  Claim 1.

Under the conditions of (20) and (21),

$\displaystyle\frac{1}{2}\sigma_{r}(\bm{M}_{\star})\leq\sigma_{r}^{2}\left(\bm{Z}_{\star{i}}^{k}\right)$
$\displaystyle\leq\sigma_{1}^{2}\left(\bm{Z}_{\star{i}}^{k}\right)\leq
2\sigma_{1}(\bm{M}_{\star}).$ (31)

###### Proof.

There we only give the proof of (31) for $i=1$, as it follows a similar route
when $i=2$. From (21)

$\displaystyle\sigma_{1}(\bm{Z}_{\star{1}}^{k})$
$\displaystyle\leq\big{(}1+6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{\infty}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\big{)}\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}$
$\displaystyle=(1+\frac{640\varepsilon_{0}}{11\eta^{\prime}})\sigma_{1}(\bm{M}_{\star})^{\frac{1}{2}}\leq\sqrt{2}\sigma_{1}(\bm{M}_{\star})^{\frac{1}{2}},$

where the last inequality is from that $\varepsilon_{0}$ is small enough such
that $\frac{640\varepsilon_{0}}{11\eta^{\prime}}\leq\sqrt{2}-1$ and
$\eta^{\prime}$ is some constant. From (20)

$\displaystyle\sigma_{r}(\bm{Z}_{\star{1}}^{k})$
$\displaystyle\geq\big{(}1-6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{\infty}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\big{)}\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}$
$\displaystyle=(1-\frac{640\varepsilon_{0}}{11\eta^{\prime}})\geq\frac{\sqrt{2}}{2}\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},$

as long as $\varepsilon_{0}$ is a small constant such that
$\frac{640\varepsilon_{0}}{11\eta^{\prime}}\leq 1-\frac{\sqrt{2}}{2}$. ∎

Now we begin to establish (22). According to the definition of the distance
metric (12)

$\displaystyle{\mbox{dist}^{2}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})$
$\displaystyle\leq\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|_{F}^{2}+\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{2}}^{k}\|_{F}^{2}$
$\displaystyle={\mbox{dist}^{2}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})-2\eta\operatorname{{Re}}\langle\nabla
f^{(k)}(\bm{Z}^{k}),\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\rangle$
$\displaystyle\quad+2\eta^{2}\|\nabla f^{(k)}(\bm{Z}^{k})\|_{F}^{2},$ (32)

where the equality results from $\tilde{\bm{Z}}^{k+1}=\bm{Z}^{k}-\eta\nabla
f^{(k)}(\bm{Z}^{k})$ and we set
$\bm{\Delta_{1}}^{k}=\bm{Z}^{k}-\bm{Z}_{\star{1}}^{k}$,
$\bm{\Delta_{2}}^{k}=\bm{Z}^{k}-\bm{Z}_{\star{2}}^{k}$. Next we build a lower
bound for $\operatorname{{Re}}\langle\nabla
f^{(k)}(\bm{Z}^{k}),\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\rangle$ and an
upper bound for $\|\nabla f^{(k)}(\bm{Z}^{k})\|_{F}^{2}$ in step 1 and step 2.
Finally, we give the upper bound of (32) in step 3.

Step 1: bounding $\operatorname{{Re}}\left\langle\nabla
f^{(k)}(\bm{Z}^{k}),\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\right\rangle$. We
reformulate this term into two parts to bound.

$\displaystyle\operatorname{{Re}}\left\langle\nabla
f^{(k)}(\bm{Z}^{k}),\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\right\rangle$
$\displaystyle=\operatorname{{Re}}\langle(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})\bar{\bm{Z}}^{k}$
$\displaystyle\quad+\hat{p}^{-1}\mathcal{G}\mathcal{P}_{{\hat{\Omega}}_{k}}\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})\bar{\bm{Z}}^{k},\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\rangle$
$\displaystyle=\operatorname{{Re}}\langle\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star},\bm{\Delta_{1}}^{k}(\bm{Z}^{k})^{T}+\bm{Z}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle+\operatorname{{Re}}$
$\displaystyle\quad\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}),\bm{\Delta_{1}}^{k}(\bm{Z}^{k})^{T}+\bm{Z}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle=I_{1}+I_{2},$

where the first equality results from
$(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})=(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}^{k}(\bm{Z}^{k})^{T})$,
and the second equality results from the symmetry of
$\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}$. In the following derivations, we
invoke the decomposition as

$\displaystyle\bm{\Delta_{1}}^{k}(\bm{Z}^{k})^{T}+\bm{Z}^{k}(\bm{\Delta_{2}}^{k})^{T}$
$\displaystyle=\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+2\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T},$
$\displaystyle\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}$
$\displaystyle=\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}.$

Thus for $I_{1}$, we have

$\displaystyle I_{1}$
$\displaystyle=\operatorname{{Re}}\langle\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T},$
$\displaystyle\quad\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+2\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle=\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}+2\|\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\quad+3\operatorname{{Re}}\langle\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T},\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle\geq\frac{3}{4}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}{+}\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}-7\|\bm{\Delta_{1}}^{k}\|_{F}^{2}\|\bm{\Delta_{2}}^{k}\|_{F}^{2}$
$\displaystyle\geq\frac{3}{4}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}{+}\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\quad-14\varepsilon_{0}^{2}\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2}),$

where the third inequality results from
$3\operatorname{{Re}}\langle\bm{A},\bm{B}\rangle=\|\frac{1}{2}\bm{A}+3\bm{B}\|_{F}^{2}-\frac{1}{4}\|\bm{A}\|_{F}^{2}-9\|\bm{B}\|_{F}^{2}\geq-\frac{1}{4}\|\bm{A}\|_{F}^{2}-9\|\bm{B}\|_{F}^{2}$.
The last inequality follows from

$\displaystyle\mathrm{max}\\{\|\bm{\Delta_{1}}^{k}\|_{F},\|\bm{\Delta_{2}}^{k}\|_{F}\\}{\leq}{\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}){\leq}2\varepsilon_{0}\sqrt{\sigma_{r}(\bm{M}_{\star})}.$
(33)

For $I_{2}$, we have

$\displaystyle I_{2}$
$\displaystyle=\operatorname{{Re}}\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),$
$\displaystyle\quad\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+2\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle=\operatorname{{Re}}\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),$
$\displaystyle\quad\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle+3\operatorname{{Re}}\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle+2\operatorname{{Re}}\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle=I_{2,1}+3I_{2,2}+2I_{2,3}.$

It is easily to see
$\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\in
T$. Invoking Lemma 7 when $\hat{m}\gtrsim O\left(\varepsilon_{0}^{-2}\mu
r\log(n)\right)$, we have

$\displaystyle I_{2,1}$
$\displaystyle\geq-\|\mathcal{P}_{T}\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}\mathcal{P}_{T}\|\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\geq-\varepsilon_{0}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}.$

We have the following results from Claim 1

$\displaystyle\|\bm{Z}_{\star{1}}^{k}\|_{2,\infty}\leq\|\bm{U}_{\star}\|_{2,\infty}\|\bm{Z}_{\star{1}}^{k}\|\leq
2\sqrt{\frac{\mu r\sigma}{n}},$
$\displaystyle\|\bm{\Delta_{1}}^{k}\|_{2,\infty}\leq\|\bm{Z}^{k}\|_{2,\infty}+\|\bm{Z}_{\star{1}}^{k}\|_{2,\infty}\leq
4\sqrt{\frac{\mu r\sigma}{n}},$

and similarly, the bounds hold for $\|\bm{Z}_{\star{2}}^{k}\|_{2,\infty}$ and
$\|\bm{\Delta_{2}}^{k}\|_{2,\infty}$.

combining the previous results and invoking Lemma 8

$\displaystyle I_{2,2}$
$\displaystyle=\operatorname{{Re}}\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}\rangle$
$\displaystyle\quad+\operatorname{{Re}}\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle\geq-
c_{3}\sqrt{\frac{n^{2}\log(n)}{\hat{m}}}(\|\bm{\Delta_{1}}^{k}\|_{2,\infty}\|\bm{\Delta_{2}}^{k}\|_{F}\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{Z}_{\star{2}}^{k}\|_{2,\infty}$
$\displaystyle\quad+\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{2,\infty}\|\bm{Z}_{\star{1}}^{k}\|_{2,\infty}\|\bm{\Delta_{2}}^{k}\|_{F})$
$\displaystyle\overset{(a)}{\geq}-16c_{3}\sqrt{\frac{\mu^{2}r^{2}(1+\varepsilon_{0})^{2}\log(n)\sigma_{1}^{2}(\bm{M}_{\star})}{\hat{m}(1-\varepsilon_{0})^{2}}}\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{F}$
$\displaystyle\overset{(b)}{\geq}-\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{F},$

$(a)$ results from
$\sigma=\frac{\sigma_{1}(\bm{M}^{0})}{1-\varepsilon_{0}}\leq\frac{1+\varepsilon_{0}}{1-\varepsilon_{0}}\sigma_{1}(\bm{M}_{\star})$,
where $\varepsilon_{0}$ is a small constant. $(b)$ follows from
$\hat{m}\gtrsim O(\varepsilon_{0}^{-2}\mu^{2}\kappa^{2}r^{2}\log(n))$.

Similarly invoking Lemma 8, one can obtain

$\displaystyle
I_{2,3}\geq-\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{F}.$

Therefore

$\displaystyle I_{2}$
$\displaystyle{\geq}{-}\varepsilon_{0}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}{+}\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}{-}5\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{F}$
$\displaystyle{\geq}{-}\varepsilon_{0}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}{+}\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}{-}\frac{5}{2}\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2}),$

and combine $I_{1}$, $I_{2}$ to obtain

$\displaystyle\operatorname{{Re}}\langle\nabla
f^{(k)}(\bm{Z}),\bm{\Delta_{1}}^{k}+$
$\displaystyle\bm{\Delta_{2}}^{k}\rangle\geq(\frac{3}{4}-\varepsilon_{0})\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle-(14\varepsilon_{0}^{2}+\frac{5}{2}\varepsilon_{0})\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2}).$

We point out a result towards the first term in the RHS of the previous
inequality.

$\displaystyle\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle=\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}\|_{F}^{2}+\|\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}+2\operatorname{{Re}}\langle\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T},\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle$
$\displaystyle\overset{(a)}{\geq}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}\|_{F}^{2}+\|\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\geq\sigma_{r}^{2}(\bm{Z}_{\star{1}}^{k})\|\bm{\Delta_{2}}^{k}\|_{F}^{2}+\sigma_{r}^{2}(\bm{Z}_{\star{2}}^{k})\|\bm{\Delta_{1}}^{k}\|_{F}^{2}$
$\displaystyle\overset{(b)}{\geq}\frac{1}{2}\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2}),$
(34)

where $(a)$ results from a corollary of the first-order optimality condition
(13):
$(\bm{Z}_{\star{1}}^{k})^{H}\bm{\Delta_{1}}^{k}=(\bm{\Delta_{2}}^{k})^{T}{\overline{\bm{Z}_{\star{2}}^{k}}}$,
which is

$\displaystyle\langle\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T},\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\rangle=\langle(\bm{Z}_{\star{1}}^{k})^{H}\bm{\Delta_{1}}^{k},(\bm{\Delta_{2}}^{k})^{T}{\overline{\bm{Z}_{\star{2}}^{k}}}\rangle\geq
0,$

and $(b)$ comes from Claim 1.

Step 2: bounding $\|\nabla f^{(k)}(\bm{Z}^{k})\|_{F}^{2}$. We bound this term
by controlling two parts $I_{3}$ and $I_{4}$.

$\displaystyle\|\nabla f^{(k)}(\bm{Z}^{k})\|_{F}^{2}$
$\displaystyle{=}\|\big{(}(\bm{Z}^{k}(\bm{Z}^{k})^{T}{-}\bm{M}_{\star}){+}\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}{-}\bm{M}_{\star})\big{)}\bar{\bm{Z}}^{k}\|_{F}^{2}$
$\displaystyle\leq
2\|(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})\bar{\bm{Z}}^{k}\|_{F}^{2}$
$\displaystyle\quad+2\|\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})\bar{\bm{Z}}^{k}\|_{F}^{2}=2(I_{3}^{2}+I_{4}^{2}),$

where the first equality results from the fact
$(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})=(\mathcal{I}-\mathcal{G}\mathcal{G}^{*})(\bm{Z}^{k}(\bm{Z}^{k})^{T})$.

$\displaystyle I_{3}^{2}$
$\displaystyle\leq\|\bm{Z}^{k}\|^{2}\|\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}\|_{F}^{2}$
$\displaystyle\leq
2\|\bm{Z}^{k}\|^{2}\left(\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}+\|\bm{\Delta_{1}}^{k}\|_{F}^{2}\|\bm{\Delta_{2}}^{k}\|_{F}^{2}\right),$
$\displaystyle I_{4}$
$\displaystyle=\|\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star})\bar{\bm{Z}}^{k}\|_{F}$
$\displaystyle=\max_{\|\tilde{\bm{X}}\|_{F}=1}|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}),\tilde{\bm{X}}(\bm{Z}^{k})^{T}\rangle|.$

Note that

$\displaystyle|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\bm{M}_{\star}),\tilde{\bm{X}}(\bm{Z}^{k})^{T}\rangle|$
$\displaystyle=|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}+\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),$
$\displaystyle\quad\tilde{\bm{X}}(\bm{\Delta_{2}}^{k})^{T}+\tilde{\bm{X}}(\bm{Z}_{\star{2}}^{k})^{T}\rangle|$
$\displaystyle\leq|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\tilde{\bm{X}}(\bm{Z}_{\star{2}}^{k})^{T}\rangle|$
$\displaystyle\quad+|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\tilde{\bm{X}}(\bm{\Delta_{2}}^{k})^{T}\rangle|$
$\displaystyle\quad+|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\tilde{\bm{X}}(\bm{\Delta_{2}}^{k})^{T}\rangle|$
$\displaystyle\quad+|\langle\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}(\bm{\Delta_{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}),\tilde{\bm{X}}(\bm{Z}_{\star{2}}^{k})^{T}\rangle|$
$\displaystyle=I_{4,1}+I_{4,2}+I_{4,3}+I_{4,4}.$

As
$\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\in
T$ and $\tilde{\bm{X}}(\bm{Z}_{\star{2}}^{k})^{T}\in T$, we invoke Lemma 7 to
obtain the upper bound of $I_{4,1}$ when $\|\tilde{\bm{X}}\|_{F}=1$, similarly
as bounding $I_{2,1}$

$\displaystyle I_{4,1}$
$\displaystyle\leq\|\mathcal{P}_{T}\mathcal{G}(\hat{p}^{-1}\mathcal{P}_{{\hat{\Omega}}_{k}}-\mathcal{I})\mathcal{G}^{*}\mathcal{P}_{T}\|$
$\displaystyle\quad\cdot\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}\cdot\|\tilde{\bm{X}}(\bm{Z}_{\star{2}}^{k})^{T}\|_{F}$
$\displaystyle\leq\varepsilon_{0}\sqrt{2\sigma_{1}(\bm{M}_{\star})}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}.$

Similarly as bounding $I_{2,2}$, when $\hat{m}\gtrsim
O(\varepsilon_{0}^{-2}\mu^{2}\kappa^{2}r^{2}\log(n))$ and
$\|\tilde{\bm{X}}\|_{F}=1$, invoke Lemma 8 to obtain

$\displaystyle I_{4,2}$ $\displaystyle\leq
c_{3}\sqrt{\frac{n^{2}\log(n)}{\hat{m}}}(\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{Z}_{\star{2}}^{k}\|_{2,\infty}\|\tilde{\bm{X}}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{2,\infty}$
$\displaystyle\quad+\|\bm{Z}_{\star{1}}^{k}\|_{2,\infty}\|\bm{\Delta_{2}}^{k}\|_{F}\|\tilde{\bm{X}}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{2,\infty})$
$\displaystyle\leq\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}+\|\bm{\Delta_{2}}^{k}\|_{F}),$
$\displaystyle I_{4,3}$ $\displaystyle\leq
c_{3}\sqrt{\frac{n^{2}\log(n)}{\hat{m}}}\|\bm{\Delta_{1}}^{k}\|_{2,\infty}\|\bm{\Delta_{2}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{2,\infty}\|\tilde{\bm{X}}\|_{F}$
$\displaystyle\leq\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})\|\bm{\Delta_{2}}^{k}\|_{F},$
$\displaystyle I_{4,4}$ $\displaystyle\leq
c_{3}\sqrt{\frac{n^{2}\log(n)}{\hat{m}}}\|\bm{\Delta_{1}}^{k}\|_{F}\|\bm{\Delta_{2}}^{k}\|_{2,\infty}\|\bm{Z}_{\star{2}}^{k}\|_{2,\infty}\|\tilde{\bm{X}}\|_{F}$
$\displaystyle\leq\varepsilon_{0}\sigma_{r}(\bm{M}_{\star})\|\bm{\Delta_{1}}^{k}\|_{F}.$

Combine $I_{4,1}$, $I_{4,2}$, $I_{4,3}$ and $I_{4,4}$ to obtain

$\displaystyle I_{4}^{2}$ $\displaystyle\leq
4\varepsilon_{0}^{2}\sigma_{1}(\bm{M}_{\star})\Big{(}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\quad+4\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2})\Big{)}.$

Consequently, combine $I_{3}^{2}$ and $I_{4}^{2}$ to obtain

$\displaystyle\|\nabla f^{(k)}(\bm{Z}^{k})\|_{F}^{2}$
$\displaystyle\overset{(a)}{\leq}\sigma_{1}(\bm{M}_{\star})\Big{(}(10+8\varepsilon_{0}^{2})\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\quad+10\|\bm{\Delta_{1}}^{k}\|_{F}^{2}\|\bm{\Delta_{2}}^{k}\|_{F}^{2}+32\varepsilon_{0}^{2}\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2})\Big{)}$
$\displaystyle\overset{(b)}{\leq}\sigma_{1}(\bm{M}_{\star})\Big{(}(10+8\varepsilon_{0}^{2})\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}$
$\displaystyle\quad+52\varepsilon_{0}^{2}\sigma_{r}(\bm{M}_{\star})(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2})\Big{)},$

where step $(a)$ comes from the fact that

$\displaystyle\|\bm{Z}^{k}\|^{2}$
$\displaystyle\leq\left(\|\bm{\Delta_{1}}^{k}\|+\|\bm{Z}_{\star{1}}^{k}\|\right)^{2}\leq
2.5\sigma_{1}(\bm{M}_{\star}),$

where we invoke the results (33) and Claim 1 for $\varepsilon_{0}$ is a small
constant. Step $(b)$ follows from (33) directly.

Finally, we give the upper bound of (32) by combining step 1 and step 2.

Step 3: bounding (32). Combining the above bounds for
$\operatorname{{Re}}\left\langle\nabla
f^{(k)}(\bm{Z}^{k}),\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\right\rangle$ and
$\|\nabla f^{(k)}(\bm{Z}^{k})\|_{F}^{2}$, we have

$\displaystyle-2\eta\operatorname{{Re}}\left\langle\nabla
f^{(k)}(\bm{Z}^{k}),\bm{\Delta_{1}}^{k}+\bm{\Delta_{2}}^{k}\right\rangle+2\eta^{2}\|\nabla
f^{(k)}(\bm{Z}^{k})\|_{F}^{2}$
$\displaystyle\leq\frac{\alpha}{\sigma_{1}(\bm{M}_{\star})}\|\bm{\Delta_{1}}^{k}(\bm{Z}_{\star{2}}^{k})^{T}+\bm{Z}_{\star{1}}^{k}(\bm{\Delta_{2}}^{k})^{T}\|_{F}^{2}+\frac{\beta}{\kappa}(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2})$
$\displaystyle\overset{(a)}{\leq}\frac{\alpha}{2\kappa}(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2})+\frac{\beta}{\kappa}(\|\bm{\Delta_{1}}^{k}\|_{F}^{2}+\|\bm{\Delta_{2}}^{k}\|_{F}^{2})$
$\displaystyle\overset{(b)}{\leq}-\frac{11\eta^{\prime}}{50\kappa}{\mbox{dist}^{2}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}),$

where $\eta=\frac{\eta^{\prime}}{\sigma_{1}(\bm{M}_{\star})}$,
$\alpha=4{\eta^{\prime}}^{2}(5+4\varepsilon_{0}^{2})-\eta^{\prime}(1.5-2\varepsilon_{0})$
and
$\beta=\eta^{\prime}(28\varepsilon_{0}^{2}+5\varepsilon_{0})+104{\eta^{\prime}}^{2}\varepsilon_{0}^{2}$.
For $\eta^{\prime}\leq\frac{1}{54}$, we list the following numerical results
for $\varepsilon_{0}$ is a small enough constant

$\displaystyle\alpha$
$\displaystyle=4{\eta^{\prime}}^{2}(5+4\varepsilon_{0}^{2})-\eta^{\prime}(1.5-2\varepsilon_{0})\leq
0,$ $\displaystyle\frac{\alpha}{2}+\beta$
$\displaystyle=2{\eta^{\prime}}^{2}(5+56\varepsilon_{0}^{2})+\eta^{\prime}(28\varepsilon_{0}^{2}+6\varepsilon_{0}-0.75)\leq-\frac{11\eta^{\prime}}{50}.$

Therefore step $(a)$ comes from (34) and $\alpha\leq 0$, step $(b)$ follows
from $\frac{\alpha}{2}+\beta\leq-\frac{11\eta^{\prime}}{50}$.

Finally, we obtain an upper bound of (32)

$\displaystyle{\mbox{dist}^{2}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})$
$\displaystyle\leq\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|_{F}^{2}+\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{2}}^{k}\|_{F}^{2}$
$\displaystyle\leq(1-\frac{11\eta^{\prime}}{50\kappa}){\mbox{dist}^{2}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}),$
(35)

and further establish (22)

$\displaystyle{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})\leq(1-\frac{11\eta^{\prime}}{100\kappa}){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}).$

## Appendix D

###  D-A Proof of the induction hypotheses from the $k$-th to the $(k+1)$-th
step

As (19), (20), and (21) hold for the $k$-th step, we have the following linear
convergence result by invoking Lemma 3

$\displaystyle{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})\leq(1-\frac{11\eta^{\prime}}{100\kappa}){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}),$
(36)

with probability at least $1-c_{2}n^{-2}$ when $\hat{m}\gtrsim
O\left(\varepsilon_{0}^{-2}\mu^{2}\kappa^{2}r^{2}\log(n)\right)$,
$\eta=\frac{\eta^{\prime}}{\sigma_{1}(\bm{M}_{\star})}$, $\eta^{\prime}$ and
$\varepsilon_{0}$ are some appropriate constants. We establish (19) for the
$(k+1)$-th step by proving the following non-expansiveness property of the
projection in terms of the distance metric:

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{k+1},\bm{Z}_{\star})\leq{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star}),$

where $\bm{Z}^{k+1}=\mathcal{P}_{\mathcal{C}}(\tilde{\bm{Z}}^{k+1})$.

Set $\tilde{\bm{Z}}_{\star{1}}^{k+1}=\bm{Z}_{\star}\tilde{\bm{P}}^{k+1}$,
$\tilde{\bm{Z}}_{\star{2}}^{k+1}=\bm{Z}_{\star}(\tilde{\bm{P}}^{k+1})^{-T}$
where $\tilde{\bm{P}}^{k+1}$ is the optimal solution as defined in (14).
According to (36) and the definition of
${\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})$

$\displaystyle\|\tilde{\bm{Z}}^{k+1}-\tilde{\bm{Z}}_{\star{1}}^{k+1}\|\leq{\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star})\leq(1-\frac{11\eta^{\prime}}{100\kappa}){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}).$
(37)

One simple corollary of (35) is that

$\displaystyle\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|$
$\displaystyle\leq(1-\frac{11\eta^{\prime}}{100\kappa}){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}).$
(38)

Invoking triangle inequality and combining (37),(38), and (21), we obtain that

$\displaystyle\sigma_{1}(\tilde{\bm{Z}}_{\star{1}}^{k+1})$
$\displaystyle\leq\|\bm{Z}_{\star{1}}^{k}\|+\|\tilde{\bm{Z}}_{\star{1}}^{k+1}-\tilde{\bm{Z}}^{k+1}\|+\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|$
$\displaystyle\leq\Big{[}1+6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{\infty}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\Big{]}\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}$
$\displaystyle=(1+\frac{640\varepsilon_{0}}{11\eta^{\prime}})\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}\leq\sqrt{2}\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}},$

when $\varepsilon_{0}$ is a small constant such that
$\frac{640\varepsilon_{0}}{11\eta^{\prime}}\leq\sqrt{2}-1$.

As $\tilde{\bm{Z}}_{\star{1}}^{k+1}=\bm{Z}_{\star}\tilde{\bm{P}}^{k+1}$, one
can obtain

$\displaystyle\|\tilde{\bm{Z}}_{\star{1}}^{k+1}\|_{2,\infty}\leq\|\bm{U}_{\star}\|_{2,\infty}\cdot\sigma_{1}(\tilde{\bm{Z}}_{\star{1}}^{k+1})\leq
2\sqrt{\frac{\mu r\sigma}{n}}.$

Thus $\tilde{\bm{Z}}_{\star{1}}^{k+1}\in\mathcal{C}$ and one can also obtain
$\tilde{\bm{Z}}_{\star{2}}^{k+1}\in\mathcal{C}$ , following a similar route.

combining the definition of the distance metric and the non-expansiveness of
convex projection $\mathcal{P}_{\mathcal{C}}$, we have

$\displaystyle{\mbox{dist}_{P}}(\bm{Z}^{k+1},\bm{Z}_{\star})$
$\displaystyle{\leq}\sqrt{\|\bm{Z}^{k+1}-\tilde{\bm{Z}}_{\star{1}}^{k+1}\|_{F}^{2}+\|\bm{Z}^{k+1}-\tilde{\bm{Z}}_{\star{2}}^{k+1}\|_{F}^{2}}$
$\displaystyle{\leq}\sqrt{\|\tilde{\bm{Z}}^{k+1}-\tilde{\bm{Z}}_{\star{1}}^{k+1}\|_{F}^{2}+\|\tilde{\bm{Z}}^{k+1}-\tilde{\bm{Z}}_{\star{2}}^{k+1}\|_{F}^{2}}$
$\displaystyle={\mbox{dist}_{P}}(\tilde{\bm{Z}}^{k+1},\bm{Z}_{\star}).$ (39)

The things left are to prove (20) and (21) for the $(k+1)$-th step. We only
give the proof of bounding $\bm{Z}_{\star{1}}^{k+1}$ in details, because
bounding $\bm{Z}_{\star{2}}^{k+1}$ follows a similar route.

From (36) and (39), one can obtain that

$\displaystyle\|\bm{Z}_{\star{1}}^{k+1}-\bm{Z}^{k+1}\|$
$\displaystyle\leq{\mbox{dist}_{P}}(\bm{Z}^{k+1},\bm{Z}_{\star})\leq(1-\frac{11\eta^{\prime}}{100\kappa}){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}),$
(40) $\displaystyle\|\bm{Z}^{k+1}-\tilde{\bm{Z}}_{\star{1}}^{k+1}\|$
$\displaystyle\leq(1-\frac{11\eta^{\prime}}{100\kappa}){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}).$
(41)

Invoking the triangle inequality for several times, and combining the results
(37), (38), (40), and (41), we have

$\displaystyle\sigma_{1}(\bm{Z}_{\star{1}}^{k+1})$
$\displaystyle=\|\bm{Z}_{\star{1}}^{k+1}\|$
$\displaystyle\leq\|\bm{Z}_{\star{1}}^{k}\|+\|\bm{Z}_{\star{1}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|$
$\displaystyle\leq\sigma_{1}(\bm{Z}_{\star{1}}^{k})+\|\bm{Z}_{\star{1}}^{k+1}-\bm{Z}^{k+1}\|+\|\bm{Z}^{k+1}-\tilde{\bm{Z}}_{\star{1}}^{k+1}\|$
$\displaystyle\quad+\|\tilde{\bm{Z}}_{\star{1}}^{k+1}-\tilde{\bm{Z}}^{k+1}\|+\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|$
$\displaystyle\leq\Big{[}1+6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{k+1}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\Big{]}\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}.$

Besides, invoke Weyl’s Theorem to obtain that

$\displaystyle\sigma_{r}\left(\bm{Z}_{\star{1}}^{k+1}\right)$
$\displaystyle\geq\sigma_{r}\left(\bm{Z}_{\star{1}}^{k}\right)-\|\bm{Z}_{\star{1}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|$
$\displaystyle\geq\sigma_{r}\left(\bm{Z}_{\star{1}}^{k}\right)-\Big{(}\|\bm{Z}_{\star{1}}^{k+1}-\bm{Z}^{k+1}\|+\|\bm{Z}^{k+1}-\tilde{\bm{Z}}_{\star{1}}^{k+1}\|$
$\displaystyle\quad+\|\tilde{\bm{Z}}_{\star{1}}^{k+1}-\tilde{\bm{Z}}^{k+1}\|+\|\tilde{\bm{Z}}^{k+1}-\bm{Z}_{\star{1}}^{k}\|\Big{)}$
$\displaystyle\geq\Big{[}1-6.4\kappa^{-1}\varepsilon_{0}\sum_{t=0}^{k+1}(1-\frac{11\eta^{\prime}}{100\kappa})^{t}\Big{]}\sigma_{r}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}.$

Therefore, (20) and (21) for the $(k+1)$-step are established.

###  D-B Proof of Remark 2

We establish a relationship between the recovery error and our distance
metric, which is

$\displaystyle\|\bm{x}^{k}-\bm{x}\|_{2}$
$\displaystyle\leq\|\bm{y}^{k}-\bm{y}\|_{2}=\|\mathcal{G}^{*}(\bm{Z}^{k}(\bm{Z}^{k})^{T}-\mathcal{G}\bm{y})\|_{2}$
$\displaystyle\leq\|\bm{Z}^{k}{(\bm{Z}^{k})}^{T}-\bm{M}_{\star}\|_{F}$
$\displaystyle=\|(\bm{Z}^{k}-\bm{Z}^{k}_{\star{1}})(\bm{Z}^{k}_{\star{2}})^{T}+\bm{Z}^{k}(\bm{Z}^{k}-\bm{Z}^{k}_{\star{2}})^{T})\|_{F}$
$\displaystyle\leq(\|\bm{Z}^{k}_{\star{2}}\|+\|\bm{Z}^{k}\|){\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})$
$\displaystyle\lesssim\sigma_{1}(\bm{M}_{\star})^{\frac{1}{2}}{\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star}).$

The last inequality invokes results from the induction hypotheses (19) and
(21), which are
$\|\bm{Z}^{k}_{\star{2}}\|\lesssim\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}$
and thus
$\|\bm{Z}^{k}\|\leq\|\bm{Z}^{k}-\bm{Z}^{k}_{\star{2}}\|+\|\bm{Z}^{k}_{\star{2}}\|\lesssim\sigma_{1}\left(\bm{M}_{\star}\right)^{\frac{1}{2}}$.
Therefore, $\|\bm{x}^{k}-\bm{x}\|_{2}\leq O(\varepsilon)$ when
${\mbox{dist}_{P}}(\bm{Z}^{k},\bm{Z}_{\star})\leq\varepsilon$.

###  D-C Proof of Lemma 1

From the definitions of $\mbox{dist}^{2}_{Q}(\bm{Z},\bm{Z}_{\star})$ and
$\mbox{dist}^{2}_{P}(\bm{Z},\bm{Z}_{\star})$

$\displaystyle 2\mbox{dist}^{2}_{Q}(\bm{Z},\bm{Z}_{\star})$
$\displaystyle\overset{(a)}{\geq}\inf_{\begin{subarray}{c}\bm{Q}\in\mathbb{C}^{r\times
r},\\\
\text{invertible}\end{subarray}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\right\|_{F}^{2}+\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}^{-T}\right\|_{F}^{2}$
$\displaystyle=\mbox{dist}^{2}_{P}(\bm{Z},\bm{Z}_{\star}),$

where $(a)$ results from that the domain becomes larger and the first part is
proven.

Then we prove the second part. For any invertible $\bm{P}$ such that
$\sqrt{\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}\right\|_{F}^{2}+\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}^{-T}\right\|_{F}^{2}}\leq\varepsilon\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}$,
one has

$\displaystyle\|\bm{Z}_{\star}(\bm{P}-\bm{P}^{-T})\|_{F}$
$\displaystyle\leq\|\bm{Z}-\bm{Z}_{\star}\bm{P}\|_{F}+\|\bm{Z}-\bm{Z}_{\star}\bm{P}^{-T}\|_{F}$
$\displaystyle\leq 2\varepsilon\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}},$

where we invoke the triangle inequality. From the relation that
$\|\bm{A}\bm{B}\|_{F}\geq\sigma_{r}(\bm{B})\|\bm{A}\|_{F}$, we have

$\displaystyle\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}\|\bm{P}-\bm{P}^{-T}\|_{F}\leq\|\bm{Z}_{\star}(\bm{P}-\bm{P}^{-T})\|_{F},$

therefore $\|\bm{P}-\bm{P}^{-T}\|_{F}\leq 2\varepsilon$, and under this
condition, $\mbox{dist}^{2}_{P}(\bm{Z},\bm{Z}_{\star})$ can be reformulated as

$\displaystyle\mbox{dist}^{2}_{P}(\bm{Z},\bm{Z}_{\star}){=}\inf_{\begin{subarray}{c}\|\bm{P}-\bm{P}^{-T}\|_{F}\leq
2\varepsilon\\\ \bm{P}\in\mathbb{C}^{r\times
r},\text{invertible}\end{subarray}}\left\|\bm{Z}{-}\bm{Z}_{\star}\bm{P}\right\|_{F}^{2}{+}\left\|\bm{Z}{-}\bm{Z}_{\star}\bm{P}^{-T}\right\|_{F}^{2}.$

When $\varepsilon=0$, one has $\bm{P}=\bm{P}^{-T}$ and this means
$\bm{P}\bm{P}^{T}=\bm{I}$ where $\bm{P}$ is a complex orthogonal matrix, then
we obtain

$\displaystyle\mbox{dist}^{2}_{P}(\bm{Z},\bm{Z}_{\star})$
$\displaystyle{=}\inf_{\begin{subarray}{c}\|\bm{P}-\bm{P}^{-T}\|_{F}=0\\\
\bm{P}\in\mathbb{C}^{r\times
r},\text{invertible}\end{subarray}}\left\|\bm{Z}{-}\bm{Z}_{\star}\bm{P}\right\|_{F}^{2}{+}\left\|\bm{Z}{-}\bm{Z}_{\star}\bm{P}^{-T}\right\|_{F}^{2}$
$\displaystyle=2\mbox{dist}^{2}_{Q}(\bm{Z},\bm{Z}_{\star}),$

that is, the distance metric $\mbox{dist}_{P}^{2}(\bm{Z},\bm{Z}_{\star})$ is
equivalent to the distance metric
$\mbox{dist}_{Q}^{2}(\bm{Z},\bm{Z}_{\star})$. According to the linear
convergence of $\mbox{dist}_{P}(\bm{Z},\bm{Z}_{\star})$ in Theorem 1, we have
the linear convergence of $\varepsilon$ to zero. As a result,
$\mbox{dist}_{P}^{2}(\bm{Z},\bm{Z}_{\star})$ reduces to
$\mbox{dist}_{Q}^{2}(\bm{Z},\bm{Z}_{\star})$ asymptotically as $\varepsilon\to
0$.

## Appendix E

###  E-A Proof of Lemma 4

Obviously, one can establish the equivalence between the following two
minimization problems

$\displaystyle\inf_{\begin{subarray}{c}\bm{P}\in\mathbb{C}^{r\times r}\\\
\text{ invertible
}\end{subarray}}\sqrt{\|\bm{Z}-\bm{Z}_{\star}\bm{P}\|_{F}^{2}+\|\bm{Z}-\bm{Z}_{\star}\bm{P}^{-T}\|_{F}^{2}}$
(42) $\displaystyle=\inf_{\begin{subarray}{c}\bm{R}\in\mathbb{C}^{r\times
r}\\\ \text{ invertible
}\end{subarray}}\sqrt{\|\bm{Z}-\bm{Z}_{\star}\hat{\bm{P}}\bm{R}\|_{F}^{2}+\|\bm{Z}-\bm{Z}_{\star}\hat{\bm{P}}^{-T}\bm{R}^{-T}\|_{F}^{2}}.$
(43)

If the minimizer of the second optimization problem (43) is attained at some
$\bm{R}$, then $\hat{\bm{P}}\bm{R}$ must be the minimizer of the first problem
(42).

As $\hat{\bm{Z}}_{\star{1}}=\bm{Z}_{\star}\hat{\bm{P}}$ and
$\hat{\bm{Z}}_{\star{2}}=\bm{Z}_{\star}\hat{\bm{P}}^{-T}$, we have

$\displaystyle\inf_{\begin{subarray}{c}\bm{R}\in\mathbb{C}^{r\times r}\\\
\text{ invertible
}\end{subarray}}\sqrt{\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\bm{R}\|_{F}^{2}+\|\bm{Z}-\hat{\bm{Z}}_{\star{2}}\bm{R}^{-T}\|_{F}^{2}}$
$\displaystyle\leq\sqrt{\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\|_{F}^{2}+\|\bm{Z}-\hat{\bm{Z}}_{\star{2}}\|_{F}^{2}}.$

For any $\bm{R}$ such that
$\sqrt{\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\bm{R}\|_{F}^{2}+\|\bm{Z}-\hat{\bm{Z}}_{\star{2}}\bm{R}^{-T}\|_{F}^{2}}\leq\sqrt{\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\|_{F}^{2}+\|\bm{Z}-\hat{\bm{Z}}_{\star{2}}\|_{F}^{2}}$,
we have

$\displaystyle\|\hat{\bm{Z}}_{\star{1}}(\bm{I}-\bm{R})\|_{F}$
$\displaystyle\leq\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\|_{F}+\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\bm{R}\|_{F}$
$\displaystyle\leq 2\varepsilon\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}},$

where we use the fact
$\sqrt{\|\bm{Z}-\hat{\bm{Z}}_{\star{1}}\|_{F}^{2}+\|\bm{Z}-\hat{\bm{Z}}_{\star{2}}\|_{F}^{2}}\leq\varepsilon\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}$
in the last inequality. When
$\sigma_{r}(\hat{\bm{Z}}_{\star{1}})\geq\frac{\sqrt{2}}{2}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}$,
we further obtain that

$\displaystyle\|\hat{\bm{Z}}_{\star{1}}(\bm{I}-\bm{R})\|_{F}\geq\frac{\sqrt{2}}{2}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}\|\bm{I}-\bm{R}\|_{F}.$

Thus we obtain that $\|\bm{I}-\bm{R}\|_{F}\leq 2\sqrt{2}\varepsilon$. From
Weyl’s inequality we see that $\sigma_{r}(\bm{R})\geq
1-2\sqrt{2}\varepsilon\geq\frac{1}{2}$ as long as
$\varepsilon\leq\frac{\sqrt{2}}{8}$, which denotes $\bm{R}$ is invertible.

Therefore the minimization problem (43) is equivalent to

$\displaystyle\min_{\begin{subarray}{c}\bm{R}\in\mathbb{C}^{r\times r}\\\
\text{ invertible }\end{subarray}}$
$\displaystyle\sqrt{\|\bm{Z}-\bm{Z}_{\star}\hat{\bm{P}}\bm{R}\|_{F}^{2}+\|\bm{Z}-\bm{Z}_{\star}\hat{\bm{P}}^{-T}\bm{R}^{-T}\|_{F}^{2}}$
$\displaystyle\mathrm{s.t.}~{}~{}~{}~{}$
$\displaystyle\|\bm{I}-\bm{R}\|_{F}\leq 2\sqrt{2}\varepsilon,$

which is a continuous optimization problem over a compact set of invertible
matrices. The minimizer exists from Weierstrass extreme value theorem.
Similarly, we can prove the existence of minimizer if we have
$\sigma_{r}(\hat{\bm{Z}}_{\star{2}})\geq\frac{\sqrt{2}}{2}\sigma_{r}(\bm{M}_{\star})^{\frac{1}{2}}$.

###  E-B Proof of Lemma 5

We introduce the following special form of orthogonal Procrustes problem
firstly

$\displaystyle\min_{\bm{Q}\in\mathcal{O}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\right\|_{F}.$

Let the SVD of $\operatorname{{Re}}(\bm{Z}_{\star}^{H}\bm{Z})$ being
$\operatorname{{Re}}(\bm{Z}_{\star}^{H}\bm{Z})=\bm{Q}_{1}\bm{\Lambda}\bm{Q}_{2}^{T}$,
this problem has a closed-form solution as
$\bm{Q}_{\bm{z}}=\bm{Q}_{1}\bm{Q}_{2}^{T}.$

Let $\bm{F}=\begin{bmatrix}\bm{Z}\\\ \bar{\bm{Z}}\end{bmatrix}$,
$\bm{F}_{\star}=\begin{bmatrix}\bm{Z}_{\star}\\\
\bar{\bm{Z}}_{\star}\end{bmatrix}$,
$\bm{\Delta}=\bm{F}\bm{Q}_{\bm{z}}^{T}-\bm{F}_{\star}$, and we know
$\operatorname{{Re}}(\bm{Z}_{\star}^{H}\bm{Z})=\frac{1}{2}\bm{F}_{\star}^{H}\bm{F}=\bm{Q}_{1}\bm{\Lambda}\bm{Q}_{2}^{T}$.
As $\bm{Q}_{\bm{z}}{=\bm{Q}_{1}\bm{Q}_{2}^{T}}$,
$\bm{F}_{\star}^{H}\bm{F}\bm{Q}_{\bm{z}}^{T}$ is a positive semi-definite real
matrix and this means
$\bm{F}_{\star}^{H}\bm{\Delta}=\bm{F}_{\star}^{H}\bm{F}\bm{Q}_{\bm{z}}^{T}-\bm{F}_{\star}^{H}\bm{F}_{\star}=\bm{Q}_{\bm{z}}\bm{F}^{H}\bm{F}_{\star}-\bm{F}_{\star}^{H}\bm{F}_{\star}=\bm{\Delta}^{H}\bm{F}_{\star}$.
One can establish

$\displaystyle\left\|\bm{F}\bm{F}^{H}-\bm{F}_{\star}\bm{F}_{\star}^{H}\right\|_{F}^{2}$
$\displaystyle=\left\|\bm{F}_{\star}\bm{\Delta}^{H}+\bm{\Delta}\bm{F}_{\star}^{H}+\bm{\Delta}\bm{\Delta}^{H}\right\|_{F}^{2}$
$\displaystyle=\operatorname{tr}(2\bm{F}_{\star}^{H}\bm{F}_{\star}\bm{\Delta}^{H}\bm{\Delta}+(\bm{\Delta}^{H}\bm{\Delta})^{2}+(\bm{F}_{\star}^{H}\bm{\Delta})^{2}+(\bm{\Delta}^{H}\bm{F}_{\star})^{2}$
$\displaystyle\quad+2\bm{F}_{\star}^{H}\bm{\Delta}\bm{\Delta}^{H}\bm{\Delta}+2\bm{\Delta}^{H}\bm{F}_{\star}\bm{\Delta}^{H}\bm{\Delta})$
$\displaystyle\overset{(a)}{=}\operatorname{tr}(2\bm{F}_{\star}^{H}\bm{F}_{\star}\bm{\Delta}^{H}\bm{\Delta}+(\bm{\Delta}^{H}\bm{\Delta}+\sqrt{2}\bm{F}_{\star}^{H}\bm{\Delta})^{2}$
$\displaystyle\quad+(4-2\sqrt{2})\bm{F}_{\star}^{H}\bm{\Delta}\bm{\Delta}^{H}\bm{\Delta})$
$\displaystyle=\operatorname{tr}(2(\sqrt{2}-1)\bm{F}_{\star}^{H}\bm{F}_{\star}\bm{\Delta}^{H}\bm{\Delta}+(\bm{\Delta}^{H}\bm{\Delta}+\sqrt{2}\bm{F}_{\star}^{H}\bm{\Delta})^{2}$
$\displaystyle\quad+(4-2\sqrt{2})\bm{F}_{\star}^{H}\bm{F}\bm{Q}_{\bm{z}}^{T}\bm{\Delta}^{H}\bm{\Delta})$
$\displaystyle\overset{(b)}{\geq}4(\sqrt{2}-1)\|(\bm{F}\bm{Q}_{\bm{z}}^{T}-\bm{F}_{\star})\bm{\Sigma}_{\star}^{\frac{1}{2}}\|_{F}^{2}$
$\displaystyle\geq
8(\sqrt{2}-1)\sigma_{r}(\bm{M}_{\star})\underset{\bm{Q}\in\mathcal{O}}{\operatorname{min}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\right\|_{F}^{2}$
$\displaystyle\overset{(c)}{\geq}4(\sqrt{2}-1)\sigma_{r}(\bm{M}_{\star}){\mbox{dist}^{2}_{P}}(\bm{Z},\bm{Z}_{\star}).$
(44)

Step $(a)$ results from
$\bm{F}_{\star}^{H}\bm{\Delta}=\bm{\Delta}^{H}\bm{F}_{\star}$. Step $(b)$
follows from the fact
$\operatorname{tr}(\bm{F}_{\star}^{H}\bm{F}\bm{Q}_{\bm{z}}^{T}\bm{\Delta}^{H}\bm{\Delta})\geq
0$ ($\bm{\Delta}^{H}\bm{\Delta}$ and
$\bm{F}_{\star}^{H}\bm{F}\bm{Q}_{\bm{z}}^{T}$ is positive semi-definite) and
the fact
$(\bm{\Delta}^{H}\bm{\Delta}+\sqrt{2}\bm{F}_{\star}^{H}\bm{\Delta})^{H}=\bm{\Delta}^{H}\bm{\Delta}+\sqrt{2}\bm{F}_{\star}^{H}\bm{\Delta}$.
Step $(c)$ invokes the fact that

$\displaystyle{\mbox{dist}^{2}_{P}}(\bm{Z},\bm{Z}_{\star})$
$\displaystyle=\inf_{\begin{subarray}{c}\bm{P}\in\mathbb{C}^{r\times r}\\\
\text{ invertible
}\end{subarray}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}\right\|_{F}^{2}+\left\|\bm{Z}-\bm{Z}_{\star}\bm{P}^{-T}\right\|_{F}^{2}$
$\displaystyle\leq\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}_{\bm{z}}\right\|_{F}^{2}+\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}_{\bm{z}}^{-T}\right\|_{F}^{2}$
$\displaystyle=2\underset{\bm{Q}\in\mathcal{O}}{\operatorname{min}}\left\|\bm{Z}-\bm{Z}_{\star}\bm{Q}\right\|_{F}^{2}.$

We establish the following inequality from some technical derivations

$\displaystyle\left\|\bm{Z}\bm{Z}^{H}-\bm{Z}_{\star}\bm{Z}_{\star}^{H}\right\|_{F}^{2}+\left\|\bar{\bm{Z}}\bar{\bm{Z}}^{H}-\bar{\bm{Z}}_{\star}\bar{\bm{Z}}_{\star}^{H}\right\|_{F}^{2}$
$\displaystyle\quad-2\left\|\bm{Z}\bm{Z}^{T}-\bm{Z}_{\star}\bm{Z}_{\star}^{T}\right\|_{F}^{2}$
$\displaystyle=\left(\left\|\bm{Z}\bm{Z}^{H}\right\|_{F}^{2}+\left\|\bar{\bm{Z}}\bar{\bm{Z}}^{H}\right\|_{F}^{2}-2\left\|\bm{Z}\bm{Z}^{T}\right\|_{F}^{2}\right)$
$\displaystyle\quad+\left(\left\|\bm{Z}_{\star}\bm{Z}_{\star}^{H}\right\|_{F}^{2}+\left\|\bar{\bm{Z}}_{\star}\bar{\bm{Z}}_{\star}^{H}\right\|_{F}^{2}-2\left\|\bm{Z}_{\star}\bm{Z}_{\star}^{T}\right\|_{F}^{2}\right)$
$\displaystyle\quad-2\left(\left\|\bm{Z}^{H}\bm{Z}_{\star}\right\|_{F}^{2}+\left\|\bar{\bm{Z}}^{H}\bar{\bm{Z}}_{\star}\right\|_{F}^{2}-2\operatorname{{Re}}\left\langle\bm{Z}\bm{Z}^{T},\bm{Z}_{\star}\bm{Z}_{\star}^{T}\right\rangle\right)$
$\displaystyle=-2\left\|\bm{Z}^{H}\bm{Z}_{\star}-\bar{\bm{Z}}^{H}\bar{\bm{Z}}_{\star}\right\|_{F}^{2}\leq
0,$

where the second equality is from
$\left\|\bm{Z}\bm{Z}^{T}\right\|_{F}^{2}=\left\|\bm{Z}\bm{Z}^{H}\right\|_{F}^{2}=\left\|\bar{\bm{Z}}\bar{\bm{Z}}^{H}\right\|_{F}^{2}=\left\|\bm{\Sigma}\right\|_{F}^{2}$
and
$\left\|\bm{Z}_{\star}\bm{Z}_{\star}^{T}\right\|_{F}^{2}=\left\|\bm{Z}_{\star}\bm{Z}_{\star}^{H}\right\|_{F}^{2}=\left\|\bar{\bm{Z}}_{\star}\bar{\bm{Z}}_{\star}^{H}\right\|_{F}^{2}=\left\|\bm{\Sigma}_{\star}\right\|_{F}^{2}$.
The previous inequality helps establish the upper-bound of
$\left\|\bm{F}\bm{F}^{H}-\bm{F}_{\star}\bm{F}_{\star}^{H}\right\|_{F}^{2}$,
which is

$\displaystyle\left\|\bm{F}\bm{F}^{H}-\bm{F}_{\star}\bm{F}_{\star}^{H}\right\|_{F}^{2}$
$\displaystyle=2\left\|\bm{Z}\bm{Z}^{T}-\bm{Z}_{\star}\bm{Z}_{\star}^{T}\right\|_{F}^{2}+\left\|\bm{Z}\bm{Z}^{H}-\bm{Z}_{\star}\bm{Z}_{\star}^{H}\right\|_{F}^{2}$
$\displaystyle\quad+\left\|\bar{\bm{Z}}\bar{\bm{Z}}^{H}-\bar{\bm{Z}}_{\star}\bar{\bm{Z}}_{\star}^{H}\right\|_{F}^{2}$
$\displaystyle\leq
4\left\|\bm{Z}\bm{Z}^{T}-\bm{Z}_{\star}\bm{Z}_{\star}^{T}\right\|_{F}^{2}=4\left\|\bm{M}-\bm{M}_{\star}\right\|_{F}^{2}.$
(45)

combining (44) and (45), we complete the proof of Lemma 5.

###  E-C Proof of Lemma 8

Let $\bm{H}_{a}$ denotes the $n_{s}\times n_{s}$ Hankel basis matrix with
entries in the $a$-th skew diagonal as $\frac{1}{\sqrt{w_{a}}}$. Since
$p=\frac{m}{n}$, we have

$\displaystyle\left\langle\mathcal{G}(p^{-1}\mathcal{P}_{{\hat{\Omega}}}-\mathcal{I})\mathcal{G}^{*}(\bm{A}\bm{B}^{T}),\bm{C}\bm{D}^{T}\right\rangle$
$\displaystyle=\sum_{k=1}^{m}\Big{(}\frac{n}{m}\langle\bm{A}\bm{B}^{T},\bm{H}_{a_{k}}\rangle\langle\bm{H}_{a_{k}},\bm{C}\bm{D}^{T}\rangle$
$\displaystyle\quad-\frac{1}{m}\sum_{a=0}^{n-1}\langle\bm{A}\bm{B}^{T},\bm{H}_{a}\rangle\langle\bm{H}_{a},\bm{C}\bm{D}^{T}\rangle\Big{)}=\sum_{k=1}^{m}z_{a_{k}},$

It is obviously that $z_{a_{k}}$ are i.i.d random variables and
$\mathbb{E}[z_{a_{k}}]=0$. We first point out some results

$\displaystyle|\langle\bm{A}\bm{B}^{T},\bm{H}_{a}\rangle|$
$\displaystyle=|\langle\bm{A},\bm{H}_{a}\bar{\bm{B}}\rangle|\leq\|\bm{A}\|_{F}\|\bm{B}\|_{2,\infty},$
$\displaystyle|\langle\bm{A}\bm{B}^{T},\bm{H}_{a}\rangle|$
$\displaystyle=|\langle\bm{B}^{T},\bm{A}^{H}\bm{H}_{a}\rangle|\leq\|\bm{A}\|_{2,\infty}\|\bm{B}\|_{F}.$

where we use the fact that
$\|\bm{H}_{a}\bar{\bm{B}}\|_{F}\leq\|\bm{B}\|_{2,\infty}$ and
$\|\bm{A}^{H}\bm{H}_{a}\|_{F}\leq\|\bm{A}\|_{2,\infty}$. Similar results hold
for $|\langle\bm{H}_{a},\bm{C}\bm{D}^{T}\rangle|$. Thus we have

$\displaystyle|z_{a_{k}}|\leq\frac{2n}{m}\max_{a}|\langle\bm{A}\bm{B}^{T},\bm{H}_{a}\rangle\langle\bm{H}_{a},\bm{C}\bm{D}^{T}\rangle|\leq\frac{2n}{m}L,$

where

$\displaystyle L=$
$\displaystyle\min\\{\left\|\bm{A}\right\|_{F}\left\|\bm{B}\right\|_{2,\infty},\left\|\bm{A}\right\|_{2,\infty}\left\|\bm{B}\right\|_{F}\\}$
$\displaystyle\cdot\min\\{\left\|\bm{C}\right\|_{F}\left\|\bm{D}\right\|_{2,\infty},\left\|\bm{C}\right\|_{2,\infty}\left\|\bm{D}\right\|_{F}\\}.$

As $z_{a_{k}}$ is a complex scalar, we have

$\displaystyle\mathbb{E}\left({z_{a_{k}}^{H}}z_{a_{k}}\right)=\mathbb{E}(z_{a_{k}}z_{a_{k}}^{H})$
$\displaystyle=\mathbb{E}(\frac{n^{2}}{m^{2}}|\langle\bm{A}\bm{B}^{T},\bm{H}_{a_{k}}\rangle|^{2}|\langle\bm{H}_{a_{k}},\bm{C}\bm{D}^{T}\rangle|^{2})$
$\displaystyle\quad-\frac{1}{m^{2}}|\sum_{a=0}^{n-1}\langle\bm{A}\bm{B}^{T},\bm{H}_{a}\rangle\langle\bm{H}_{a},\bm{C}\bm{D}^{T}\rangle|^{2}$
$\displaystyle\leq\mathbb{E}(\frac{n^{2}}{m^{2}}|\langle\bm{A}\bm{B}^{T},\bm{H}_{a_{k}}\rangle|^{2}|\langle\bm{H}_{a_{k}},\bm{C}\bm{D}^{T}\rangle|^{2}).$

Therefore

$\displaystyle\sum_{k=1}^{m}\mathbb{E}\left(z_{a_{k}}z_{a_{k}}^{H}\right)=\sum_{k=1}^{m}\mathbb{E}\left(z_{a_{k}}^{H}z_{a_{k}}\right)\leq\frac{n^{2}}{m}L^{2}.$

Invoking Bernstein’s inequality [47, Theorem 1.6], we have

$\displaystyle\mathbb{P}\left(\left|\sum_{k=1}^{m}z_{a_{k}}\right|\geq
t\right)\leq 2\exp\left(\frac{-mt^{2}/2}{n^{2}L^{2}+2nLt/3}\right).$

Consequently when $m\gtrsim O(\log(n))$,
$\left|\sum_{k=1}^{m}z_{a_{k}}\right|\leq c_{4}L\sqrt{\frac{n^{2}\log(n)}{m}}$
holds with probability at least $1-n^{-2}$.

## References

  * [1] L. Potter, E. Ertin, J. Parker, and M. Cetin, “Sparsity and compressed sensing in radar imaging,” Proc. IEEE, vol. 98, no. 6, pp. 1006–1020, 2010.
  * [2] J. A. Tropp, J. N. Laska, M. F. Duarte, J. K. Romberg, and R. G. Baraniuk, “Beyond nyquist: Efficient sampling of sparse bandlimited signals,” IEEE Trans. Inf. Theory, vol. 56, no. 1, pp. 520–544, 2010.
  * [3] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med., vol. 58, no. 6, pp. 1182–1195, 2007.
  * [4] J. L. Paredes, G. R. Arce, and Z. Wang, “Ultra-wideband compressed sensing: Channel estimation,” IEEE J. Sel. Top. Signal Process., vol. 1, no. 3, pp. 383–395, 2007.
  * [5] Y. Chen and Y. Chi, “Robust spectral compressed sensing via structured matrix completion,” IEEE Trans. Inf. Theory, vol. 60, no. 10, pp. 6576–6601, 2014\.
  * [6] J.-F. Cai, T. Wang, and K. Wei, “Spectral compressed sensing via projected gradient descent,” SIAM J. Optim., vol. 28, no. 3, pp. 2625–2653, 2018\.
  * [7] J.-F. Cai, T. Wang, and K. Wei, “Fast and provable algorithms for spectrally sparse signal reconstruction via low-rank hankel matrix completion,” Appl. Comput. Harmon. Anal., vol. 46, no. 1, pp. 94–121, 2019.
  * [8] X. Zhang, Y. Liu, and W. Cui, “Spectrally sparse signal recovery via hankel matrix completion with prior information,” IEEE Trans. Signal Process., vol. 69, pp. 2174–2187, 2021.
  * [9] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressed sensing off the grid,” IEEE Trans. Inf. Theory, vol. 59, no. 11, pp. 7465–7490, 2013.
  * [10] R. A. Horn and C. R. Johnson, Matrix Analysis.  USA: Cambridge University Press, 2nd ed., 2012.
  * [11] R. A. Horn and D. I. Merino, “The jordan canonical forms of complex orthogonal and skew-symmetric matrices,” Linear Algebra Appl., vol. 302-303, pp. 411–421, 1999.
  * [12] F. T. Luk and S. Qiao, “Using complex-orthogonal transformations to diagonalize a complex symmetric matrix,” in Opt. Photon., 1997.
  * [13] V. D. Didenko and V. A. Chernetskii, “The riemann boundary problem with a complex orthogonal matrix,” Math. Notes Acad. Sci. USSR, vol. 23, pp. 220–227, Mar. 1978.
  * [14] E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489 – 509, 2006.
  * [15] D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, no. 4, pp. 1289 –1306, 2006.
  * [16] J. Shen, J.-S. Chen, H.-D. Qi, and N. Xiu, “A penalized method of alternating projections for weighted low-rank hankel matrix optimization,” Math. Program. Comput., vol. 14, no. 3, pp. 417–450, 2022.
  * [17] Q. Zheng and J. Lafferty, “Convergence analysis for rectangular matrix completion using Burer-Monteiro factorization and gradient descent,” arXiv:1605.07051, 2016.
  * [18] W. Xu and S. Qiao, “A fast symmetric SVD algorithm for square hankel matrices,” Linear Algebra Appl., vol. 428, no. 2-3, pp. 550–563, 2008\.
  * [19] F. Andersson, M. Carlsson, and P.-A. Ivert, “A fast alternating projection method for complex frequency estimation,” arXiv:1107.2028, 2011.
  * [20] K. Wei, J.-F. Cai, T. F. Chan, and S. Leung, “Guarantees of riemannian optimization for low rank matrix recovery,” SIAM J. Matrix Anal. Appl., vol. 37, no. 3, pp. 1198–1222, 2016.
  * [21] K. Wei, J.-F. Cai, T. F. Chan, and S. Leung, “Guarantees of riemannian optimization for low rank matrix completion,” Inverse Probl. Imag., vol. 14, no. 2, pp. 233–265, 2020.
  * [22] X. Wu, Z. Yang, J.-F. Cai, and Z. Xu, “Spectral super-resolution on the unit circle via gradient descent,” in IEEE Int. Conf. Acoust., Speech, Signal Process., 2023.
  * [23] S. Mao and J. Chen, “Blind super-resolution of point sources via projected gradient descent,” IEEE Trans. Signal Process., vol. 70, pp. 4649–4664, 2022.
  * [24] J. Yin, S. Peng, Z. Yang, B. Chen, and Z. Lin, “Hypergraph based semi-supervised symmetric nonnegative matrix factorization for image clustering,” Pattern Recogn., vol. 137, p. 109274, May 2023.
  * [25] Z. He, S. Xie, R. Zdunek, G. Zhou, and A. Cichocki, “Symmetric nonnegative matrix factorization: Algorithms and applications to probabilistic clustering,” IEEE Trans. Neural Netw., vol. 22, no. 12, pp. 2117–2131, 2011.
  * [26] Y. Chen and M. J. Wainwright, “Fast low-rank estimation by projected gradient descent: General statistical and algorithmic guarantees,” arxiv:1509.03025, 2015.
  * [27] P. Comon, G. Golub, L.-H. Lim, and B. Mourrain, “Symmetric tensors and symmetric tensor rank,” SIAM J. Matrix Anal. Appl, vol. 30, pp. 1254–1279, Jan. 2008.
  * [28] J. Brachat, P. Comon, B. Mourrain, and E. Tsigaridas, “Symmetric tensor decomposition,” Linear Algebra Appl., vol. 433, pp. 1851–1872, Dec. 2010\.
  * [29] J. Cai, M. Baskaran, B. Meister, and R. Lethin, “Optimization of symmetric tensor computations,” in IEEE High Perform. Extreme Comput. Conf., pp. 1–7, 2015.
  * [30] E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Found. Comput. Math., vol. 9, no. 6, pp. 717–772, 2009\.
  * [31] W. Liao and A. Fannjiang, “MUSIC for single-snapshot spectral estimation: Stability and super-resolution,” Appl. Comput. Harmon. Anal., vol. 40, no. 1, pp. 33–67, 2016.
  * [32] Y. Cherapanamjeri, K. Gupta, and P. Jain, “Nearly optimal robust matrix completion,” in Proc. Int. Conf. Mach. Learn., pp. 797–805, 2017.
  * [33] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix completion using alternating minimization,” in Proc. 45th Annu. ACM symp. Theory Comput., pp. 665–674, ACM, 2013.
  * [34] S. Zhang, Y. Hao, M. Wang, and J. H. Chow, “Multichannel hankel matrix completion through nonconvex optimization,” IEEE J. Sel. Top. Signal Process., vol. 12, no. 4, pp. 617–632, 2018.
  * [35] S. Zhang and M. Wang, “Correction of corrupted columns through fast robust hankel matrix completion,” IEEE Trans. Signal Process., vol. 67, no. 10, pp. 2580–2594, 2019.
  * [36] A. M. Chebotarev and A. E. Teretenkov, “Singular value decomposition for the takagi factorization of symmetric matrices,” Appl. Math. Comput., vol. 234, pp. 380–384, 2014.
  * [37] K. D. Ikramov, “Takagi’s decomposition of a symmetric unitary matrix as a finite algorithm,” Comp. Math. Math. Phys., vol. 52, no. 1, pp. 1–3, 2012\.
  * [38] S. Tu, R. Boczar, M. Simchowitz, M. Soltanolkotabi, and B. Recht, “Low-rank solutions of linear matrix equations via procrustes flow,” in Proc. Int. Conf. Mach. Learn., pp. 964–973, 2016.
  * [39] C. Ma, K. Wang, Y. Chi, and Y. Chen, “Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion, and blind deconvolution,” Found. Comput. Math., vol. 20, no. 3, pp. 451–632, 2019.
  * [40] R. Vershynin, High-Dimensional Probability.  Cambridge University Press, 2018.
  * [41] J. B. Sanson, P. M. Tome, D. Castanheira, A. Gameiro, and P. P. Monteiro, “High-resolution delay-doppler estimation using received communication signals for OFDM radar-communication system,” IEEE Trans. Veh. Technol., vol. 69, pp. 13112–13123, nov 2020.
  * [42] H. Nikookar and R. Prasad, “Multicarrier transmission with nonuniform carriers in a multipath channel,” in Proc. Int. Conf. Univ. Pers. Commun., IEEE, 1996.
  * [43] Q. Yang, Z. Wang, and Q. Huang, “An efficient non-uniform multi-tone system based on ramanujan sums,” in Proc. 2nd Joint Int. Inf. Tech. Mech. Electron. Eng. Conf., Atlantis Press, 2017.
  * [44] C. R. Berger, B. Demissie, J. Heckenbach, P. Willett, and S. Zhou, “Signal processing for passive radar using ofdm waveforms,” IEEE J. Sel. Top. Signal Process., vol. 4, no. 1, pp. 226–238, 2010.
  * [45] C. Ma, Y. Li, and Y. Chi, “Beyond procrustes: Balancing-free gradient descent for asymmetric low-rank matrix sensing,” IEEE Trans. Signal Process., vol. 69, pp. 867–877, 2021.
  * [46] B. Recht, “A simpler approach to matrix completion,” J. Mach. Learn. Res., vol. 12, p. 3413–3430, 2011.
  * [47] J. A. Tropp, “User-friendly tail bounds for sums of random matrices,” Found. Comput. Math., vol. 12, no. 4, pp. 389–434, 2012.