Provable Low Rank Phase Retrieval

TL;DR

This paper introduces a provably correct and efficient algorithm, AltMinLowRaP, for solving the Low Rank Phase Retrieval problem, enabling accurate recovery of low-rank matrices from phaseless measurements with theoretical guarantees.

Contribution

The paper presents the first provably correct and fast algorithm for Low Rank Phase Retrieval, with rigorous guarantees and applicability to dynamic settings.

Findings

Algorithm achieves $ ext{epsilon}$-accuracy with high probability.

Recovery guaranteed when $mq \, \geq C n r^4 \log(1/\epsilon)$.

Time complexity is $C m q n r \log^2(1/\epsilon)$.

Abstract

We study the Low Rank Phase Retrieval (LRPR) problem defined as follows: recover an $n \times q$ matrix $X^*$ of rank $r$ from a different and independent set of $m$ phaseless (magnitude-only) linear projections of each of its columns. To be precise, we need to recover $X^*$ from $y_k := |A_k{}' x^*_k|, k=1,2,\dots, q$ when the measurement matrices $A_k$ are mutually independent. Here $y_k$ is an $m$ length vector, $A_k$ is an $n \times m$ matrix, and $'$ denotes matrix transpose. The question is when can we solve LRPR with $m \ll n$? A reliable solution can enable fast and low-cost phaseless dynamic imaging, e.g., Fourier ptychographic imaging of live biological specimens. In this work, we develop the first provably correct approach for solving this LRPR problem. Our proposed algorithm, Alternating Minimization for Low-Rank Phase Retrieval (AltMinLowRaP), is an AltMin based solution and hence is also provably fast (converges geometrically). Our guarantee shows that AltMinLowRaP solves LRPR to $\epsilon$ accuracy, with high probability, as long as $m q \ge C n r^4 \log(1/\epsilon)$, the matrices $A_k$ contain i.i.d. standard Gaussian entries, and the right singular vectors of $X^*$ satisfy the incoherence assumption from matrix completion literature. Here $C$ is a numerical constant that only depends on the condition number of $X^*$ and on its incoherence parameter. Its time complexity is only $ C mq nr \log^2(1/\epsilon)$. Since even the linear (with phase) version of the above problem is not fully solved, the above result is also the first complete solution and guarantee for the linear case. Finally, we also develop a simple extension of our results for the dynamic LRPR setting.

Tables

Table 1. TABLE I : Comparing our result (first result for LRPR) with the first and the best results for (bounded time) non-convex algorithms for the three related problems to ours - LRMC, PR, Sparse PR. Here “best” refers to best sample complexity. We treat κ , μ 𝜅 𝜇 \kappa,\mu as constants. 𝑿 ∗ superscript 𝑿 {\bm{X}^{*}} is an n × q 𝑛 𝑞 n\times q rank- r 𝑟 r matrix; 𝒙 ∗ superscript 𝒙 \bm{x}^{*} is an n 𝑛 n -length vector signal; and m 𝑚 m is the number of samples needed per signal (per column of 𝑿 ∗ superscript 𝑿 {\bm{X}^{*}} ). LRPR, Sparse PR, and Standard PR results assume iid Gaussian measurements, while LRMC results assume iid Bernoulli model on observed entries.

Problem	Global	Assumptions	Sample Complexity	Time Complexity per signal
	Measurements?			(with $m$ = its lower bound)
LRPR (first)	No	right incoherence,	$m\ge C\frac{n}{q}{r}^4\log \left(1/ϵ\right)$,	$Cn\left(n/q\right)r^5{\log }^3\left(1/ϵ\right)$
(our work)		${𝑿}^{\ast }$ has rank $r$	$m\ge C\max (r,\log n,\log q)$
LRMC (first) [25]	No	left & right incoherence,	$m\ge C\frac{\max (n,q)}{q}{r}^{4.5}\log \left(1/ϵ\right)$	$C\left(n/q\right)r^{6.5}\log n{\log }^2\left(1/ϵ\right)$
		${𝑿}^{\ast }$ has rank $r$
LRMC (best) [26]	No	left & right incoherence	$m\ge C\frac{\max (n,q)}{q}{r}^2{\log }^{2}n{\log }^2\left(1/ϵ\right)$	$C\left(n/q\right)r^3{\log }^n{\log }^2\left(1/ϵ\right)$
		${𝑿}^{\ast }$ has rank $r$
Sparse PR (first) [6]	Yes	${𝒙}^{\ast }$ is $s$-sparse in canonical basis,	$m\ge C{s}^2\log n\log \left(1/ϵ\right)$	$Cn{s}^3\log n{\log }^2\left(1/ϵ\right)$
		min nonzero entry lower bounded
Sparse PR (best) [17, 18]	Yes	${𝒙}^{\ast }$ is $s$-sparse in canonical basis	$m\ge C{s}^2\log n$	$Cn{s}^2\log n\log \left(1/ϵ\right)$
Standard PR (first) [6]	Yes	None	$m\ge Cn{\log }^{3}n\log \left(1/ϵ\right)$	$C{n}^2{\log }^{3}n\log \left(1/ϵ\right)$
Standard PR (best) [9, 10]	Yes	None	$m\ge Cn$	$C{n}^2\log \left(1/ϵ\right)$

Table 2. TABLE II : mat-dist ​ ( 𝑿 ^ , 𝑿 ∗ ) mat-dist ^ 𝑿 superscript 𝑿 \text{mat-dist}(\hat{\bm{X}},{\bm{X}^{*}}) and time comparison for the mouse and plane videos. We generate the measurements using the CDP model and consider two different number of settings. Notice that AltMinLowRaP is slightly better than LRPR2 but is slower than in the simulated data experiments.

Algorithm	$\text{mat-dist}(\widehat{𝑿},{𝑿}^{\ast })$ (Running Time in sec)
	m = $5n$ (mouse)	m = $5n$ (plane)	m = $2n$ (mouse)
RWF [10]	$0.65$ ($0.35s$)	$0.65$ ($0.35$s)	$1.36$ ($0.25$s)
LRPR2 [20] ($r=25$)	$0.48$ ($81.8$s)	$0.10$ ($122.3$s)	$0.61$ ($31.0$s)
AltMinLowRaP ($r=25$)	$0.39$ ($297.6$s)	$0.09$ ($467.6$s)	$0.52$ ($122.1$s)
AltMinLowRaP ($r=15$)	$0.57$ ($277.2$s)	$0.15$ ($418.9$s)	$0.60$ ($113.6$s)
AltMinLowRaP ($r=10$)	$0.70$ ($262.7$s)	$0.23$ ($409.0$s)	$0.72$ ($97.4$s)

Table 3. TABLE III : Low rank versus Sparse PR: We compare with a recent state of the art algorithm for provable sparse phase retrieval. mat-dist ​ ( 𝑿 ^ , 𝑿 ∗ ) mat-dist ^ 𝑿 superscript 𝑿 \text{mat-dist}(\hat{\bm{X}},{\bm{X}^{*}}) and time comparison for the mouse and plane videos using real valued CDP measurements. We conclude that low-rank is a better model than wavelet-sparse for slowly changing videos.

Algorithm	Video
	$m=5n$ (plane)	$m=10n$ (plane)	$m=5n$ (sparsified-plane)
CoPRAM [18]	$2.113$	$1.019$	$0.3104$
RWF	$0.6531$	$0.4134$	$0.6514$
AltMinLowRaP ($r=15$)	$0.111$	$0.109$	$0.1427$

Equations

$${\mathbf{y}}_{ik}:=|⟨{\mathbf{a}}_{ik}{\mathbf{e}}_k{}^{\prime },{\mathbf{X}}^{\ast }⟩|=|⟨{\mathbf{a}}_{ik},{\mathbf{x}}_k^{\ast }⟩|,i\in [m],k\in [q],$$

$${\mathbf{y}}_k:=|{\mathbf{A}}_k{}^{\prime }{\mathbf{x}}_k^{\ast }|,k=1,2,\dots ,q$$

$${\tilde{\mathbf{B}}}^{\ast }:={\mathbf{\Sigma }}^{\ast }{\mathbf{B}}^{\ast }.$$

$$\underset{k}{\max }\Vert {\mathbf{b}}_k^{\ast }{\Vert }^2\le {\mu }^2\frac{r}{q},$$

$$\Vert {\mathbf{x}}_k^{\ast }{\Vert }^2=\Vert {\tilde{\mathbf{b}}}_k^{\ast }{\Vert }^2\le {{\sigma }_{\max }^{\ast }}^2{\mu }^2\frac{r}{q}={\kappa }^2{{\sigma }_{\min }^{\ast }}^2{\mu }^2\frac{r}{q}\le {\kappa }^2{\mu }^2\frac{\Vert {\mathbf{X}}^{\ast }{\Vert }_F^2}{q}$$

$$\sin \Theta ({\mathbf{U}}_1,{\mathbf{U}}_2):=\Vert (\mathbf{I}-{\mathbf{U}}_1{\mathbf{U}}_1^{\prime }){\mathbf{U}}_2\Vert .$$

$$\mathrm{dist}({\mathbf{x}}^{\ast },\hat{\mathbf{x}}):=\underset{\theta \in [-\pi ,\pi ]}{\min }\Vert {\mathbf{x}}^{\ast }-e^{-j\theta }\hat{\mathbf{x}}\Vert$$

$$\text{mat-dist}(\hat{\mathbf{X}},{\mathbf{X}}^{\ast }{)}^2:=\sum _{k=1}^q{\mathrm{dist}}^2({\mathbf{x}}_k^{\ast },{\hat{\mathbf{x}}}_k).$$

$$\sum _{k=1}^q\Vert {\mathbf{y}}_k-|{\mathbf{A}}_k{}^{\prime }\mathbf{U}{\mathbf{b}}_k|{\Vert }^2$$

$${\mathbf{Y}}_U=\frac{1}{mq}\sum _{k=1}^q\sum _{i=1}^m{\mathbf{y}}_{ik}^2{\mathbf{a}}_{ik}{\mathbf{a}}_{ik}^{\prime }\text{mathbbm}{1}_{\left\{{\mathbf{y}}_{ik}^2\le C_Y\frac{1}{mq}{\sum }_{ik}{\mathbf{y}}_{ik}^2\right\}}.$$

$$\sin \Theta ({\mathbf{U}}^{\ast },{\mathbf{U}}^T)\le ϵ,\text{mat-dist}({\hat{\mathbf{X}}}^T,{\mathbf{X}}^{\ast })\le ϵ\Vert {\mathbf{X}}^{\ast }{\Vert }_F$$

$$\sin \Theta ({\mathbf{U}}^{\ast },{\mathbf{U}}^t)\le 0.7^t{\delta }_{\mathrm{init}},t=0,1,2,\dots ,T$$

$${\mathbf{y}}_{ik}:=⟨{\mathbf{a}}_{ik},{\mathbf{x}}_k^{\ast }⟩,i=1,2,\dots ,m,k=1,2,\dots ,q.$$

$$\sin \Theta ({\mathbf{U}}_{\mathrm{init}},{\mathbf{U}}^{\ast })\le {\delta }_{\mathrm{init}}.$$

$${\mathbf{g}}_k^t:=({\mathbf{U}}^t{)}^{\prime }{\mathbf{x}}_k^{\ast }\text{and}{\mathbf{e}}_k^t:=(\mathbf{I}-{\mathbf{U}}^t{{\mathbf{U}}^t}^{\prime }){\mathbf{x}}_k^{\ast }.$$

$${\mathbf{g}}_k^t=({\mathbf{U}}^t{)}^{\prime }{\mathbf{x}}_k^{\ast }=(({\mathbf{U}}^t{)}^{\prime }{\mathbf{U}}^{\ast }){\tilde{\mathbf{b}}}_k^{\ast }$$

$$\mathrm{dist}\left({\mathbf{g}}_k^t,{\hat{\mathbf{b}}}_k^t\right)$$

$$\text{mat-dist}({\mathbf{G}}^t,{\hat{\mathbf{B}}}^t)$$

$$\mathrm{dist}({\hat{\mathbf{x}}}_k^t,{\mathbf{x}}_k^{\ast })$$

$$\sin \Theta ({\mathbf{U}}^{t+1},{\mathbf{U}}^{\ast })\le 0.7{\delta }_t:={\delta }_{t+1}.$$

$${\mathbf{Y}}_U=\frac{1}{mq}\sum _{ik}|{\mathbf{a}}_{ik}{}^{\prime }{\mathbf{x}}_k^{\ast }{|}^2{\mathbf{a}}_{ik}{\mathbf{a}}_{ik}^{\prime }\text{mathbbm}{1}_{\left\{|{\mathbf{a}}_{ik}{}^{\prime }{\mathbf{x}}_k^{\ast }{|}^2\le \frac{9{\mu }^2{\kappa }^2}{mq}{\sum }_{ik}|{\mathbf{a}}_{ik}{}^{\prime }{\mathbf{x}}_k^{\ast }{|}^2\right\}}$$

$${\mathbf{Y}}_{-}({ϵ}_1)=$$

$$\frac{1}{mq}\sum _{ik}|{\mathbf{a}}_{ik}{}^{\prime }{\mathbf{x}}_k^{\ast }{|}^2{\mathbf{a}}_{ik}{\mathbf{a}}_{ik}^{\prime }\text{mathbbm}{1}_{\left\{|{\mathbf{a}}_{ik}{}^{\prime }{\mathbf{x}}_k^{\ast }{|}^2\le \frac{9{\mu }^2{\kappa }^2(1-{ϵ}_1)}{q}\Vert {\mathbf{X}}^{\ast }{\Vert }_F^2\right\}}.$$

$$\mathbb{E}\left[{\mathbf{Y}}_{-}({ϵ}_1)\right]=\frac{1}{q}\left\{\sum _k{\beta }_{1,k}^{-}{\mathbf{x}}_k^{\ast }{{\mathbf{x}}_k^{\ast }}^{\prime }+\left(\sum _k{\beta }_{2,k}^{-}\Vert {\mathbf{x}}_k^{\ast }{\Vert }^2\right)\mathbf{I}\right\}$$

$${\beta }_{1,k}^{-}({ϵ}_1):=\mathbb{E}[({\xi }^4-{\xi }^2)\text{mathbbm}{1}_{{\xi }^2\le (1-{ϵ}_1){\gamma }_k}]$$

$${\beta }_{2,k}^{-}({ϵ}_1):=\mathbb{E}[{\xi }^2\text{mathbbm}{1}_{{\xi }^2\le (1-{ϵ}_1){\gamma }_k}],$$

$${\gamma }_k:=\frac{9{\mu }^2{\kappa }^2\Vert {\mathbf{X}}^{\ast }{\Vert }_F^2}{q\Vert {\mathbf{x}}_k^{\ast }{\Vert }^2},$$

$$\mathbb{E}\left[{\mathbf{Y}}_{-}({ϵ}_1)\right]=\frac{1}{q}[{\mathbf{U}}^{\ast }(\sum _k{\beta }_{1,k}^{-}{\tilde{\mathbf{b}}}_k^{\ast }{\tilde{\mathbf{b}}}_k^{\ast }{}^{\prime }){\mathbf{U}}^{\ast }{}^{\prime }+(\sum _k{\beta }_{2,k}^{-}\Vert {\tilde{\mathbf{b}}}_k^{\ast }{\Vert }^2)\mathbf{I}]$$

$$\sin \Theta (\mathbf{U},{\mathbf{U}}^{\ast })\le \frac{\Vert \mathbf{D}-\hat{\mathbf{D}}\Vert }{{\lambda }_r(\mathbf{D})-{\lambda }_{r+1}(\mathbf{D})-\Vert \mathbf{D}-\hat{\mathbf{D}}\Vert }.$$

$$\sin \Theta ({\mathbf{U}}_{\mathrm{init}},{\mathbf{U}}^{\ast })$$

$$\le \frac{\Vert {\mathbf{Y}}_U-\mathbb{E}\left[{\mathbf{Y}}_{-}({ϵ}_1)\right]\Vert }{{\lambda }_r(\mathbb{E}\left[{\mathbf{Y}}_{-}({ϵ}_1)\right])-{\lambda }_{r+1}(\mathbb{E}\left[{\mathbf{Y}}_{-}({ϵ}_1)\right])-\Vert {\mathbf{Y}}_U-\mathbb{E}\left[{\mathbf{Y}}_{-}({ϵ}_1)\right]\Vert }$$

$${\lambda }_r(\mathbb{E}\left[{\mathbf{Y}}_{-}\right])-{\lambda }_{r+1}(\mathbb{E}\left[{\mathbf{Y}}_{-}\right])$$

$$\ge (\underset{k}{\min }{\beta }_{1,k}^{-})\frac{{{\sigma }_{\min }^{\ast }}^2}{q}$$

$${\mathbf{Y}}_{-}({ϵ}_1)⪯{\mathbf{Y}}_U⪯{\mathbf{Y}}_{+}({ϵ}_1)$$

Full text

Since even the linear (with phase) version of the above problem is not fully solved, the above result is also the first complete solution and guarantee for the linear case. Finally, we also develop a simple extension of our results for the dynamic LRPR setting.

I Introduction

In recent years, there has been a resurgence of interest in the classical phase retrieval (PR) problem [2, 3]. The original PR problem involved recovering an $n$-length signal $\bm{x}^{*}$ from the magnitudes of its Discrete Fourier Transform (DFT) coefficients. Its generalized version, studied in recent literature, replaces DFT by inner products with any arbitrary design vectors, $\bm{a}_{i}$. Thus, the goal is to recover $\bm{x}^{*}$ from $\bm{y}_{i}:=|\langle\bm{a}_{i},\bm{x}^{*}\rangle|$, $i=1,2,\dots,m$. These are commonly referred to as phaseless linear projections of $\bm{x}^{*}$. While practical PR methods have existed for a long time, e.g., see [2, 3], the focus of the recent work has been on obtaining correctness guarantees for these and newer algorithms. This line of work includes convex relaxation methods [4, 5] as well as non-convex methods [6, 7, 8, 9, 10, 11, 12, 13]. It is easy to see that, without extra assumptions, PR requires $m\geq n$. The best known guarantees – see [9] and follow-up works – prove exact recovery with high probability (whp) with order-optimal number of measurements/samples: $m=Cn$; and with time complexity $Cmn\log(1/\epsilon)$ that is nearly linear in the problem size. Here and below, $C$ is reused often to refer to a constant more than one. Most guarantees for PR assume that $\bm{a}_{i}$’s are independent and identically distributed (iid) standard Gaussian vectors. When this is assumed, we refer to the PR problem as “standard PR”.

A natural approach to reduce the sample complexity to below $n$ is to impose structure on the unknown signal(s). In existing literature, with the exception of sparse PR which has been extensively studied, e.g., [14, 15, 6, 16, 17, 18], there is little other work on structured PR. Low rank is the other common structure. This can be used in one of two ways. One is to assume that the unknown signal/image, whose phaseless linear projections are available, can be rearranged to form a low-rank matrix. This would be valid only for very specific types of images for which different image rows or columns look similar, so that the entire image matrix can be modeled as low rank. In general it is not a very practical model for images, and this is probably why this setting has not been explored in the literature. We do not consider this model here either.

A more practical, and commonly used, low-rank model in biological applications [19], is for the dynamic imaging setting. It assumes that a set, e.g., a time sequence, of signals/images is generated from a lower dimensional subspace of the ambient space. For our problem, we assume that we have a set of $m$ phaseless linear projections of each signal, with a different set of measurement vectors used for each signal. The question is can we jointly recover the signals using an $m\ll n$ and when? This setting was first studied in our recent work [20] where we called it “Low-Rank PR” (LRPR). It is a valid model whenever the set/sequence of signals is sufficiently similar (correlated). A solution to LRPR can enable fast and low-cost phaseless dynamic imaging of live biological specimens, in vitro. See Sec. I-C and [21] for a detailed motivation for studying LRPR.

I-A Low Rank PR (LRPR) Problem Setting and Notation

We study the LRPR problem described above. This was first introduced and briefly studied in [20] where we developed two algorithms, evaluated them experimentally, and provided a guarantee for the initialization step of one of them. The goal is to to recover an $n\times q$ matrix ${\bm{X}^{*}}:=[\bm{x}^{*}_{1},\bm{x}^{*}_{2},\dots,\bm{x}^{*}_{k},\dots,\bm{x}^{*}_{q}]$, of rank $r$, from measurements

[TABLE]

when all the $\bm{a}_{ik}$’s are mutually independent. For proving guarantees we assume also that they are iid standard Gaussian and real-valued. Here and below the notation $[m]:=\{1,2,\dots,m\}$, and $\bm{e}_{k}$ refers to the $k$-th column of $\bm{I}_{q}$ (identity matrix of size $q\times q$). We are interested in the low rank setting when $r\ll\min(n,q)$.

By defining the $m$-length vector $\bm{y}_{k}:=[\bm{y}_{1,k},\bm{y}_{2,k},\dots,\bm{y}_{m,k}]^{\prime}$ and the $n\times m$ matrix $\bm{A}_{k}:=[\bm{a}_{1,k},\bm{a}_{2,k},\dots,\bm{a}_{m,k}]$, and letting $|\bm{z}|$ denote element-wise magnitude of a vector, the above measurement model can also be rewritten as

[TABLE]

where ′ denotes vector or matrix transpose.

The requirement that the measurement vectors used for imaging different $\bm{x}^{*}_{k}$’s be different and independent is what allows us to show that $m\ll n$ suffices. To understand this point in a simple fashion, consider the $r=1$ setting and suppose that $\bm{x}^{*}_{k}=\bm{x}^{*}_{1}$ (all columns are equal). We would then have $mq$ iid Gaussian measurements of $\bm{x}^{*}_{1}$ and hence $mq\geq Cn$ would suffice [9]. If $q=n$, this means just a constant number of measurements $m$ per column (signal) suffices. For $r>1$ but small, we will show that we can extend this idea to show that, when $q=n$ (or is larger), the required value of $m$ depends only on the value of $r$ and not on $n$. On the other hand, if $\bm{a}_{ik}=\bm{a}_{i}$ for all $k$, then, in the above $r=1$ example, only the first $m$ measurements are useful (the others are just repeats of these). Thus we will still need $m\geq Cn$ in this case. This $\bm{a}_{ik}=\bm{a}_{i}$ case, and its linear version, is what has been studied extensively in the literature [22, 23]. In this case, $m$ needs to be at least $(n+q)r$.

Let $\bm{X}^{*}\overset{\mathrm{SVD}}{=}{\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}$ denote its singular value decomposition (SVD) so that ${\bm{U}}^{*}\in\mathbb{R}^{n\times r}$, $\bm{B}^{*}\in\mathbb{R}^{r\times q}$, and $\bm{\Sigma}^{*}\in\mathbb{R}^{r\times r}$ is a diagonal matrix. Observe that this notation is a little non-standard, if the SVD was ${\bm{U}}^{*}\bm{\Sigma}^{*}{\bm{V}^{*}}^{\prime}$, we are letting $\bm{B}^{*}:={\bm{V}^{*}}^{\prime}$. Thus, columns of ${\bm{U}}^{*}$ and rows of $\bm{B}^{*}$ are orthonormal. We use ${\sigma_{\max}^{*}},{\sigma_{\min}^{*}}$ to denote the maximum, minimum singular values of ${\bm{X}^{*}}$ and $\kappa={\sigma_{\max}^{*}}/{\sigma_{\min}^{*}}$ to denote its condition number. Finally, we define

[TABLE]

We use the above non-standard notation for SVD because our solution approach will recover columns of $\tilde{\bm{B}}^{*}$, $\tilde{\bm{b}}^{*}_{k}$, individually by solving an $r$-dimensional standard PR problem (it is more intuitive to talk about recovery of column vectors than of rows). With the above notation, the QR decomposition of an estimate of $\tilde{\bm{B}}^{*}$, denoted $\hat{\bm{B}}$, will be written as $\hat{\bm{B}}\overset{\mathrm{QR}}{=}{\bm{R}}_{B}\bm{B}$ with $\bm{B}$ being an $r\times q$ matrix with orthonormal rows (or equivalently $\hat{\bm{B}}^{\prime}\overset{\mathrm{QR}}{=}\bm{B}^{\prime}({\bm{R}}_{B})^{\prime}$).

*Right Incoherence. * Observe that we have global measurements of each column, but not of the entire matrix. Thus, in order to correctly recover ${\bm{X}^{*}}$ with small $m$, we need an assumption that allows for correct “interpolation” across the rows. One way to ensure this is to borrow the right incoherence (incoherence or denseness of right singular vectors) assumption from matrix completion literature [24, 25]. In our notation, this means that we need to assume that

[TABLE]

with $\mu\geq 1$ being a constant. Clearly, this implies that

[TABLE]

for each $k$. If we assume $\kappa$ is a constant, up to constant factors, (4) also implies (3). Thus, up to constant factors, requiring right incoherence is the same as requiring that the maximum energy of any signal $\bm{x}^{*}_{k}$ is within constant factors of the average.

I-A1 Notation

We use $\|.\|$ to denote the $l_{2}$-norm of a vector or the induced 2-norm matrix and $\|.\|_{F}$ to denote the Frobenius norm. We use $\mathbbm{1}_{\mathrm{statement}}$ to denote the indicator function; it takes the value one if $\mathrm{statement}$ is true and is zero otherwise. A tall matrix with orthonormal columns is referred to as a “basis matrix”. For two basis matrices ${\bm{U}}_{1},{\bm{U}}_{2}$, we define the subspace error (distance) as

[TABLE]

This measures the sine of the largest principal angle between the two subspaces. We often use terms like “estimate ${\bm{U}}$” when the goal is to really estimate its column span, $\mathrm{Span}({\bm{U}})$. The phase-invariant distance is defined as

[TABLE]

For our guarantees, we work with real valued vectors and matrices, and in this case this simplifies to $\mathrm{dist}(\bm{x}^{*},\bm{\hat{x}})=\min(\|\bm{x}^{*}-\bm{\hat{x}}\|,\|\bm{x}^{*}+\bm{\hat{x}}\|)$. Define the corresponding distance between two matrices as

[TABLE]

We reuse the letters $c$, $C$ to denote different numerical constants in each use, with the convention $C\geq 1$ and $c<1$.

I-B Our Contributions and their Significance and Novelty

This work provides the first provably correct solution, AltMinLowRaP (Alternating Minimization for Low Rank PR), for Low Rank PR. AltMinLowRaP is a fast alternating minimization (AltMin) based solution approach with a carefully designed spectral initialization. We can prove that AltMinLowRaP converges geometrically to an $\epsilon$-accurate solution as long as (i) right incoherence stated in (3) holds, and (ii) the total number of available measurements, $mq$, is at least a constant (that depends on $\kappa,\mu$) times $nr^{4}\log(1/\epsilon)$. Its time complexity is order $mqnr\log^{2}(1/\epsilon)$, but if we replace $mq$ by its lowest allowed value, then the time complexity becomes $O(n^{2}r^{5}\log^{3}(1/\epsilon))$. If $q=n$ (or is larger), ignoring log factors, this implies that only about $r^{4}$ (or lesser) samples per signal suffice when using AltMinLowRaP. Moreover, when using these many samples, the time complexity per signal is only about $Cnr^{5}$. On the other hand, standard PR approaches (to recover each signal $\bm{x}^{*}_{k}$ individually) necessarily need $m\geq Cn$ samples, and order $Cn^{2}$ time, per signal [9, 10]. In the regime of small $r$, e.g., $r=\log n$, our result provides a significant sample, and time, complexity improvement over standard PR. Moreover, in this regime, our sample complexity is also only a little worse than its order optimal value of $(n+q)r$.We demonstate the practical power of AltMinLowRaP in Fig. 1, Fig. 2, and the other experiments described later.

The key insight that helps obtain the above reduction in sample complexity is the following observation: for both the initialization and the update steps for ${\bm{U}}^{*}$, conditioned on ${\bm{X}^{*}}$, we have access to $mq$ mutually independent measurements. These are not identically distributed (because the different $\bm{b}^{*}_{k}$’s could have different distributions), however, we can carefully use the right incoherence assumption to show that the distributions are “similar enough” so that concentration inequalities can be applied jointly for all the $mq$ samples. We also briefly study the dynamic setting, Phaseless Subspace Tracking, which allows the underlying signal subspace to change with time in a piecewise constant fashion.

To our best knowledge, even the linear version of our problem – recover ${\bm{X}^{*}}$ from $\bm{y}_{k}:=\bm{A}_{k}{}^{\prime}\bm{x}^{*}_{k},k=1,2,\dots,q$ – does not have any provably correct solutions (as we explain in Sec. I-D). Thus, our work also provides the first provable solution for this linear version. Our result provides an immediate corollary for this case as well. What has been studied extensively is the $\bm{A}_{k}=\bm{A}$ version of both LRPR and the above linear version [23, 22]. These are completely different problems as explained earlier in Sec. I-A.

I-B1 Significance

The other somewhat related problem to ours, Sparse PR, is actually quite different. This is because it involves recovery from global measurements of the sparse vector (each measurement depends on the entire unknown sparse vector), where as, in our case, the measurements are not global for the entire matrix ${\bm{X}^{*}}$. It is well known that, when studying iterative (non-convex) solutions to problems, the global measurements’ setting is easier to study, and one can obtain better sample complexity guarantees for it, as compared to its non-global counterpart [27]. For example, we can compare guarantees for iterative low rank matrix sensing (LRMS) from iid Gaussian linear projections with those for low rank matrix completion (LRMC) when assuming the iid Bernoulli model on observed entries [24, 25, 28]. LRMS can be solved using an iterative algorithm with nearly order optimal number of measurements, e.g., the approach of [28] needs $mq\geq C(n+q)r\log\max(n,q)$, while even the best iterative LRMC guarantee (under the iid Bernoulli measurement model) [26] needs $mq\geq C(n+q)r^{2}\log^{2}(nq)\log^{2}(1/\epsilon)$.

In this sense, the problem closest to ours that is extensively studied is LRMC. Of course LRMC involves recovery from completely local but linear measurements of ${\bm{X}^{*}}$, while LRPR involves recovery from nonlinear but column-wise global measurements. For this reason, for LRPR, in the regime of $q$ significantly larger than $n$, the required sample complexity $m$ is very small. As an example, suppose that $q\geq nr^{4}$, then we only need $m\geq C\max(r,\log q,\log(1/\epsilon))$. But this does not happen for LRMC.

We provide a comparison with the first and the best guarantees for non-convex (iterative) solutions for LRMC, sparse PR, and standard PR in Table I. These and other works are discussed in detail in Sec. I-D. As can be seen from the table, the first guarantee for iterative solutions to many problems is often sub-optimal (either needs more samples or more assumptions) compared to the best one that appeared later. Moreover, in the practical regime of $r$ being order $\log n$ or smaller, our LRPR sample complexity is as good or better than that of the best LRMC guarantee.

I-B2 Novelty

In the absence of relevant existing work for even solving the linear version of our problem, except for a convex solution for the $\bm{a}_{ik}=\bm{a}_{i}$ case (which is a significantly different problem), developing and analyzing our approach was not a straightforward extension of existing ideas. For example, the AltMinLowRaP algorithm itself is not just alternating standard PR over ${\bm{U}}$ and $\bm{B}$. The PR problem for recovering ${\bm{U}}^{*}$ given an estimate of $\tilde{\bm{B}}^{*}$ is significantly different from standard PR; see Sec. II-A1 below.

For the above reasons, it is also not possible to directly modify proof techniques from existing work. We borrow some ideas from LRMC [25] and standard PR results [10, 9]. But the major difference is that concentration bounds need to applied differently than for either of these problems. (i) The LRMC guarantees use results for Bernoulli random variables (which is a much more well-developed literature that has also been studied in the context of random graphs). In our setting, the random variables are not Bernoulli and not even bounded. Hence we rely on the sub-exponential Bernstein inequality [29, Theorem 2.8.1] and the fact that the product of two sub-Gaussian random variables is a sub-exponential [29, Lemma 2.7.7]. A second difference is that LRMC results do not need to deal with the phase error term. (ii) Standard PR results do have a phase error term and do deal with unbounded random variables using results from [29]. But they do not have to prove concentration using a set of $mq$ measurements that are not identically distributed and, on first glance, may not even be “similar enough” to get a useful result. The “similarity” that is needed is of the following form: the maximum sub-exponential norm of any of the $mq$ random variables being summed is not much larger than its average value. For each term, we have to carefully exploit the right incoherence assumption to show that this holds.

I-C Motivation for studying Low Rank PR (LRPR)

Low rank is a commonly used model in many dynamic biomedical imaging applications since (i) such images cannot change too much from one frame to the next, and (ii) these images are taken in controlled settings and so there are no fast changing foreground occlusions to worry about111Occlusions by moving objects or persons in the foreground are a common feature in computer vision problems such as surveillance or autonomous vehicles etc; for such videos a sparse + low-rank model is more appropriate. For example, it is an important part of many practically useful fast compressive dynamic MRI solutions, e.g., see [19], and follow-up works [30, 31]222These follow-up works exploit both low-rank of the entire sequence and wavelet sparsity of each image to further reduce the number of measurements needed in practice. This is the so-called “sparse and low-rank” model which is very different from sparse+low-rank model where the sparse component models occlusions by foreground moving objects.. In a similar fashion, a low sample complexity solution to LRPR can enable fast or low-cost dynamic phaseless imaging in applications such as solar imaging when the sun’s surface properties gradually change over time [32], or Fourier ptychographic imaging of live biological specimens and other dynamic scenes [33, 21]. Suppose the scene resolution is $n$ and the total number of captured frames is $q$. If the dynamics is approximated to be linear and slow changing, with most of the change being explained by $r$ linearly independent factors, then the matrix formed by stacking the vectorized image frames next to each other can be modeled as a rank-$r$ matrix plus small modeling error. In typical settings, $r\ll\min(n,q)$ is a valid assumption, making the unknown images’ matrix approximately low-rank.

In all the above applications, measurement acquisition is either expensive or slow. For example, Fourier ptychography is a technique for super-resolution in which each of a set of low resolution cameras measures the magnitude of a different band-pass filtered version of the target high-resolution image. To get enough measurements per image, one either needs many cameras (expensive), or one needs to move a single camera to different locations to acquire the different bands [33]. This can make the acquisition process slow. By exploiting the low-rank assumption, it is possible to get an accurate reconstruction with using fewer total measurements (fewer cameras in this example). This has been demonstrated experimentally for dynamic Fourier ptychography in our recent work [21] and its follow-up [34]. Moreover, it is indeed practically valid to assume that a different measurement matrix $\bm{A}_{k}$ is used for each different signal/image. In the ptychography example, this would correspond to using a different randomly selected subset of cameras at different times $k$. Modified cameras can also be designed that save power by switching off a different set of pixels at different time instants. We have explored both settings in [21].

Another practical point that should be mentioned is that, often, in practice, a very small value of rank $r$ suffices. For example, we used $r=20$ in all our experiments on image sequences with n = 32400 in [21]. In follow-up work [34], we show that just $r=5$ suffices for the same datasets, as long as a “modeling error correction step”333This step applies a few iterations of any standard PR approach column-wise to the output of AltMinLowRaP, in order to recover some of the “modeling error” in the low-rank assumption. is applied to the output of AltMinLowRaP.

Lastly, in comparison to sparsity or structured sparsity priors, the low rank prior is a significantly more flexible one since it does not require knowledge of the dictionary or basis in which the signal is sufficiently sparse. In Table III, we demonstrate this via a simple experiment. We compare AltMinLowRaP with the most recent provable sparse PR algorithm [18], CoPRAM, applied with assuming wavelet sparsity (which is a generic choice for any piecewise smooth image, but is not necessarily the best choice for the particular image). As can be seen, AltMinLowRaP has significantly superior performance not just for the real image sequence, but also for its deliberately sparsified version. The sparsified sequence had sparsity level $s\approx 0.1n\approx 100$ and we provided CoPRAM with this ground truth. AltMinLowRaP used just $r=15$ for all three results and still had much lower reconstruction error than CoPRAM. Details of this experiment are provided in Sec. IV-B. Moreover, low-rank also includes certain types of dynamic sparsity models (those with fixed of very small changes in support over time) as special cases.

I-D Review of Related Work

I-D1 Linear version of our problem: Compressive PCA

While one would think that the linear (with phase) version of our problem would been extensively studied, this is not true. There have been a few algorithmic solutions for this problem in prior work [35, 36], and attempts to prove some facts theoretically. Follow-up work consists of an Asilomar 2014 paper [37] that solves the general PCA problem for any (not necessarily low rank) matrix ${\bm{X}^{*}}$, but does not discuss recovery of ${\bm{X}^{*}}$. We explain these in detail in Sec. II-C.

I-D2 Our measurement model, but with same set of $m$ measurement vectors used for all signals, and its linear version

The covariance sketching problem, e.g., see [23], assumes that measurements satisfying (1), but with $\bm{a}_{ik}=\bm{a}_{i}$, are available. One aggregates these over $k$ to get $\bm{y}_{i}:=\sum_{k}\bm{y}_{ik}=\bm{a}_{i}^{\prime}(\sum_{k}\bm{x}^{*}_{k}\bm{x}^{*}_{k}{}^{\prime})\bm{a}_{i}=\bm{a}_{i}^{\prime}{\bm{X}^{*}}{\bm{X}^{*}}^{\prime}\bm{a}_{i}$. This aggregation is what ensures that the memory complexity of storing the measurements is order $m$ and not $mq$ (which is what we need). Also, the aggregated $\bm{y}_{i}$ is a function of ${\bm{X}^{*}}{\bm{X}^{*}}^{\prime}/q$ only in the $\bm{a}_{ik}=\bm{a}_{i}$ setting, otherwise it is a meaningless quantity. Assuming random zero mean iid signals $\bm{x}^{*}_{k}$, ${\bm{X}^{*}}{\bm{X}^{*}}^{\prime}/q$ is the empirical covariance matrix of a signal. The question is can we recover ${\bm{X}^{*}}{\bm{X}^{*}}^{\prime}/q$ from the scalar sketches $\bm{y}_{i},i=1,2,\dots m$ with using $m$ much smaller than $nq$, when ${\bm{X}^{*}}$ is low rank (or has other structure)? When ${\bm{X}^{*}}$ is rank $r$, the result of [23] proves that $m$ of order $(n+q)r$ suffices to estimate the empirical covariance from the aggregated measurements $\bm{y}_{i}$ if one solves an appropriately defined nuclear norm minimization problem. For solving LRPR, we need a much smaller $m$ than this. The reason is we assume independent $\bm{a}_{ik}$’s for different $k$, and we assume we have access to each individual $\bm{y}_{ik}$.

The linear version of the above problem, but with random noise added, is considered in Corollary 3 of [22] and in the remark immediately below it. In our notation, its measurement model can be written as $\bm{y}_{ik}=\langle\bm{a}_{i}\bm{e}_{k}{}^{\prime},{\bm{X}^{*}}\rangle+\bm{w}_{ik}$ where $\bm{w}_{ik}$ is iid zero mean Gaussian noise with variance $\nu^{2}$. This result (specialized to the exact low rank case) shows that, whp, a nuclear norm minimization based solution will recover an estimate $\hat{\bm{X}}$ of ${\bm{X}^{*}}$ that satisfies $\|{\bm{X}^{*}}-\hat{\bm{X}}\|_{F}^{2}\leq C\ \nu^{2}r(n+q)/m.$ In this paper, the focus is on using the low rank property to achieve noise robustness. If the low rank property was not used, and one attempted to recover the columns individually, the recovery error bound would scale as $\nu^{2}nq/m$ which is much larger. This paper also studies other settings of recovering an approximately low rank matrix from linear measurements.

I-D3 Tangentially related work

Some other tangentially related work includes: (i) computing the approximate rank $r$ approximation of any matrix (need not be low rank) from its random sketches [38, 39] (sketched SVD); (ii) compressed covariance estimation using different sketching matrices for each data vector, but without the low-rank assumption [40]; and (iii) a generalization of low-rank covariance sketching [41]: this attempts to recover an $n\times r$ matrix ${\bm{U}}^{*}$ from measurements $\bm{y}_{i}=\|\bm{a}_{i}^{\prime}{\bm{U}}^{*}\|^{2}$ with $r\ll n$. When $r=1$, this is the standard PR problem. In the general case, this is related to covariance sketching described above, but not to our problem.

I-D4 Linear low-rank matrix recovery – LRMS and LRMC

Low-rank matrix recovery problems with linear measurements that have been extensively studied can be split into two kinds - those with “global measurements” and those without. “Global measurements” means that each measurement contains information about the entire structured quantity-of-interest, here the low-rank matrix. Such problems are called “affine rank minimization problems” or “low-rank matrix sensing” (LRMS) and involve recovery of ${\bm{X}^{*}}$ from $\bm{y}_{i}=\langle\bm{A}_{i},{\bm{X}^{*}}\rangle$ with $\bm{A}_{i}$ being dense matrices (typically iid Gaussian), see for example, [42, 28, 43, 44, 25]. More recent work studies the case of $\bm{A}_{i}=\bm{a}_{i}\bm{a}_{i}{}^{\prime}$ [45, 46]. Low-rank Matrix Completion (LRMC) is the completely local measurements’ setting that involves recovering ${\bm{X}^{*}}$ from measurements of a randomly (iid Bernoulli) selected subset of its entries [24, 47, 25, 48, 26, 12] . Thus $\bm{A}_{i}$’s are one-sparse matrices.

A precursor to LRMS is compressive sensing (CS) of sparse signals. This has the same property as LRMS, it involves recovering a sparse $\bm{x}^{*}$ from $\bm{y}_{i}=\langle\bm{a}_{i},\bm{x}^{*}\rangle$ with $\bm{a}_{i}$ being dense (sub-)Gaussian random vectors. Similar to CS, even for LRMS, it is typically possible to prove a simple (sparse or low-rank) restricted isometry property which simplifies the rest of the analysis. Our problem setting is different from, and more difficult than, LRMS. There are no “global measurements” of the entire ${\bm{X}^{*}}$ and, moreover, the measurements’ phase/sign is unknown. In this sense it is closer to LRMC than to LRMS. But, unlike LRMC, we do have column-wise global measurements. This is why, for our problem only incoherence of right singular vectors suffices, while LRMC needs incoherence of both left and right singular vectors.

In summary, the problem closest to ours that is well studied is LRMC. The first iterative solution to LRMC was [47]. However, its guarantee does not bound the required number of algorithm iterations, and thus its time complexity cannot be bounded. The first iterative LRMC solution with bounded time complexity, AltMinComplete [25], needs a sample complexity of about $C\kappa^{4}\mu^{2}nr^{4.5}\log(1/\epsilon)$ and assumes sample-splitting (a different independent set of measurements is used at each iteration). The most recent work on LRMC [12] removes the sample-splitting requirement and has bounded time complexity, but its sample complexity was $Cnr^{3}\log^{6}(n)$. The best iterative LRMC solution in terms of sample complexity [26] needs $Cnr^{2}\log^{2}n\log^{2}\kappa/\epsilon$ samples but needs sample-splitting. In the practical regime $r\in O(\log n)$, clearly our sample complexity for LRPR is comparable to the best LRMC guarantee [26]. For all values of $r$, it is slightly better than the first bounded time iterative LRMC solution [25]. Time-wise the AltMin algorithm for LRMC is faster than our AltMinLowRaP algorithm for LRPR (see Table I). This is because the LS problem to be solved in each AltMin step of LRMC involves a matrix with a large number of zeros but this is not the case for LRPR.

I-D5 Sparse PR

Sparse PR is a somewhat related problem to ours since it involves PR with a different type of structural assumption on the signals. But, as noted earlier, there is a major difference. Sparse PR involves recovery from global measurements of the entire sparse vector. It can be understood as the phaseless version of Compressive Sensing with random Gaussian measurements. The global measurements’ setting is typical easier than its non-global counterparts. Provably correct sparse PR approaches include convex relaxation approaches such as $\ell_{1}$-PhaseLift [14]; older combinatorial methods [49]; and a series of fast iterative approaches: (i) AltMinSparse [6], (ii) Sparse Truncated Amplitude Flow (SPARTA) [16], (iii) Thresholded WF [17] and CoPRAM [18]. The first two fast nonconvex approaches – AltMinSparse and SPARTA – needed to assume a lower bound on the minimum nonzero entry of $\bm{x}$. In follow-up work on Thresholded WF and then on CoPRAM, this extra assumption was removed. All four results need at least order $s^{2}\log n$ measurements. A summary of comparison of our work with LRMC, sparse PR, and standard PR is provided in Table I.

I-E Organization

We present our algorithm and guarantee along with a detailed discussion of the novelty of our proof techniques in Sec. II given next. The overall proof is given in Sec. III. The lemmas introduced in Sec. III are proved in Appendix A. Numerical experiments are provided in Sec. IV. We develop extensions to phaseless subspace tracking in Sec. V. We conclude in Sec. VI with a detailed discussion of ongoing and future work.

II Low Rank PR: Algorithm and Guarantee

II-A AltMinLowRaP algorithm

The complete algorithm is summarized in Algorithm 1. We explain its main idea next and then explain the details.

II-A1 Main idea

AltMinLowRaP minimizes the following

[TABLE]

alternatively over ${\bm{U}},\bm{B}$ with the constraint that ${\bm{U}}$ is a basis matrix. To initialize, we develop a spectral initialization for $\mathrm{Span}({\bm{U}}^{*})$ explained below. At a top level, the alternating minimization (AltMin) can be understood as alternating PR: minimize (5) over $\bm{B}$ keeping ${\bm{U}}$ fixed at its current value and then vice versa. But there are important differences between the two PR problems and how they can be solved.

• •

Given an estimate of $\mathrm{Span}({\bm{U}}^{*})$, denoted ${\bm{U}}$, and assuming that ${\bm{U}}$ contains orthnormal columns and is independent of the measurement vectors, the recovery of each $\tilde{\bm{b}}^{*}_{k}$ is an $r$-dimensional “standard PR” problem. We can use either of [9, 10] to solve it. The estimate that we get, denoted $\bm{\hat{b}}_{k}$, is actually an estimate of $\bm{g}_{k}:={\bm{U}}{}^{\prime}{\bm{U}}^{*}\tilde{\bm{b}}^{*}_{k}$ which is a rotated version of $\tilde{\bm{b}}^{*}_{k}$. If $\sin\Theta({\bm{U}},{\bm{U}}^{*})\leq\delta$, then, by the noisy PR result of [10], we can see that $\|\bm{\hat{b}}_{k}-\bm{g}_{k}\|\leq C\delta\|\bm{b}^{*}_{k}\|$.

• •

Given a previous estimate of $\tilde{\bm{B}}^{*}$, the update of ${\bm{U}}^{*}$, or equivalently of its vectorized version, ${\bm{U}}^{*}_{vec}$, is a significantly non-standard PR problem for two reasons. First, the “measurement vectors” for this PR problem are no longer independent or identically distributed. Second, and more importantly, by using the previous estimates of ${\bm{U}}^{*}$ and of $\tilde{\bm{b}}^{*}_{k}$, with accuracy level $\delta$, we can get an estimate, $\bm{\hat{x}}_{k}={\bm{U}}\bm{\hat{b}}_{k}$, of $\bm{x}^{*}_{k}$ with the same accuracy level. With this, we can also get an estimate of the phase/sign of the measurements, $\bm{c}_{ik}:=\mathrm{phase}(\bm{a}_{ik}{}^{\prime}\bm{x}^{*}_{k})$ with the same accuracy level. As a result, obtaining a new estimate of ${\bm{U}}^{*}_{vec}$ becomes a much simpler Least Squares (LS) problem rather than a PR problem.

• •

We stress here that an argument similar to the above does not apply when recovering $\tilde{\bm{b}}^{*}_{k}$’s. The reason is, with a new estimate of ${\bm{U}}^{*}$, denoted ${\bm{U}}^{+}$, the previous estimate of $\tilde{\bm{b}}^{*}_{k}$ becomes useless: (i) it is close to ${\bm{U}}^{\prime}{\bm{U}}^{*}\tilde{\bm{b}}^{*}_{k}$, and not to ${\bm{U}}^{+}{}^{\prime}{\bm{U}}^{*}\tilde{\bm{b}}^{*}_{k}$; and (ii) we only estimate the span of the columns of ${\bm{U}}^{*}$ accurately, so ${\bm{U}}$ or ${\bm{U}}^{+}$ are close to ${\bm{U}}^{*}$ (and hence to each other) only in the subspace error $\sin\Theta$ (and not in spectral or Frobenius norm). As a result an estimate of the form ${\bm{U}}^{+}\bm{\hat{b}}_{k}$ cannot be shown to be close to $\bm{x}^{*}_{k}$.444To understand this point easily, suppose both ${\bm{U}}$ and ${\bm{U}}^{+}$ are perfect estimates of ${\bm{U}}^{*}$ in terms of the subspaces they span. Suppose ${\bm{U}}^{+}={\bm{U}}^{*}{\bm{R}}_{2}$ and ${\bm{U}}={\bm{U}}^{*}{\bm{R}}_{1}$ where ${\bm{R}}_{1},{\bm{R}}_{2}$ could be any rotation matrices. Then it is easy to see that $\|\bm{x}^{*}_{k}-{\bm{U}}^{+}\bm{\hat{b}}_{k}\|\geq\|({\bm{U}}-{\bm{U}}^{+})\bm{\hat{b}}_{k}\|-\|\bm{x}^{*}_{k}-{\bm{U}}\bm{\hat{b}}_{k}\|\geq\|{\bm{U}}^{*}({\bm{R}}_{1}-{\bm{R}}_{2})\bm{\hat{b}}_{k}\|-\delta\|\bm{x}^{*}_{k}\|$. Since ${\bm{R}}_{1},{\bm{R}}_{2}$ can be any rotation matrices, e.g, one could have ${\bm{R}}_{2}=-{\bm{R}}_{1}$; in this case, the above error is lower bounded by $2\|\bm{\hat{b}}_{k}\|-\delta\|\bm{x}^{*}_{k}\|\geq(2-C\delta)\|\bm{x}^{*}_{k}\|$. The last inequality follows since $\|\bm{\hat{b}}_{k}\|=\|\bm{\hat{x}}_{k}\|$ and $\mathrm{dist}(\bm{x}^{*}_{k},\bm{\hat{x}}_{k})\leq C\delta\|\bm{x}^{*}_{k}\|$ by Lemma 3.2. The error can thus be even higher than using a zero vector to estimate $\bm{x}^{*}_{k}$. .

II-A2 Details

A different way to understand AltMinLowRaP is to split it into a three-way AltMin problem over ${\bm{U}}^{*}$, $\tilde{\bm{b}}^{*}_{k}$’s, and $\bm{c}_{ik}$’s. This discussion assumes “sample-splitting”: a new set of $mq$ measurements is used for each update of ${\bm{U}}^{*}$ and another new set for each update of $\tilde{\bm{B}}^{*}$. Thus the total number of measurements used is $2mq$ times the number of iterations.

At each new iteration, we first obtain a new estimate of ${\bm{U}}^{*}$ using previous estimates of $\tilde{\bm{b}}^{*}_{k}$’s and of the measurements’ phases $\bm{c}_{ik}$’s. This is a LS problem, see line 10 of Algorithm 1. The output of the LS step may not have orthonormal columns; this is easily resolved by a QR decomposition step after it (line 11). 2. 2.

Given a new estimate of ${\bm{U}}^{*}$, we recover each $\tilde{\bm{b}}^{*}_{k}$, by solving easy individual $r$-dimensional standard PR problems. These are easy because $m$ of order $r$ suffices. 3. 3.

Given a good estimate of $\tilde{\bm{b}}^{*}_{k}$ and of ${\bm{U}}^{*}$, we can get an equally good estimate of $\bm{x}^{*}_{k}$ and hence of the signs/phases $\bm{c}_{ik}$’s.

Consider the PR step to update $\tilde{\bm{b}}^{*}_{k}$’s. Observe that we can rewrite $\bm{y}_{ik}$ as $\bm{y}_{ik}=|\langle\bm{a}_{ik}{},{\bm{U}}^{*}\tilde{\bm{b}}^{*}_{k}\rangle|=|\langle({{\bm{U}}^{*}}^{\prime}\bm{a}_{ik}),\tilde{\bm{b}}^{*}_{k}\rangle|$. If ${\bm{U}}^{*}$ were known, we would have a noise-free standard PR problem. If, at the $t$-the iteration, instead, we have a good estimate of $\mathrm{Span}({\bm{U}}^{*})$, denoted ${\bm{U}}^{t}$, we can still recover the $\tilde{\bm{b}}^{*}_{k}$’s by solving a noisy version of the same problem. Due to sample splitting, ${\bm{U}}^{t}$ is independent of $\bm{a}_{ik}$’s and so the design vectors $({{\bm{U}}^{t}}^{\prime}\bm{a}_{ik})$ are still iid standard Gaussian. Any standard PR solution can be used, here we use Reshaped Wirtinger flow (RWF) [10]. The noise seen by RWF is proportional to $\sin\Theta({\bm{U}}^{t},{\bm{U}}^{*})$. The error in the output of RWF cannot be lower than this value [10]. Thus, one needs to use just enough iterations of RWF so that the error at the end of the final RWF iteration is proportional to $\sin\Theta({\bm{U}}^{t},{\bm{U}}^{*})$. Since we prove geometric convergence of $\sin\Theta({\bm{U}}^{t},{\bm{U}}^{*})$, we can let $T_{RWF,t}$ grow linearly with $t$.

II-A3 Initialization

To obtain the initialization, we develop a careful modification of the truncated spectral initialization idea from [9, 20]. First assume that $r$ is known. We initialize $\hat{\bm{U}}$ as the top $r$ left singular vectors of the following matrix:

[TABLE]

where $C_{Y}$ is a constant that decides the truncation threshold (which measurements are too large in magnitude compared to the mean energy of the measurements and should be discarded). For our guarantee, we set it equal to $9\kappa^{2}\mu^{2}$. In practice a good value can be chosen by experimentation and cross-validation or by using the ideas in [50]. An alternative approach is to set it differently for each $k$ as done in [20]. This approach does not require knowledge of $\kappa$ or $\mu$, but it results in a worse lower bound on just $m$. We discuss the effect of this choice in Remark 2.2 and also in Sec. II-D3.

To understand why the above approach works, first consider the above matrix with the indicator function removed. Then it is not hard to see that its expected value equals $(1/q)[{\bm{U}}^{*}({\bm{\Sigma}^{*}}^{2}){{\bm{U}}^{*}}^{\prime}+2\mathrm{trace}({\bm{\Sigma}^{*}}^{2})\bm{I}]$, and so its span of top $r$ singular vectors equals $\mathrm{Span}({\bm{U}}^{*})$. Hence, with large enough $mq$, the same should approximately hold for the original matrix. However, when using $\bm{Y}_{U}$ with the indicator function removed, a few “bad” measurements (those with very large magnitude $\bm{y}_{ik}^{2}$ compared to their empirical mean over $i,k$) can heavily bias its value. To mitigate this effect, and get a good initialization in spite of it, we will need a larger value of $mq$. Using the indicator function helps truncate the summation to only sum over the “good” measurements, and as a result a smaller value of $mq$ suffices. Mathematically, this helps ensure that $\bm{Y}_{U}$ is close to a matrix that can be written as $\sum_{ik}\bm{w}_{ik}\bm{w}_{ik}{}^{\prime}$ with $\bm{w}_{ik}$’s being iid sub-Gaussian vectors (instead of sub-exponential in the case without truncation) [9].

We can also use $\bm{Y}_{U}$ to correctly estimate $r$ whp by using the fact that, when $m$ and $q$ are large, the gap between its $r$-th and $(r+1)$-th singular value is close to ${\sigma_{\min}^{*}}^{2}/q$. With this idea, we estimate $r$ as given in the first step of Algorithm 1. As explained in [20], another way of estimating the rank is to set $\hat{r}=\arg\max_{j}(\sigma_{j}(\bm{Y}_{U})-\sigma_{j+1}(\bm{Y}_{U}))$. This approach does not require knowledge of any model parameters. Hence it is easily applicable for real data (even without training samples being available). However, it works under the assumption that consecutive nonzero singular values of $\bm{X}^{*}$ are close (do not have significant gap), see [20, Corollary 3.7] for one precise statement of this claim.

The main idea of the initialization step explained above was first developed in our previous work [20], we explain the difference later in Sec. II-D3. Also, as pointed out by an anonymous reviewer, a matrix that is related to $\bm{Y}_{U}$ was used in earlier work [35, 36] to try to solve what can be called the linear version of our problem; see Sec. II-C. We were not aware of this work when developing our approach. Our approach was developed independently in [20] as a modification of the truncated spectral initialization idea from PR literature [9] and then modified here.

We summarize the complete algorithm in Algorithm 1. As is commonly done in existing literature, e.g., see [25, 6], in order to obtain a provable guarantee in a simple fashion, we assume sample-splitting. Since we prove geometric convergence of the iterates, this increases the required sample complexity by a factor of only $\log(1/\epsilon)$. In our empirical evaluations, we reuse the same set of measurements.

II-B Guarantee

We have the following guarantee.

Theorem 2.1**.**

Consider Algorithm 1. Assume that the $\bm{y}_{ik}$’s satisfy (1) with $\bm{a}_{ik}$ being iid standard Gaussian; and ${\bm{X}^{*}}$ is an $n\times q$ rank-$r$ matrix that satisfies right-incoherence with parameter $\mu$. Set $T:=C\log(1/\epsilon)$, $T_{RWF,t}=C(\log r+\log\kappa+t(\log(0.7)/\log(1-c)))$, $\omega=1.3{\sigma_{\min}^{*}}^{2}/q$, and $C_{Y}=9\kappa^{2}\mu^{2}$ in (6). Assume that, for the initialization step and for each new update, we use a new set of $m$ measurements with $m$ satisfying $mq\geq C\kappa^{12}\mu^{4}\cdot nr^{4}$ and $m\geq C\max(r,\log q,\log n)$. Then, with probability (w.p.) at least $1-Cn^{-10}$,

[TABLE]

and $\mathrm{dist}(\bm{\hat{x}}_{k}^{T},\bm{x}_{k}^{*})\leq\epsilon\|\bm{x}_{k}^{*}\|$ for each $k$. Moreover, after the $t$-th iteration,

[TABLE]

where $\delta_{\mathrm{init}}=\frac{c}{\kappa^{2}r}$. Similar bounds also hold on the error in estimating $\bm{x}^{*}_{k}$s. The time complexity is $mqnr\log^{2}(1/\epsilon)$.

Proof: We prove this in Sec. III.

Remark 2.2**.**

If we are willing to tolerate another lower bound on just $m$ of $m\geq Cr^{4}$, we can improve the dependence on $\kappa$, $\mu$ to $\kappa^{8}\mu^{2}$. This will require the following change: define $\bm{Y}_{U}$ as done in [20]: use $9\sum_{i=1}^{m}\bm{y}_{ik}^{2}/m$ as the threshold inside the indicator function. A second advantage of using this is that it allows one to use $C_{Y}=9$ instead of $C_{Y}$ that depends on $\kappa,\mu$. The disadvantage is of course a more stringent lower bound on $m$.

Remark 2.3**.**

To understand the time complexity, observe that the most expensive step at each iteration is the update of ${\bm{U}}^{*}$. This requires solving an LS problem of recovering an $nr$ length vector from $mq$ measurements. One can solve this by conjugate gradient descent with a cost of $mq\cdot nr\cdot\log(1/\epsilon)$ [6]. This times the total number of iterations, $T=C\log(1/\epsilon)$ gives the complexity of the algorithm after initialization. Consider the initialization step. Observe that the top $r$ singular vectors of $\bm{Y}_{U}$ are also the top $r$ singular vectors of an $mq\times n$ matrix $\bm{G}$ whose columns are given by $\bm{a}_{ik}\bm{y}_{ik}$ times the indicator function used in (6). Thus $\bm{Y}_{U}=\bm{G}\bm{G}^{\prime}$. Since computing the $r$-SVD of an $a\times b$ matrix to $\delta$ accuracy needs time of order $abr\log(1/\delta)$ [51], the SVD needed for initialization can be computed in time $Cmq\cdot n\cdot r\cdot\log(1/\delta_{\mathrm{init}})=Cmqnr\log r$ where $\delta_{\mathrm{init}}$ is the error level to which the initialization needs to be accurate. As explained later, $\delta_{\mathrm{init}}=c/r$ suffices.

Theorem 2.1 implies that one can achieve geometric convergence as long as the sample complexity $m_{\mathrm{tot}}:=(2T+1)m$ satisfies $m_{\mathrm{tot}}q\geq C\kappa^{12}\mu^{4}nr^{4}\log(1/\epsilon)$ along with $m_{\mathrm{tot}}\geq C\max(r,\log q,\log n)\log(1/\epsilon)$. The second lower bound is very small and essentially redundant555We need this lower bound because we recover the $q$ $\tilde{\bm{b}}^{*}_{k}$’s individually by solving a standard PR problem for each. This step works correctly w.p. at least $1-2q\exp(r-cm)$. except when $q\geq Cnr^{4}$.

Notice that the LRPR sample complexity is significantly better than that of standard (unstructured) PR methods which necessarily need $m=Cn$ samples per signal (matrix column). For fixed $m$ and $q$, LRPR time complexity is about $r$ times worse than that of standard PR. But, if we use the smallest value of $mq$ needed by each method to get an $\epsilon$-accurate estimate, AltMinLowRaP is actually faster when $r$ is small: its needs time of order $n^{2}r^{5}\log(1/\epsilon)$ while standard PR methods need time of oder $n^{2}q$. We demonstrate this fact experimentally in Fig 2.

The minimum number of samples needed to recover an $n\times q$ matrix of rank $r$ is $(q+n)r$. Thus, in general, our sample complexity is $r^{3}\log(1/\epsilon)$ times worse than its order-optimal value. As noted earlier, in problem settings like ours, where the measurements are not global, non-convex algorithms typically do need more than the order-optimal number of samples. Since neither our problem nor its linear version have any complete provable guarantees for correct recovery in existing work, LRMC is the closest problem to ours that has been extensively studied and that also uses non-global measurements. Sparse PR is the other somewhat related problem to ours but it is easier because it involves recovery from global measurements of a sparse vector. Table I provides a comparison of our guarantee with the first and best results for both problems. More details are in Sec. I-D. As can be seen from the table, our sample complexity compares favorably with that for the first non-convex solution for LRMC that has bounded time complexity (bounds the required number of iterations) [25]. In the practical regime of $r$ being order $\log n$, it even compares with the best iterative LRMC result [26]. Also, the the first guarantee for iterative solutions to both problems is sub-optimal (needs more samples or more assumptions) compared to the best one that appeared later.

We should reiterate that, (i) in many practical applications, a small value of $r$ suffices, e.g., we used $r=5$ for images with $n=32400$ in [34] followed by a few iterations of model error correction via standard PR; and (ii) low rank is a more flexible model for dynamic imaging than sparsity. Also see Tables II and III.

II-C Linear version of our problem: Compressive PCA

Consider the linear (with phase) version of our problem – recover a low rank matrix ${\bm{X}^{*}}$, or its column span, from

[TABLE]

This problem can also be understood as “compressive PCA” although compressive PCA typically allows ${\bm{X}^{*}}$ to also be only approximately low rank [35, 36, 37]. Clearly, both our algorithm and our guarantee, Theorem 2.1, directly apply to this simpler special case as long as the $\bm{a}_{ik}$’s are iid Gaussian. Even when phases are available one could use the magnitude measurements and solve our harder problem instead. We can state the following corollary.

Corollary 2.4**.**

Consider the problem of recovering an $n\times q$ rank-$r$ matrix ${\bm{X}^{*}}$ from $\bm{y}_{ik}:=\langle\bm{a}_{ik},\bm{x}^{*}_{k}\rangle,i=1,2,\dots,m,\ k=1,2,\dots,q$. As long as right incoherence given in (3) holds, one needs $m\geq C\frac{n}{q}r^{4}\log(1/\epsilon)$ to recover ${\bm{X}^{*}}$ and its column span to $\epsilon$ accuracy (precisely defined in Theorem 2.1).

Thus our work also provides a simple and fast AltMin solution that provably converges geometrically to the compressive PCA solution. Of course, when the phase/sign is known, a simpler version of our algorithm should work and a better guarantee should be obtainable. We discuss this further in the conclusions’ section.

In terms of existing related work, [35] provided an approach for solving the above problem when ${\bm{X}^{*}}$ is only approximately low-rank. A follow-up paper [36] studied its modification where $\bm{a}_{ik}$’s are sparse random vectors, e.g., “sparse Bernoulli” as they call it (each entry takes values $-1,0,1$ with probabilities $1/2s,1-1/s,1/2s$). Their proposed approach is related to our initialization step: one computes an estimate of the principal subspace of the column span of ${\bm{X}^{*}}$ by computing the top $r$ singular vectors of a matrix that is similar to $\bm{Y}_{U}$ defined in (6). The difference is that the indicator function is absent, and they also sum over terms of the form $\bm{a}_{ik}\bm{y}_{ik}\bm{y}_{jk}\bm{a}_{jk}{}^{\prime}$ for all $j\neq i$ which, in this linear setting, have nonzero expected value. To be precise, they compute the top $r$ singular vectors of $\sum_{k}\bm{A}_{k}\bm{y}_{k}\bm{A}_{k}{}^{\prime}\bm{y}_{k}{}^{\prime}$. Denote such a matrix by $\bm{Y}_{U,modified}$. These works treat the columns $\bm{x}^{*}_{k}$ as random vectors and assume that they have nonzero mean and provide a simple intuitive approach to estimate the mean. This is of course only possible because the measurements are linear, and cannot be done in our setting. Both works [35, 36] show that the expected value of $\bm{Y}_{U,modified}$ is equal to $c_{1}\bm{I}+c_{2}\mathbb{E}[{\bm{X}^{*}}{\bm{X}^{*}}^{\prime}/q]$, thus, the span of top r eigenvectors of the expected value equals the desired subspace: principal subspace of the population covariance, $\mathbb{E}[{\bm{X}^{*}}{\bm{X}^{*}}^{\prime}/q]$. They also claim that as $q\rightarrow\infty$, $\bm{Y}_{U,modified}$ converges to its expected value, but do not provide a rate of convergence. Finally, they also upper bound the expected value of the Frobenius norm of the error between $\bm{Y}_{U,modified}$ and its expected value for any value of $q$. Since this result also only bounds the expected value of the error, it cannot be used to get any useful information about the sample complexity $m$ that is required even for the initialization step. Moreover, of course these works do not provide any guarantee on how to recover the entire matrix ${\bm{X}^{*}}$.

A later work in Asilomar 2014 [37] attempted to solve a generalization of the above problem: it provided an approach (which is again somewhat related to only our initialization step) and a guarantee for recovering the top $r$ singular vectors of a general matrix ${\bm{X}^{*}}$. Specialized to the exactly rank $r$ case, their result proves that, if each column of ${\bm{X}^{*}}$ is bounded, and if $m\geq\frac{1}{\epsilon}\max\left(\frac{n}{\sqrt{q}}\sqrt{r}\kappa,\frac{n}{q}r^{2}\kappa^{2}\right)$, one can obtain an $\epsilon$-accurate recovery of the subspace $\mathrm{Span}({\bm{U}}^{*})$. This is a much weaker result than ours: (i) its sample complexity $m$ depends on $1/\epsilon$ instead of on $\log(1/\epsilon)$; and (ii) $m$ needs to grow as $n/\sqrt{q}$ instead of as $n/q$. The reason it is weaker of course, is because they used a single step approach instead of an iterative algorithm.

II-D Discussion of Proof Techniques and Reason for Worse Dependence on $\kappa,r$

II-D1 Proof techniques

Since even the linear version of our problem has not been studied theoretically (except for a convex solution for the $\bm{a}_{ik}=\bm{a}_{i}$ case which is a significantly different problem), it is not possible to directly modify proof techniques from existing work. We do borrow some ideas from LRMC [25] or from standard PR results [10, 9]. But, as explained earlier in Sec. I-B2, concentration bounds need to applied in a significantly different way than for either of these problems. The algebra for obtaining an expression that bounds the subspace recovery error between $\hat{\bm{U}}^{t}$ and ${\bm{U}}^{*}$ uses the overall approach of [25]. After this, the details are different because, for LRMC, the random variables are Bernoulli, while in our case they are unbounded sub-Gaussian or sub-exponential. Also, in LRPR, the sign/phase are unknown and this introduces an extra term that needs to be bounded, see $\mathrm{Term2}$ defined in Lemma 3.9. To bound this we first use Cauchy-Schwarz to bound it by a product of two terms. One term is easy to deal with. For the second term, we borrow a lemma from the RWF paper [10] on standard PR, but the rest of our approach is different because of the need to prove concentration of a set of $mq$ measurements that are not identically distributed (we need to carefully exploit right concerence to ensure that they are “similar enough” to apply the concentration bounds for sums of products of sub-Gaussians).

Our initialization step uses the truncated spectral initialization approach, first introduced for PR in [9]. The proof for it also uses the overall approach of [9] but there are many important differences in proving the concentration bounds (see above). Our initialization is also significantly different from that of LRMC or LRMS. In these cases, the phase is known and thus one can come up with a matrix whose expected value equals ${\bm{X}^{*}}$. This is not possible for LRPR. The matrix whose top $r$ singular vectors we compute has expected value $2\beta_{2}{\bm{X}^{*}}{\bm{X}^{*}}^{\prime}+\beta_{1}\bm{I}$. This point has important implications for the sample complexity dependence on $\kappa,r$, we discuss this next.

II-D2 Worse dependence on $\kappa$, $r$

Our sample complexity is comparable to that of [25] in terms of its dependence on $n$ and $r$. However, our result has a worse dependence on $\kappa$ because we have access only to phaseless measurements. Because of this, (1) our initialization step needs to use the matrix $\bm{Y}_{U}$ and find its top $r$ eigenvectors in order to get an initial estimate of the column span ${\bm{U}}^{*}$. For simplicity, consider $\bm{Y}_{U}$ without the truncation (without the indicator function). Then, its expected value is ${\bm{U}}^{*}{\bm{\Sigma}^{*}}^{2}{\bm{U}}^{*}{}^{\prime}+2\mathrm{trace}({\bm{\Sigma}^{*}}^{2})\bm{I}$. Notice that the condition number of the first term of this matrix (the term of interest) is $\kappa^{2}$. Because of this, when analyzing this step, we end up with a dependence of $mq$ on $\kappa^{8}\mu^{4}nr^{2}/\delta_{\mathrm{init}}^{2}$ (see Claim 3.1 and Sec. III-B where this claim is proved). Here $\delta_{\mathrm{init}}$ is the subspace error bound for the initialization step. Instead, the expected value of the matrix used for initialization of AltMinComplete [25] has expected value equal to ${\bm{X}^{*}}$ and thus its condition number is just $\kappa$. (2) A second issue is as follows: because of magnitude-only measurements, we are having to deal with phase error (sign error) in each LS step that updates the estimate of ${\bm{U}}^{*}$. In bounding the phase error term – $\mathrm{Term2}$ defined in Lemma 3.9 in Sec. III-C – we need to use the Cauchy-Schwarz inequality; see proof of Lemma 3.12. Because of this, when using the bound on $\mathrm{Term2}$ from Lemma 3.12 to prove the main descent claim, Claim 3.4, we end up with a bound of the form $\sin\Theta({\bm{U}}^{t+1},{\bm{U}}^{*})\leq C\delta_{t}(\epsilon_{1}+\sqrt{\delta_{t}+\epsilon_{2}})\sqrt{r}\kappa$ where $\delta_{t}$ is the bound on the subspace error from the previous step and $\epsilon_{i}$ are quantities used in the concentration bounds for Term1 and Term2 respectively. Thus, to ensure that the $(t+1)$-th step error is below $0.7\delta_{t}$ (decays geometrically), we need $\delta_{t}\leq c/(\kappa^{2}r)$ for each iteration $t$. Since we show $\delta_{t}\leq\delta_{\mathrm{init}}$ for all $t$, this is ensured if $\delta_{\mathrm{init}}=c/(\kappa^{2}r)$. We also need $\epsilon_{1}\leq c/kappa\sqrt{r}$ and $\epsilon_{2}\leq c/(\kappa^{2}r)$. This, along with the fact that the initialization step needs $mq\geq\kappa^{8}\mu^{4}nr^{2}/\delta_{\mathrm{init}}^{2}$, implies that our sample complexity per iteration becomes $C\kappa^{12}\mu^{4}nr^{4}$.

If Cauchy-Schwarz were not used, we would only need $mq\geq C\kappa^{10}\mu^{4}r^{3}$. Similarly, if we somehow did not use the loose bound $\|\bm{X}\|\leq\|\bm{X}\|_{F}$ for rank $r$ matrices at a few different places in the proof, we could remove another factor of $r$.

II-D3 Our previous work

The only other work that also studies our problem is our previous work [20]. This introduced a series of heuristics and evaluated them experimentally. It also provided a guarantee for the initialization step of one of them. If we compare their main result (their Theorem 3.2) with ours, it required the following lower bound on just $m$: $m\geq C\max(\sqrt{n},r^{4})/\epsilon^{2}$ in addition to a lower bound on $mq$ that also depends on $1/\epsilon^{2}$. We remove the $1/\epsilon^{2}$ dependence by analyzing the complete algorithm.

The first two requirements on just $m$ are also significantly relaxed in our work because we study a significantly modified version of our previous algorithm. The most important algorithmic difference is that, both for initialization and for later iterations, we recover $\tilde{\bm{b}}^{*}_{k}$’s by solving the standard PR problem either fully (or, at least for enough iterations so that the error in recovering $\tilde{\bm{b}}^{*}_{k}$’s is of the same order as the subspace error in the estimate of ${\bm{U}}^{*}$). This is what allows us to replace the strong requirement $m\geq C\sqrt{n}$ that [20] needed by just $m\geq Cr$. This is also what enables us to get a complete guarantee for the entire algorithm. The algorithm in [20] used only one iteration of AltMinPhase [6] for obtaining a new estimate of $\tilde{\bm{b}}^{*}_{k}$’s using a new estimate of ${\bm{U}}^{*}$. With this, it was not possible to show that the recovery error of $\tilde{\bm{b}}^{*}_{k}$’s is of the same order as that of ${\bm{U}}^{*}$.

Our approach for initializing ${\bm{U}}^{*}$ is taken from [20], but with a simple, but important, difference: the threshold in the indicator function used for defining $\bm{Y}_{U}$ in (6) now takes an average over all $mq$ measurements (instead of over only the $m$ measurements of the $k$-th column in [20]). This simple change allows us to use concentration over all the $mq$ measurements (and design vectors) in every step of deriving the initialization guarantee for ${\bm{U}}^{*}$. This is what helps us eliminate the requirement of $m\geq Cr^{4}$ on just $m$ that was needed in [20].

III Proof of Theorem 2.1

The proof borrows ideas from past works – [9, 20] (for initialization of ${\bm{U}}^{*}$), [25] (the overall approach for getting a subspace error bound given in the Appendix), [52] for careful $\epsilon$-net arguments for unit Frobenius norm matrices, and [10] (for recovering $\tilde{\bm{b}}^{*}_{k}$’s, and in one step of trying to show that the phase error is small). We cite the relevant reference again where it is used. In Sec. III-A next, we provide the two main claims (one for initialization and one for the descent), the two other auxiliary lemmas needed for proving Theorem 2.1 and the theorem’s proof using these. In Sec. III-B, we give the key lemmas needed for proving the initialization claim and also prove it. The same is done for the descent claim in Sec. III-C. Each of these subsections also provides the main ideas (intuition) used for proving the lemmas. We then prove all but one of the lemmas from this entire section in Appendix A. Lemma 3.9 which uses the overall approach of [25] is proved in Appendix B.

III-A Overall lemmas and proof of Theorem 2.1

Claim 3.1** (Rank estimation and Initialization of ${\bm{U}}^{*}$).**

Let ${\bm{U}}_{\mathrm{init}}=\hat{\bm{U}}^{0}$. Pick a $\delta_{\mathrm{init}}<0.25$. Assume $mq\geq\kappa^{8}\mu^{4}nr^{2}/\delta_{\mathrm{init}}^{2}$. Set the rank estimation threshold $\omega=1.3{\sigma_{\min}^{*}}^{2}/q$ (we can actually set the multiplier to any number between $0.025$ and $1.5$). Then, w.p. at least $1-6n^{-10}$, the rank is correctly estimated and

[TABLE]

Define

[TABLE]

It is easy to see that $\bm{x}^{*}_{k}={\bm{U}}^{t}\bm{g}_{k}^{t}+\bm{e}_{k}^{t}$ and so $\bm{y}_{ik}=|(({\bm{U}}^{t}){}^{\prime}\bm{a}_{ik})^{\prime}\bm{g}_{k}^{t}+\bm{a}_{ik}{}^{\prime}\bm{e}_{k}^{t}|$. Thus, we have a noisy PR problem to solve with the noise magnitude proportional to $\|\bm{e}_{k}^{t}{}\|$ with $\|\bm{e}_{k}^{t}{}\|\leq\sin\Theta({\bm{U}}^{t},{\bm{U}}^{*})\|\bm{x}^{*}_{k}\|$. We use RWF to solve it. RWF provides an estimate, $\bm{\hat{b}}_{k}^{t}$, of $\bm{g}_{k}^{t}$. Observe that

[TABLE]

is just a rotated version of $\tilde{\bm{b}}^{*}_{k}$. We show in the next lemma that, whp, the error in the RWF estimate, $\mathrm{dist}(\bm{g}_{k}^{t},\bm{\hat{b}}_{k}^{t})$, is proportional to $\sin\Theta({\bm{U}}^{t},{\bm{U}}^{*})$; and the same is true for the error in $\bm{\hat{x}}_{k}^{t}:={\bm{U}}^{t}\bm{\hat{b}}_{k}^{t}$.

Lemma 3.2** (Recovery of $\tilde{\bm{b}}^{*}_{k}$’s).**

At iteration $t$, assume that $\sin\Theta({\bm{U}}^{*},{\bm{U}}^{t})\leq\delta_{t}$. Pick a $\delta_{b}<1$. If $m\geq Cr$, and if we set $T_{RWF,t}=C\log\delta_{t}/\log(1-c)$, then, w.p. at least $1-2q\exp\left(-c\delta_{b}^{2}m\right)$, the following is true for each $k=1,2,\cdots,q$

[TABLE]

with $C=\sqrt{1+\delta_{b}}+1$.

Thus, if $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$, then the above bounds hold w.p. at least $1-n^{-10}$.

From above, $\bm{\hat{b}}_{k}^{t}$ is close to $\bm{g}_{k}^{t}$ (which is a rotated version of $\tilde{\bm{b}}^{*}_{k}$) for each $k$. We thus expect $\hat{\bm{B}}^{t}$, or equivalently $\bm{B}^{t}$, to also satisfy the incoherence assumption. We show next that this is indeed true if $\delta_{t}$ is small enough. Recall that $\hat{\bm{B}}^{t}\overset{\mathrm{QR}}{=}{\bm{R}}_{B}\bm{B}^{t}$.

Lemma 3.3** (Incoherence of ${\bm{B}^{*}}$ implies incoherence of $\bm{B}^{t}$).**

Pick a $\delta_{b}<1/10$ and assume that $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$. At iteration $t$, assume that $\sin\Theta({\bm{U}}^{*},{\bm{U}}^{t})\leq\delta_{t}$ with $\delta_{t}\leq\frac{0.25}{C\sqrt{r}\kappa}$. If ${\bm{B}^{*}}$ is $\mu$-incoherent, then, w.p. at least $1-n^{-10}$, $\bm{B}^{t}$ is $\hat{\mu}$-incoherent with $\hat{\mu}=C\kappa\mu$.

Finally, the next claim shows that the LS step to update ${\bm{U}}$ reduces its error by a factor of $0.7$ at each iteration. Its proof relies on the previous two lemmas and the fact that $\bm{\hat{x}}_{k}^{t}$ close to $\bm{x}^{*}_{k}$ implies that, with large probability, the phases (signs) of $(\bm{a}_{ik}{}^{\prime}\bm{\hat{x}}_{k}^{t})$ and $(\bm{a}_{ik}{}^{\prime}\bm{x}^{*}_{k})$ are equal too.

Claim 3.4** (Descent Lemma).**

At iteration $t$, assume that $\sin\Theta({\bm{U}}^{*},{\bm{U}}^{t})\leq\delta_{t}$ and $\delta_{t}\leq\delta_{\mathrm{init}}\leq c/r\kappa^{2}$. If $mq\geq C\kappa^{6}\mu^{2}nr^{3}$ and $m\geq C\max(r,\log n,\log q)$ then w.p. at least $1-C\exp(-nr)-n^{-10}$,

[TABLE]

Proof of Theorem 2.1.

The $\sin\Theta({\bm{U}}^{t},{\bm{U}}^{*})$ bounds are an immediate consequence of Claims 3.1 and 3.4, along with setting $\delta_{\mathrm{init}}=c/\kappa^{2}r$. Claim 3.1 needs $mq\geq C\kappa^{8}\mu^{4}nr^{2}/\delta_{\mathrm{init}}^{2}=C\kappa^{8}\mu^{4}nr^{4}$ along with $m\geq C\max(r,\log n,\log q)$ while Claim 3.4 needs $mq\geq C\kappa^{6}\mu^{2}nr^{3}$ and $m\geq C\max(r,\log n,\log q)$. The first lower bound on $mq$ dominates. Thus the required sample complexity is $m_{tot}q\geq C\kappa^{8}\mu^{4}\cdot nr^{4}\log(1/\epsilon)$ and $m_{tot}\geq C\max(r,\log n,\log q)\log(1/\epsilon)$.

The other bounds of Theorem 2.1 follow by Lemma 3.2.∎

We prove Claims 3.1 and 3.4 next in Sec. III-B and III-C. The proof of Lemmas 3.2 and 3.3 and of the lemmas needed for proving these two claims is postponed to Appendix A.

III-B Proof of Claim 3.1

In this section, we let $\bm{a}_{ik}:=\bm{a}_{ik}^{(0)}$ and $\bm{y}_{ik}:=\bm{y}_{ik}^{(0)}$.

The overall idea for proving this is inspired by the approach in [20] which itself borrows ideas from [9]. But there are many important differences because we define $\bm{Y}_{U}$ differently in this work, see (6): the threshold in the indicator function now takes an average over all $mq$ measurements (instead of over only the $m$ measurements of the $k$-th column as in [20]). This simple change enables us to get a significantly improved result. It lets us use concentration over all the $mq$ measurements (and design vectors) in each of the three steps of the proof. This is what helps eliminate the lower bound $m\geq Cr^{4}$ on just $m$ that was needed in [20]. However, this also means that the proofs are much more involved (more quantities now vary with $k$).

Recall the expression for $\bm{Y}_{U}$ from earlier, and define matrices $\bm{Y}_{-}(\epsilon_{1})$ and $\bm{Y}_{+}(\epsilon_{1})$ as

[TABLE]

Define $\bm{Y}_{+}(\epsilon_{1})$ similarly but with $(1-\epsilon_{1})$ replaced by $(1+\epsilon_{1})$ in the indicator function. We will show that $\bm{Y}_{U}$ is sandwiched between $\bm{Y}_{-}$ and $\bm{Y}_{+}$. This, along with showing that $\bm{Y}_{-}$ and $\bm{Y}_{+}$ are close, will help us show that $\bm{Y}_{U}$ is close to $\bm{Y}_{-}$ and, hence, also to its expected value. After this, use of the $\sin\theta$ theorem will give us the desired bound.

Adapting the approach of [9, 20],

[TABLE]

where

[TABLE]

and $\xi$ is a scalar standard Gaussian random variable. The expression for $\mathbb{E}\left[\bm{Y}_{+}(\epsilon_{1})\right]$ is similar but with $(1-\epsilon_{1})$ replaced by $(1+\epsilon_{1})$ in the expression for $\beta_{1,k}^{+}$, $\beta_{2,k}^{+}$.

Observe that $\mathbb{E}\left[\bm{Y}_{-}(\epsilon_{1})\right]$ can be simplified as

[TABLE]

Thus, the span of its top $r$ eigenvectors (same as singular vectors) equals $\mathrm{Span}({\bm{U}}^{*})$. Hence, we can use the $\sin\Theta$ theorem [53] stated below in a fashion similar to [20] (Sec 6).

Lemma 3.5** (Davis-Kahan $\sin\Theta$ theorem).**

Given two symmetric matrices $\bm{D}$ and $\hat{\bm{D}}$. Let ${\bm{U}}^{*}$ (${\bm{U}}$) be the matrix of top eigenvectors of $\bm{D}$ ($\hat{\bm{D}}$). If $\lambda_{r}(\bm{D})-\lambda_{r+1}(\bm{D})-\|\bm{D}-\hat{\bm{D}}\|>0$, then

[TABLE]

Using Lemma 3.5 with $\hat{\bm{D}}=\bm{Y}_{U}$ and $\bm{D}=\mathbb{E}\left[\bm{Y}_{-}\right]$,

[TABLE]

Moreover,

[TABLE]

Now we just need to upper bound $\|\bm{Y}_{U}-\mathbb{E}\left[\bm{Y}_{-}(\epsilon_{1})\right]\|$ and lower bound $\min_{k}\beta_{1,k}^{-}$. Both these follow by combining the three lemmas given next and triangle inequality.

Lemma 3.6**.**

We have that, w.p. at least $1-\exp(-\epsilon_{1}^{2}\frac{mq}{\mu^{2}\kappa^{2}})$,

[TABLE]

and so $\|\bm{Y}_{U}-\bm{Y}_{-}(\epsilon_{1})\|\leq\|\bm{Y}_{+}(\epsilon_{1})-\bm{Y}_{-}(\epsilon_{1})\|$.

Lemma 3.7**.**

Let $\bm{Y}_{+}=\bm{Y}_{+}(\epsilon_{1})$ and $\bm{Y}_{-}=\bm{Y}_{-}(\epsilon_{1})$. We have

[TABLE]

and, assuming $\epsilon_{1}<0.01$,

[TABLE]

Lemma 3.8**.**

*We have that,

w.p. at least $1-2\exp\left(n\log 9-c\epsilon_{2}^{2}mq\right)$,*

[TABLE]

We get the exact same claim also for $\|\bm{Y}_{+}-\mathbb{E}\left[\bm{Y}_{+}\right]\|$.

We prove the above lemmas in Sec. A-B.

Using triangle inequality, $\|\bm{Y}_{U}-\mathbb{E}[\bm{Y}_{-}]\|\leq\|\bm{Y}_{U}-\bm{Y}_{-}\|+\|\bm{Y}_{-}-\mathbb{E}[\bm{Y}_{-}]\|$. Moreover, using Lemma 3.6, $\|\bm{Y}_{U}-\bm{Y}_{-}\|\leq\|\bm{Y}_{+}-\bm{Y}_{-}\|$. Using these and again using triangle inequality, $\|\bm{Y}_{U}-\mathbb{E}[\bm{Y}_{-}]\|\leq 2\|\bm{Y}_{-}-\mathbb{E}[\bm{Y}_{-}]\|+\|\bm{Y}_{+}-\mathbb{E}[\bm{Y}_{+}]\|+\|\mathbb{E}[\bm{Y}_{+}]-\mathbb{E}[\bm{Y}_{-}]\|$. Thus, combining bounds from the above lemmas and setting $\epsilon_{1}=\epsilon_{2}=\frac{\delta_{\mathrm{init}}}{C(\kappa^{2}\mu^{2})\kappa^{2}r}$ for a $\delta_{\mathrm{init}}<1$, we conclude that

w.p. $1-2\exp\left(n-\frac{c\delta_{\mathrm{init}}^{2}mq}{\kappa^{8}\mu^{4}r^{2}}\right)-2\exp\left(-\frac{c\delta_{\mathrm{init}}^{2}mq}{\kappa^{8}\mu^{4}r^{2}}\right)$,

[TABLE]

Using $\min_{k}\beta_{1,k}^{-}\geq 1.5$, from Lemma 3.7, and (10), since $\delta_{\mathrm{init}}<1$,

[TABLE]

III-B1 Proof that rank is correctly estimated

Consider the rank estimation step. This requires lower bounding $\lambda_{r}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})$ and upper bounding $\lambda_{r+1}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})$. Both bounds follow using (i) (10), along with Weyl’s inequality, and (ii) the lower bound on $\beta^{-}:=\min_{k}\beta_{1,k}^{-}$ from Lemma 3.7: $\beta^{-}>1.5$.

We have $\lambda_{r}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})\geq\lambda_{r}(\mathbb{E}\left[\bm{Y}_{-}\right])-\lambda_{n}(\mathbb{E}\left[\bm{Y}_{-}\right])-2\|\bm{Y}_{-}-\mathbb{E}\left[\bm{Y}_{-}\right]\|\geq\beta^{-}\lambda_{r}({\bm{U}}^{*}{\Sigma^{*}}^{2}{{\bm{U}}^{*}}^{\prime})-2\|\bm{Y}_{U}-\mathbb{E}\left[\bm{Y}_{-}\right]\|\geq 1.5({\sigma_{\min}^{*}}^{2}/q)-0.25\delta_{\mathrm{init}}{\sigma_{\min}^{*}}^{2}/q=(1.5-0.25\delta_{\mathrm{init}}){\sigma_{\min}^{*}}^{2}/q>1.4({\sigma_{\min}^{*}}^{2}/q)$ as long as $\delta_{\mathrm{init}}<0.1$.

Also, $\lambda_{r+1}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})\leq\lambda_{r+1}(\mathbb{E}\left[\bm{Y}_{-}\right])-\lambda_{n}(\mathbb{E}\left[\bm{Y}_{-}\right])+2\|\bm{Y}_{U}-\mathbb{E}\left[\bm{Y}_{-}\right]\|=2\|\bm{Y}_{-}-\mathbb{E}\left[\bm{Y}_{-}\right]\|\leq 0.25\delta_{\mathrm{init}}({\sigma_{\min}^{*}}^{2}/q)<0.025({\sigma_{\min}^{*}}^{2}/q)$ as long as $\delta_{\mathrm{init}}<0.1$.

In summary, as long as $\delta_{\mathrm{init}}<0.1$, $\lambda_{r}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})\geq 1.4{\sigma_{\min}^{*}}^{2}/q$ and $\lambda_{r+1}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})\leq 0.025{\sigma_{\min}^{*}}^{2}/q$. Thus by setting the threshold $\omega=C{\sigma_{\min}^{*}}^{2}/q$ with $C$ being any constant between 0.025 and 1.4, we can ensure that the rank is correctly estimated whp.

III-C Proof of Claim 3.4

In this section, we remove the superscript t except where essential. Also, we let $\bm{a}_{ik}:=\bm{a}_{ik}^{(T+t)}$ and $\bm{y}_{ik}:=\bm{y}_{ik}^{(T+t)}$.

We first use the overall approach of [25] to get the following deterministic bound on the subspace error of the $(t+1)$-th estimate of ${\bm{U}}^{*}$, ${\bm{U}}^{t+1}$. The proof requires some messy algebra and hence we give it in Appendix B.

Lemma 3.9**.**

We have

[TABLE]

where $\mathrm{MainTerm}:=$

[TABLE]

[TABLE]

[TABLE]

is the space of all $n\times r$ matrices with unit Frobenius norm, and $\bm{c}_{ik},\hat{\bm{c}}_{ik}$ are the phases (signs) of $\bm{a}_{ik}{}^{\prime}\bm{x}^{*}_{k}$ and $\bm{a}_{ik}{}^{\prime}\bm{\hat{x}}_{k}$.

We obtain high probability bounds on the three terms above in the three lemmas that follow, Lemmas 3.10, 3.11, 3.12. All three lemmas first bound the terms for a fixed $\bm{W}$, followed by using a carefully developed epsilon-net argument to extend the bounds for all unit Frobenius norm $\bm{W}$’s. This is inspired by similar arguments in [52].

Consider a fixed $\bm{W}$. To bound $\mathrm{Term1}$, we first show $\mathbb{E}[\mathrm{Term1}]=0$. Next, we use Lemmas 3.2 and 3.3 to show that $\|{\bm{X}^{*}}\bm{B}^{\prime}\bm{b}_{k}-\bm{x}^{*}_{k}\|\leq C\delta_{t}\|{\bm{X}^{*}}\|_{F}\max(\|\bm{b}^{*}_{k}\|,\|\bm{b}_{k}\|)\leq C\delta_{t}\|{\bm{X}^{*}}\|_{F}\hat{\mu}\sqrt{r/q}$ and $\|{\bm{X}^{*}}(\bm{B}^{\prime}\bm{B}-\bm{I})\|\leq C\delta_{t}\|{\bm{X}^{*}}\|_{F}$. Finally, we use these two facts and the concentration bound of Lemma A.1 (bounds sums of products of sub-Gaussian random variables) to show that, if $mq$ is large enough, whp, $|\mathrm{Term1}|\leq Cm\epsilon_{1}\delta_{t}\|{\bm{X}^{*}}\|_{F}$ for any $\delta_{t}<0.1$. This is followed by a careful epsilon-net argument to extend the bound for all unit Frobenius norm $\bm{W}$’s.

To bound $\mathrm{Term2}$ for a fixed $\bm{W}$, we first use Cauchy-Schwarz. This implies that

[TABLE]

We explain how to upper bound $\mathrm{Term3}(\bm{W})$ in the next paragraph. Consider $\mathrm{Term22}$. Notice that $(\bm{c}_{ik}\hat{\bm{c}}_{ik}-1)^{2}$ takes only two values - zero or four. It is zero when the signs are equal, else it is four. To start bounding $\mathbb{E}[\mathrm{Term22}]$, we can use Lemma 1 of [10]. This shows that the probability that the signs are unequal is upper bounded by a term that is directly proportional to the ratio $\mathrm{dist}^{2}(\bm{x}^{*}_{k},\bm{\hat{x}}_{k})/(\bm{a}_{ik}{}^{\prime}\bm{x}^{*}_{k})^{2}$. The probability bound is thus large when this ratio is large and small otherwise. Moreover, it is easy to see that $\mathbb{E}[(\bm{a}_{ik}{}^{\prime}\bm{x}^{*}_{k})^{2}]=\|\bm{x}^{*}_{k}\|^{2}$ and by Lemma 3.2, whp, $\mathrm{dist}^{2}(\bm{x}^{*}_{k},\bm{\hat{x}}_{k})\leq\delta_{t}^{2}\|\bm{x}^{*}_{k}\|^{2}$. Relying on these ideas, we can argue that, on average, $\mathrm{Term22}$ is very small: $\mathbb{E}[\mathrm{Term22}]\leq Cm\delta_{t}^{3}\|{\bm{X}^{*}}\|_{F}^{2}$. Use of concentration bound of Lemma A.1 then implies that, if $mq$ is large enough, $\mathrm{Term22}$ is bounded by $Cm(\epsilon_{2}+\delta_{t})\delta_{t}^{2}\|{\bm{X}^{*}}\|_{F}^{2}$ whp.

To upper and lower bound $\mathrm{Term3}$, notice first that $\mathbb{E}[\mathrm{Term3}]=m\|\bm{W}\bm{B}\|_{F}^{2}=m$. Also, each summand in this term is sub-exponential with sub-exponential norm bounded by $\|\bm{W}\bm{b}_{k}\|^{2}$; and $\|\bm{W}\bm{b}_{k}\|^{2}\leq\|\bm{b}_{k}\|^{2}\leq\hat{\mu}^{2}r/q$ (by Lemma 3.3), and $\sum_{ik}\|\bm{W}\bm{b}_{k}\|^{2}=m$. Using these facts and Lemma A.1, we can show that $\mathrm{Term3}$ concentrates around $m$ whp.

Lemma 3.10**.**

Pick a $\delta_{b}<1/10$ and assume that $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$. Under the conditions of Theorem 2.1, for a $\delta_{t}<1/10$, w.p. at least $1-2\exp\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)-n^{-10}$,

[TABLE]

and

[TABLE]

Lemma 3.11**.**

Pick a $\delta_{b}<1/10$ and assume that $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$. Under the conditions of Theorem 2.1 and assuming that $\sin\Theta({\bm{U}}^{*},{\bm{U}})\leq\delta_{t}$, with $\delta_{t}<1/10$, w.p. at least $1-2\exp\left(nr(\log 17)-c\frac{\epsilon_{1}^{2}mq}{\kappa^{3}\mu^{2}r}\right)-n^{-10}$.

[TABLE]

In proving the above, we also show that

[TABLE]

Lemma 3.12**.**

Pick a $\delta_{b}<1/10$ and assume that $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$. Under the conditions of Theorem 2.1 and assuming $\sin\Theta({\bm{U}}^{*},{\bm{U}})\leq\delta_{t}$ with $\delta_{t}<1/10$, w.p. at least $1-2\exp\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)-2\exp\left(-c\epsilon_{2}^{2}mq\right)-n^{-10}$,

[TABLE]

Finally, we lower bound the first denominator term of (11). This can be done by using the bound on $\|\tilde{\bm{B}}^{*}\left(\bm{I}-\bm{B}{}^{\prime}\bm{B}\right)\|_{F}$ from Lemma 3.11. This, in turn, implies a lower bound on the minimum singular value of $\bm{B}^{*}\bm{B}^{\prime}$, and hence the following.

Lemma 3.13**.**

Pick a $\delta_{b}<1/10$ and assume that $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$. Under the conditions of Theorem 2.1, if $\sin\Theta({\bm{U}}^{*},{\bm{U}})\leq\delta_{t}$ with $\delta_{t}\leq\frac{1}{4C\sqrt{r}\kappa}$, then, w.p. at least $1-n^{-10}$, $\sigma_{\min}({\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}\bm{B}^{\prime})\geq 0.9{\sigma_{\min}^{*}}.$

We prove the above lemmas in Appendix A-E.

Proof of Claim 3.4.

Combining Lemmas 3.9 and 3.13, if $\delta_{t}<c/\sqrt{r}\kappa$,

[TABLE]

Set $\delta_{b}=1/11$. Combining Lemmas 3.11, 3.12 and 3.10, and using $\|{\bm{X}^{*}}\|_{F}\leq\sqrt{r}{\sigma_{\max}^{*}}$, we conclude that, w.p. at least

$1-2\exp\left(nr(\log 17)-c\frac{\epsilon_{1}^{2}mq}{\kappa^{3}\mu^{2}r}\right)-2\exp\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)-2\exp\left(-c\epsilon_{2}^{2}mq\right)-n^{-10}$

[TABLE]

In order to ensure that $\mathrm{MainTerm}\leq 0.7\delta_{t}{\sigma_{\min}^{*}}$, we need to set $\epsilon_{3}=0.1$, $\epsilon_{1}=c/\sqrt{r}\kappa$, $\sqrt{\epsilon_{2}}=c/\sqrt{r}\kappa$, and $\sqrt{\delta_{t}}\leq c/\sqrt{r}\kappa$. Since we prove $\delta_{t}\leq\delta_{\mathrm{init}}$, the last bound is ensured if $\sqrt{\delta_{\mathrm{init}}}\leq c/\sqrt{r}\kappa$. With these settings, if $mq\geq Cnr^{2}/\epsilon_{1}^{2}=Cnr^{3}$, $mq\geq Cnr^{2}/\epsilon_{3}^{2}=Cnr^{2}$, $mq\geq Cnr/\epsilon_{2}^{2}=Cnr^{3}$, and $\delta_{\mathrm{init}}=c/r\kappa^{2}$, then

w.p. $1-n^{-10}-4\exp(-cnr)$,

[TABLE]

∎

IV Numerical Evaluation

In this section we provide detailed description of the numerical evaluation of our algorithms on synthetic and real data. All time comparisons are performed on a single Desktop Computer with Intel${}^{\textsuperscript{\textregistered}}$ Xeon E$3$-$1240$ $8$-core CPU @ $3.50$GHz and $32$GB RAM. We must mention that for all the algorithms, (a) we compare on the same system, and we do not run other memory and compute intensive programs while performing the time comparison, (b) we use the most efficient sub-routines that are provided by the respective authors and these implementations rely on the “anonymous functions” feature of MATLAB. This provides a significant speed up since the major chunk of computation time involves matrix-vector/matrix-matrix products. This ensures uniformity to a large extent for all algorithms.

IV-A Synthetic Data Experiments

We demonstrate the effectiveness of AltMinLowRaP over existing work using two synthetic experiments. For both, we generate an error versus time-taken plot as follows: for each $t=0,1,\dots,T$, we plot the matrix recovery error (at the end of that iteration) and time-taken (until the end of that iteration) on the y- and x-axes respectively. All synthetic data experiments are performed for $100$ independent trials, and for each algorithm, we plot the average error over the best $90$ trials (drop the 10 trials with the largest error). This is done because all algorithms are guaranteed to work well only with high probability.

In the first experiment, we compare the performance of several algorithms for different values of $m$. We generated data as follows: ${\bm{X}^{*}}={\bm{U}}^{*}\tilde{\bm{B}}^{*}$ where ${\bm{U}}^{*}\in\mathbb{R}^{n\times r}$ is generated by orthonormalizing a iid standard Gaussian matrix. The entries of $\tilde{\bm{B}}^{*}\in\mathbb{R}^{r\times q}$ are chosen from another iid standard Gaussian distribution. Thus, in this setting, ${\sigma_{\min}^{*}}^{2}/q\approx{\sigma_{\max}^{*}}^{2}/q\approx 1$. Measurements were generated using (1) with $\bm{a}_{ik}$’s being iid Gaussian. We compare AltMinLowRaP with LRPR2 (best heuristic from [20]), RWF [10], and with what we call projected-RWF or proj-RWF (for all $t$, after the $t$-th RWF iteration, we project the matrix $\hat{X}^{t}_{RWF}$ onto the space of rank-$r$ matrices, with $r$ known). We provide results for three cases: (i) $n=200$, $q=400$, $r=4$, and $m=0.4n$; (ii) $n=600$, $q=1000$, $r=4$, and $m=0.25n$; and (iii) $n=200$, $m=5n=1000$, $q=200$ and $r=4$. The results are summarized in Fig. 3. In the first two settings (Figs. 3(a),(b)), we show that by exploiting the low-rank structure, AltMinLowRaP is able to outperform the unstructured PR methods. Notice that, proj-RWF performs better than RWF (which does not work at all since $m\ll n$). Secondly, AltMinLowRaP is significantly better than proj-RWF, and for reasons explained earlier, it is also better than LRPR2. AltMinLowRaP and LRPR2 estimate the rank, whereas proj-RWF is provided the true rank. We observed that, in practice, the estimated rank is higher than the true rank in all the cases. In both these settings, $m\ll n$, and thus, unstructured algorithms fail. The third setting is a setting with a very large value of $m$: we are using $m=5n$ measurements. This is a case where all of AltMinLowRaP, TWF, and RWF work, but AltMinLowRaP is much slower as shown in Fig. 3(c). Additionally, for the first two cases, based on a reviewer’s comment, we also empirically evaluate a random initialization scheme for the AltMinLowRaP algorithm. We observed that although the initial estimates are approximately in the order of $10^{-2}$, the algorithm itself fails to improve this with subsequent iterations. A possible workaround to this could be through using a gradient descent based algorithm (inspired by [12]) but this requires a detailed analysis, and we will study this as part of future work..

Our second experiment illustrates the time complexity discussion given in Sec. II-B. For a given $m$, AltMinLowRaP is about $r$ times slower than the best provably correct regular (unstructured) PR methods - TWF and RWF. But, if for each algorithm, we use the minimum $m$ needed for the algorithm to achieve $\epsilon$ accuracy, then, theoretically, AltMinLowRaP should be faster if the rank $r$ is small enough. We tested this empirically as follows. We generated data as before with $n=600$, $q=1000$ and $r=4$. We implemented TWF and RWF using two values of $m$, $m=3n$ and $m=4n$. We evaluated AltMinLowRaP with using $m=n/4$. The error-at-iteration-$t$ versus time-taken-until-iteration-$t$ plot is shown in Fig. 2 for all these cases. As can be seen, using $3n$ measurements, neither of TWF or RWF works. Using $m=4n$, both work. But if we compare the time taken (x-axis value) for any value of error level $\epsilon$, both are at least 5 times slower than AltMinLowRaP ($m=n/4$). For all algorithms, we repeat the expriments for $100$ independent trials, and plot the mean taken over the best $90$ trials to illustrate the high probability results.

IV-A1 Algorithm parameters

For both experiments, AltMinLowRaP was implemented as Algorithm 1 but using the same set of measurements (does not require sample-splitting), and with the following parameters: $T_{RWF,t}$ scales linearly from 5 to 30, $\omega=1.3{\sigma_{\min}^{*}}^{2}/q\approx 1.3$ and $C_{Y}=9$. All parameters are as suggested in the theorem. For LRPR2 we used the default parameters mentioned in the documentation. We set the maximum number of outer-loop iterations, $T_{max}=10$ for both. For RWF and TWF, we used the default parameters suggested by the authors with the exception that we let the maximum number of iterations $T_{max}=300$ (to try to see if its error reduces with more iterations). However, as we plot the time-taken at the end of each iteration, this is not an unfair implementation of RWF; it only means that we have 300 data points to plot on our graph. We should point out that the recovery error for AltMinLowRaP is sensitive to the choice of $\omega$. If $\omega$ is too small (for given values of $m$ and of $C_{Y}$ used in computing $\bm{Y}_{U}$), the algorithm will significantly over-estimate the rank. This is especially problematic when $m$ is small (or $C_{Y}$ is large for a given $m$). Thus, as a thumb rule, for a lower value of $m/n$, the threshold $\omega$ should be larger. Of course if it is too large, it will underestimate the rank666This is a bigger problem when $\kappa=1$ as in the simulated data above. It is a lesser problem for real approximately low-rank data (e.g., slow changing videos) with a larger $\kappa$, since in those cases, the missed directions will be the ones with smaller singular values.

IV-B Real Videos with Simulated CDP measurements: Small $r$ suffices

We demonstrate the effectiveness of AltMinLowRaP for recovering two real video sequences (these are only approximately low-rank) from simulated Coded Diffraction Pattern (CDP) measurements. These measurements can be represented as $\bm{Y}=|\mathcal{F}(\bm{D}{\bm{X}^{*}})|$ where $\mathcal{F}(\cdot)$ is the DFT operation and the matrix $\bm{D}$ represents a diagonal mask matrix whose diagonal entries are chosen uniformly at random from $\{\pm 1,\pm\sqrt{-1}\}$ and modulate the intensity of the input. We generate CDP measurements of each frame of the video (the $k$-th frame vectorized is $\bm{x}^{*}_{k}$). We compared our algorithm with LRPR2 and RWF. We present the quantitative results in Table II and the visual comparisons in Fig. 1 (given in the beginning), and in Fig. 4. Notice that, in this case, even with $m=5n$ measurements, RWF is unable to accurately recover the video and AltMinLowRaP has a slightly better performance w.r.t. LRPR2. The algorithm parameters are set as in the synthetic data experiments, swith the exception that we now set $T_{max}=30$ for ALtMinLowRaP and LRPR2. AltMinLowRaP implementation used all the speed-up ideas for Fourier measurements explained in [20] for LRPR2 and so did LRPR2.

We tested AltMinLowRaP with three possible values of rank, $r=15,20,25$. As can be seen, even $r=15$ suffices to get a significantly better reconstruction error than RWF for $m=5n$ CDP measurements. For the plane video, $n=6912$ and $q=105$ and for the mouse video, $n=5182$ and $q=90$, and thus in both cases, $r\approx 0.003n$. If $q$ is larger, a natural idea would be to use similar parameter settings, but instead implement the tracking variant of AltMinLowRaP (Algorithm 2).

IV-C Real Videos with Simulated CDP measurements: Low-Rank versus (Wavelet) Sparse Models

To justify the low-rank assumption on videos, we compare with CoPRAM [18], a state-of-the-art, provable algorithm for compressive phase retrieval. Since the videos are not sparse in the spatial domain, as suggested in [18], we use the Haar wavelet as the sparsifying basis777We also experimented with the Daubhechies-3 wavelet as the sparsifying basis in our experiments. However, we noticed that, for the plane video, irrespective of the wavelet basis, the number of coefficients necessary to preserve $\approx 90\%$ of the energy of the video required $s\approx 0.1n$ and thus the choice of wavelet basis is not detrimental in this experiment..

As can be seen from Table III and Fig. 5, the low-rank prior gives a much better reconstruction error in all three cases in this table including the exact sparse case. Since the video is not exactly wavelet sparse, we also performed a comparison on the sparsified video, wherein, for each image frame, we truncate the wavelet coefficients such that approximately $90\%$ of the energy in each frame is preserved. We refer to this as the sparsified video in this experiment. Sparsifying the video significantly improves the performance of CoPRAM, but AltMinLowRaP ($r=15$) is still better. Even with $m=10n$ measurements, the sparse model is unable to capture the finer details in the video. We also observed that a standard complex mask does not work very well for CoPRAM and hence for this experiment, we report the results when the entries of the CDP mask, $\bm{D}$ are chosen uniformly at random from $\{\pm 1\}$. We reshaped each video frame into size $32\times 32$ since the online implementation of their code only works for small sized data. We provide the quantitative results in Table III and qualitative results for the plane video in Fig. 5.

V Phaseless Subspace Tracking

When the matrix ${\bm{X}^{*}}$ consists of a time sequence of signals $\bm{x}^{*}_{k}$, then the column-wise measurements appear one column at a time (sequentially). Hence, there is benefit in trying to develop a mini-batch algorithm that works with measurements of short batches of $\alpha$ consecutive columns. Moreover, for long data sequences, the subspace from which the data are generated could itself change with time. Detecting and being able to track such subspace changes is important for long sequences. Interestingly the algorithm that works for this purpose is a simple modification of the static case idea along with a carefully designed subspace change detection step.

V-A Problem setting

The low-rank assumption is equivalent to assuming that $\bm{x}^{*}_{k}={\bm{U}}^{*}\tilde{\bm{b}}^{*}_{k}$ where ${\bm{U}}^{*}$ specifies a fixed $r$-dimensional subspace. For long signal/image sequences, a better model (one that allows the required subspace dimension $r$ to be smaller) is to let the subspace change with time. As is common in time-series analysis, the simplest model for time-varying quantities is to assume that they are piecewise constant with time. We adopt this approach here. Moreover, in order to easily borrow ideas from the static setting, we will assume that we now have a total of $q_{\mathrm{full}}$ signals (matrix columns) and we will denote the $n\times q_{\mathrm{full}}$ matrix formed by all these columns by ${\bm{X}^{*}}_{\mathrm{full}}$. Our algorithm will operate on measurements of $\alpha$-consecutive-column sub-matrices of ${\bm{X}^{*}}_{\mathrm{full}}$.

Let $k_{0}=1$, and let $k_{j}$ denote the $j$-th subspace change time, for $j=1,2,\dots,J$ and let $k_{J+1}=q_{\mathrm{full}}$. We have the following model

[TABLE]

where ${\bm{U}}^{*}_{{\mathrm{sub}},(j)}$ is an $n\times r$ “basis matrix” for the $j$-th subspace and $\tilde{\bm{d}}^{*}_{k}$ is the coefficients’ vector at time $k$.

The goal is to track the subspaces $\mathrm{Span}({\bm{U}}^{*}_{{\mathrm{sub}},(j)})$ on-the-fly; of course, “on-the-fly” for subspace tracking means with a delay of at least $r$. Once this can be done accurately enough, it is easy to also recover the matrix columns $\bm{x}^{*}_{k}$ (by solving a simple $r$-dimensional PR problem to recover the $\tilde{\bm{d}}^{*}_{k}$’s).

The reason we use a different notation here (the subscript ${\mathrm{sub}}$ and use of $\tilde{\bm{d}}^{*}_{k}$ instead of $\tilde{\bm{b}}^{*}_{k}$) is as follows. Consider an $\alpha$-column sub-matrix formed by $\alpha$ consecutive signals. Let us call it ${\bm{X}^{*}}$ and let ${\bm{X}^{*}}\overset{\mathrm{SVD}}{=}{\bm{U}}^{*}\bm{\Sigma}^{*}{\bm{B}^{*}}$. If all the $\bm{x}^{*}_{k}$’s forming this matrix are generated from the same subspace, say ${\bm{U}}^{*}_{{\mathrm{sub}},(j)}$, then $\mathrm{Span}({\bm{U}}^{*})=\mathrm{Span}({\bm{U}}^{*}_{{\mathrm{sub}},(j)})$ and there is no need for a different notation. However, if a subspace change occurred inside this interval, then we cannot say anything simple like this. All we can say is that ${\bm{X}^{*}}=[{\bm{U}}^{*}_{{\mathrm{sub}},(j-1)}\bm{D}_{(j-1)},{\bm{U}}^{*}_{{\mathrm{sub}},(j)}\bm{D}_{(j)}]$ and so $\mathrm{Span}({\bm{U}}^{*})\subseteq\mathrm{Span}({\bm{U}}^{*}_{{\mathrm{sub}},(j-1)})\cup\mathrm{Span}({\bm{U}}^{*}_{{\mathrm{sub}},(j)})$.

V-B Basic PST algorithm and extensions

As noted earlier the PST algorithm is a simple modification of the static case algorithm (AltMinLowRaP) along with a carefully designed change detection strategy. In the static case, in each iteration, we used a set of $mq$ measurements of a single $n\times q$ matrix ${\bm{X}^{*}}$. For obtaining the guarantees, we assumed a new (independent) set of $2mq$ measurements of the same matrix ${\bm{X}^{*}}$ were used in each iteration ($mq$ for updating the estimate of ${\bm{B}^{*}}$ and another $mq$ for ${\bm{U}}^{*}$). For the tracking setting, using a mini-batch size of $\alpha$, we proceed as follows: each new update iteration uses $2m\alpha$ measurements of a new $\alpha$-consecutive-column sub-matrix of ${\bm{X}^{*}}_{\mathrm{full}}$. The input to the update iteration is the subspace estimate from the previous iteration. Under the assumption that the subspace remains constant for at least $T\alpha$ time instants after a subspace change has been detected, this approach works: with $T=C\log(1/\epsilon)$, we can show that, after $T\alpha$ time instants, we get an $\epsilon$-accurate estimate of the $j$-th subspace.

We summarize the algorithm in Algorithm 2. This toggles between a “detect” and an “update” mode. It starts in the “update” mode (described above) and remains in it for the first $T\alpha$ time instants. At this time it enters the “detect” mode. We are able to guarantee that, when the algorithm enters this detect mode, the previous subspace has been estimated to $\epsilon$ error whp. In the detect mode, the algorithm does not perform any subspace updates. This is done to simplify our analysis; it ensures that, in the interval during which the subspace change occurs, the subspace is not updated. This is what allows us to use our previous two main claims (Claims 3.1 and 3.4) without change to analyze the update mode. Practically, this is of course wasteful. We develop an improvement below.

To understand the change detection strategy, let $\hat{k}_{j}$ denote the estimated change times. Consider an $\alpha$-length interval, $\mathcal{J}_{\alpha}$, contained in $[k_{j},k_{j+1})$. Assume that an $\epsilon$-accurate estimate of the previous subspace ${\bm{U}}^{*}_{{\mathrm{sub}},(j-1)}$ has been obtained by $\hat{k}_{j-1}+T\alpha$ and that this time is before $k_{j}$. Let ${\bm{U}}_{{\mathrm{sub}},(j-1)}$ denote this estimate. Define the matrix

[TABLE]

with $\bm{Y}_{U}=\bm{Y}_{U}(\mathcal{J}_{\alpha})$. This means that $\bm{Y}_{U}$ is as defined earlier in (6) with the $k$ summation being over all $k\in\mathcal{J}_{\alpha}$ (it is using measurements for all the columns within this $\alpha$-length interval). With a little bit of work (see Lemma C.1 and its proof), one can show that, in this interval, the matrix $\bm{Y}_{U,det}:={\bm{U}}_{{\mathrm{sub}},(j-1),\perp}{}^{\prime}\bm{Y}_{U,det,big}{\bm{U}}_{{\mathrm{sub}},(j-1),\perp}{}$ is close to a matrix $\bm{E}_{det}$ whose eigenvalues satisfy

[TABLE]

On the other hand, in an $\alpha$-length interval contained in $[\hat{k}_{j-1}+T\alpha,k_{j})$,

[TABLE]

Thus, this quantity is small when the $j$-th change has not occurred (before $k_{j}$), and is large when the subspace has changed (after $k_{j}$). By using a large enough lower bound on the product $m\alpha$, the same can be shown for the difference between the maximum and minimum eigenvalues of $\bm{Y}_{U,det}$ (these are equal to the maximum and $(n-r)$-th eigenvalues of $\bm{Y}_{U,det,big}$).

Once we have an $\epsilon$-accurate estimate of the current subspace, it is straightforward to also recover the corresponding signals $\bm{x}^{*}_{k}$. This can simply be done by solving a standard PR problem to recover the coefficients vector. See last line of Algorithm 2. This borrows a similar idea from [55].

V-B1 Improved algorithm: PST-all

Notice from Theorem 5.1 that Algorithm 2 can only provably detect and track subspace changes that are larger than a small threshold. While this makes sense for detection, it should be possible to track all types of changes. By including a simple modification in Algorithm 2 (include the “update” step during the detection mode as well), we can empirically demonstrate that this is indeed true. We demonstrate this in Fig 6(a). Moreover, PST-all also removes the other limitation of basic PST (not using the detect phase samples for improving the subspace estimate). Thus, even for large changes that basic PST can detect, PST-all has better tracking performance; see Fig 6(b).

V-C Guarantee for basic PST

We can prove the following about Algorithm 2 (basic PST).

Corollary 5.1** (PST algorithm).**

Consider Algorithm 2. Pick any value of $m\geq C\max(r,\log n,\log q_{\mathrm{full}})$. For this $m$, set $\alpha=\frac{C\kappa^{12}\mu^{4}\cdot nr^{4}}{m}$. Set $T:=C\log(1/\epsilon)$, and the detection threshold $\omega_{det}=c/(\kappa^{2}r)$. Assume that $k_{j+1}-k_{j}\geq(T+3)\alpha$ and that $\sin\Theta({\bm{U}}^{*}_{{\mathrm{sub}},(j-1)},{\bm{U}}^{*}_{{\mathrm{sub}},(j)})^{2}>\frac{2c}{\kappa^{2}r}$. Then, w.p. at least $1-Cn^{-10}$,

we can detect the change with a delay of at most $2\alpha$, while ensuring no false detections, i.e., $k_{j}\leq\hat{k}_{j}\leq k_{j}+2\alpha$; 2. 2.

for any $\epsilon>0$, we can get an $\epsilon$-accurate estimate of the $j$-th subspace with a delay of at most $(T+3)\alpha$ from $k_{j}$ (when the subspace changed); 3. 3.

we have the following subspace error bounds: let ${\bm{U}}_{{\mathrm{sub}},(j)}^{(-1)}={\bm{U}}_{{\mathrm{sub}},(j-1)}:={\bm{U}}_{{\mathrm{sub}},(j-1)}^{(T)}$, and let ${\bm{U}}_{{\mathrm{sub}},(j)}^{(\ell)}$, $\ell=0,1,\dots,T$, be the $\ell$-th estimate;

[TABLE]

*Offline PST returns $\hat{\bm{X}}$ that satisfies $\text{mat-dist}(\hat{\bm{X}},{\bm{X}^{*}})\leq\epsilon$. *

We provide a proof sketch in Appendix C.

The above result shows that, if the subspace remains constant for at least $\alpha\log(1/\epsilon)$ time instants, and if the amount of subspace change (largest principal angle of subspace change) is of order $1/\sqrt{r}$ or larger, then we can both detect the change and track the changed subspace to $\epsilon$ error within a delay of order $\alpha\log 1/\epsilon$. Moreover, for only at most $3\alpha$ time instants after a change, the subspace error does not reduce and is essentially bounded by the amount of change. After this, it decays exponentially every $\alpha$ time instants.

Notice from the expression for $\alpha$ that, if we pick the smallest allowed value of $m$, then the required $\alpha$ (and hence the required delays) will be large. However, we are allowed to tradeoff $m$ and $\alpha$. If we let $m$ grow linearly with $n$, then we will only need $\alpha\approx r^{4}$, which is, in fact, close to the minimum required delay of $r$. This also matches what is seen in existing works on provable subspace tracking (ST) in other settings (e.g., robust ST, ST with missing data, or streaming PCA with missing data) [56, 55, 57]. These are able to allow close to optimal detection and tracking delays but all these assume that $m$ increases linearly with $n$. We can also pick any value of $m$ in between the two extremes of $m=Cr$ or $m=Cn$. For example, if $m=Cn/r$, then $\alpha=r^{5}$ and so on.

V-C1 Related Work

For Phaseless Subspace Tracking (PST) the only works before this work was our first huristic versions [58], and [59]. Other subspace tracking (ST) problems that have been extensively studied include dynamic compressive sensing [60] (a special case of ST where the subspace is defined by the span of a subset of $r$ vectors from a known dictionary matrix), dynamic robust PCA (or robust ST), see [56, 55] and references therein, streaming PCA with missing data [57, 61], and ST with missing data [62, 63, 64, 65, 66]. In terms of works with complete provable guarantees, there is the nearly optimal robust subspace tracking via recursive projected compressive sensing approach [56, 55, 66] and its precursors; recent papers on streaming PCA with missing data [57, 61], and older work on dynamic compressive sensing (CS) [60]. For robust ST, the problem setting itself implies $m=n/2$. In the streaming PCA case, the availability of $m=\rho n$ measurements, with $\rho<1$, is assumed. This is why both achieve close to optimal tracking delays (at least when the added unstructured noise is nearly zero). As noted earlier, our method can also achieve a delay of order $r^{4}$ if we let $m$ grow linearly with $n$.

Dynamic CS (like basic CS) is able to detect support changes (with sufficiently nonzero magnitude) immediately even with a small value of $m=Cr\log n$ measurements; here $r$ is the sparsity level (support size). This is because it is a much simpler special case of ST: in this case, one just needs to be finding the correct subset of basis vectors from a large provided set (dictionary matrix).

V-D Numerical experiments

This experiment evaluates the PST algorithm (Algorithm 2) and PST-all algorithms from Sec. V. We generate the true data for the first subspace ${\bm{X}^{*}}_{0}={\bm{U}}^{*}_{{\mathrm{sub}},(0)}{\bm{D}}_{0}^{*}$ where ${\bm{U}}^{*}_{{\mathrm{sub}},(0)}\in\mathbb{R}^{n\times r}$ with $n=300$, $r=2$ is generated by orthonormalizing the columns of a $n\times r$ iid standard normal matrix. The entries of ${\bm{D}}_{0}^{*}\in\mathbb{R}^{r\times t_{1}}$ with $t_{1}=2992$ are also generated from an i.i.d. standard normal distribution. We generate the true data from the second subspace similarly and set ${\bm{X}^{*}}_{1}={\bm{U}}^{*}_{{\mathrm{sub}},(1)}{\bm{D}}_{1}^{*}$ and we set $q=6000$. Notice that $\kappa\approx 1$. The subspace ${\bm{U}}^{*}_{{\mathrm{sub}},(1)}$ is generated using the idea of [55] as ${\bm{U}}^{*}_{{\mathrm{sub}},(1)}=e^{-\gamma\bm{M}}{\bm{U}}^{*}_{{\mathrm{sub}},(0)}$ in order to control the subspace error. Here $\bm{M}$ is a skew-symmetric matrix and $\gamma$ controls the amount of subspace change. We study two cases in which we set $\gamma=0.08,0.001$ which roughly translates to $\sin\Theta({\bm{U}}^{*}_{{\mathrm{sub}},(0)},{\bm{U}}^{*}_{{\mathrm{sub}},(1)})=0.8,0.01$. We generate the measurement matrices $\bm{A}_{k}(\in\mathbb{R}^{n\times m})\overset{i.i.d.}{\sim}\mathcal{N}(0,\bm{I})$ with $m=100$ for $i=1,\cdots q$. We then implemented PST (Algorithm 2) and PST-all. PST requires large-enough change in order to ensure good results, and PST-all which works even with small changes. We chose the algorithm parameters as follows. We set $\alpha=250$ and $L=8$. For the detection, and initialization steps of both algorithms we set $m_{\mathrm{init}}=m$. We set the threshold for detection, $\omega=0.6$ through cross-validation. The results for the two algorithms are shown in Fig. 6. Notice that for the small change case, since PST is always in the detect mode, it does not improve the estimation error whereas PST-all does. However, when the change is large enough, both algorithms converge to a small error. The results are averaged over $100$ independent trials.

VI Conclusions and Future Work

This work introduced the first simple, fast, and provably correct, algorithm for Low Rank Phase Retrieval – low-rank matrix recovery from different (mutually independent) column-wise phaseless linear projections – AltMinLowRaP. Moreover, since, even the linear version of our problem has not been studied, this work also provides the first fast and provable solution to the linear version, “Compressive” or “Sketched PCA”. AltMinLowRaP relies on a careful spectral initialization followed by alternating minimization. We showed that its required sample complexity is about $r^{3}$ times the order-optimal value of $nr$. We also developed its dynamic extension that is relevant for datasets where we would like to develop a mini-batch solution that recovers the current sub-matrix of ${\bm{X}^{*}}$ without waiting for all the measurements of all the signals (columns) to arrive.

In ongoing work we are (i) exploring how to remove the sample-splitting requirement by studying an alternating gradient descent solution, and attempting to borrow the leave-one-out ideas from [12]; (ii) how to reduce the dependence of our sample complexity on $r$ and on $\kappa$. Some thoughts are provided in Sec. II-D; and (iii) how to analyze an easy modification of AltMinLowRaP to get a better guarantee for the linear version of our problem. The algorithm modification is easy, it just involves replacing the standard PR step for recovering $\tilde{\bm{b}}^{*}_{k}$’s by a simpler LS step, and of course remove the phase/sign estimation step before updating ${\bm{U}}$. In terms of analysis, (a) we can possibly simplify the analysis of the initialization step because in the linear case, $\mathbb{E}[\bm{a}_{i}k\bm{y}_{i}k]=\bm{x}^{*}_{k}$; and (b) in the iterations, there will be no phase error term, $\mathrm{Term2}$, and hence, no need for Cauchy-Schwarz. This latter change itself will reduce the sample complexity to $nr^{3}$ instead of $nr^{4}$.

Open questions for future work include (i) study if we can exploit the right incoherence assumption in the algorithm itself, for example, by using a projected GD approach inspired by [67]; (ii) can phaseless LRMS be solved (this would be the other possible LRPR problem alluded to in the introduction), and (iii) develop a fast algorithm and a guarantee for exploiting both low-rank (as we do) and column-wise sparsity. This type of modeling been used very successfully in the MRI literature to come up with practical algorithms to reduce the sample complexity empirically, see for example, [30, 31]. It has also been studied theoretically in the linear setting [68]. (iv) Another open question is how to analyze the improved tracking algorithm PST-all that we currently only empirically evaluate. In experiments, it is clearly much better than the simpler version we analyze.

Appendix A Proofs of the Lemmas from Sec. III

A-A Simple facts for various proofs

Our proofs will use the following facts: for two arbitrary matrices $\bm{A},\bm{H}$,

$\sigma_{\max}(\bm{A}+\bm{H})\leq\sigma_{\max}(\bm{A})+\|\bm{H}\|$ 2. 2.

$\sigma_{\min}(\bm{A}+\bm{H})\geq\sigma_{\min}(\bm{A})-\|\bm{H}\|$. 3. 3.

$\sigma_{\min}(\bm{A}\bm{H})\geq\sigma_{\min}(\bm{A})\sigma_{\min}(\bm{H})$. 4. 4.

For two basis matrices, ${\bm{U}}_{1},{\bm{U}}_{2}$, $\sigma_{\min}^{2}({\bm{U}}_{1}{}^{\prime}{\bm{U}}_{2})=1-\sin\Theta^{2}({\bm{U}}_{1},{\bm{U}}_{2})$. 5. 5.

For any matrix ${\bm{P}}$, $\|{\bm{P}}\bm{M}\|_{F}\leq\|{\bm{P}}\|\|\bm{M}\|_{F}$ and $\|\bm{M}{\bm{P}}\|_{F}\leq\|{\bm{P}}\|\|\bm{M}\|_{F}$. 6. 6.

For an invertible matrix ${\bm{P}}$,

$\|\bm{M}\|_{F}=\|{\bm{P}}^{-1}{\bm{P}}\bm{M}\|_{F}\leq\|{\bm{P}}^{-1}\|\|{\bm{P}}\bm{M}\|_{F}=\frac{1}{\sigma_{\min}\left({\bm{P}}\right)}\|{\bm{P}}\bm{M}\|_{F}$.

The following lemma is a simple but useful modification of [Lemma 5.16][54]. The proof follows by combining Lemma 2.7.7. and Theorem 2.8.1. of [29]. The sub-Gaussian norm for a random vector $X$, in this lemma, can be defined as follows:

[TABLE]

where $\|z\|_{\psi_{2}}$ is the sub-Gaussian norm of a scalar $z$ [54, Definition 5.7].

Lemma A.1**.**

Let $X_{i},Y_{i}$ be independent sub-Gaussian random variables with sub-Gaussian norm $K_{X_{i}}$ and $K_{Y_{i}}$ respectively and with $\mathbb{E}[X_{i}Y_{i}]=0$. Then

[TABLE]

*When $X_{i}=Y_{i}$, this simplifies to Lemma 5.16 of [54]. *

A-B Proof of the lemmas for Claim 3.1

In this section, we let $\bm{a}_{ik}:=\bm{a}_{ik}^{(0)}$ and $\bm{y}_{ik}:=\bm{y}_{ik}^{(0)}$.

Proof of Lemma 3.6.

Observe that we will be done if we can show that, whp, $\frac{1}{mq}\sum_{ik}\left(\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}\right)^{2}$ lies in the interval $[(1-\epsilon_{1})\|\bm{X}^{*}\|_{F}^{2}/q,(1+\epsilon_{1})\|\bm{X}^{*}\|_{F}^{2}/q]$. Using Lemma A.1, with $K_{X_{ik}}=K_{Y_{ik}}=\|\bm{x}_{k}\|$, with probability more than $1-2\exp\left(\frac{-C\epsilon_{1}^{2}mq}{\mu^{2}\kappa^{2}}\right)$, we have

[TABLE]

Details for obtaining this bound: using $\sum_{k}\|\bm{x}^{*}_{k}\|^{4}\leq\max_{k}\|\bm{x}^{*}_{k}\|^{2}\sum_{k}\|\bm{x}^{*}_{k}\|^{2}$ and right incoherence,

[TABLE]

∎

Proof of Lemma 3.7.

It is easy to see that

[TABLE]

Recall $\gamma_{k}={9\|{\bm{X}^{*}}\|_{F}^{2}\ \mu^{2}\kappa^{2}}/{(q\|\bm{x}^{*}_{k}\|^{2})}$. Using $x^{3}e^{-x^{2}/2}\leq 3\sqrt{3}e^{-3/2}$, we have

[TABLE]

Similarly, using $xe^{-x^{2}/2}\leq e^{-1/2}$,

[TABLE]

Therefore,

[TABLE]

To lower bound $\beta_{1,k}^{-}$, we will use right incoherence which implies that $\gamma_{k}\geq 9$.

[TABLE]

where we used the fact that $x^{3}e^{-x^{2}/4}\leq 3.5$ for any $x$; $\gamma_{k}\geq 9$ (follows by right incoherence); and $\epsilon_{1}<0.01$. ∎

Proof of Lemma 3.8.

Let us define

[TABLE]

As argued in [20], which itself borrows the key idea from [9], we can show that the $\bm{w}_{ik}$s are sub-Gaussian random variables with sub-Gaussian norm $K=C\mu\kappa\|\bm{X}^{*}\|_{F}/\sqrt{q}$. Notice that we have defined $\bm{Y}_{U}$ differently in this paper (in order to be able to exploit concentration over $mq$) as compared to that in [20] and hence only the above argument is similar.

Observe that

$mq\|\bm{Y}_{-}-\mathbb{E}\left[\bm{Y}_{-}\right]\|=\max_{\bm{z}:\|\bm{z}\|=1}|\bm{z}^{\prime}\sum_{ik}(\bm{w}_{ik}\bm{w}_{ik}{}^{\prime}-\mathbb{E}[\bm{w}_{ik}\bm{w}_{ik}{}])\bm{z}|$.

First consider a fixed unit vector $\bm{z}$. Observe that $\bm{z}^{\prime}\bm{w}_{ik}$ is sub-Gaussian with sub-Gaussian norm $K=C\mu\kappa\|\bm{X}^{*}\|_{F}/\sqrt{q}$. Thus, using Lemma A.1 with $t=\epsilon_{2}m\mu^{2}\kappa^{2}\|\bm{X}^{*}\|_{F}^{2}$, and $K_{X_{ik}}=K_{Y_{ik}}=\frac{\mu\kappa\|\bm{X}^{*}\|_{F}}{\sqrt{q}}$, we can conclude that w.p. at least $\geq 1-2\exp\left(-c\epsilon_{2}^{2}mq\right)$,

[TABLE]

After this, we can use a standard epsilon-net argument to extend the bound to all unit vectors $\bm{z}$. With it, we can conclude that, w.p. at least $1-2\exp\left(n\log 9-c\epsilon_{2}^{2}mq\right)$,

[TABLE]

∎

A-C Clarifying the sign inconsistency issue

Recall that we had defined $\bm{g}_{k}^{t}:=({\bm{U}}^{t})^{\prime}\bm{x}^{*}_{k}$ in (7). Since the solution of phase retrieval always comes with a phase (sign) ambiguity, at each iteration $t$, for each $k$, the output of RWF, $\bm{\hat{b}}_{k}^{t}$, may be closer to either $\bm{g}_{k}^{t}$ or $-\bm{g}_{k}^{t}$. This is what decides whether $\mathrm{dist}(\bm{g}_{k}^{t},\bm{\hat{b}}_{k}^{t})$ equals $\|\bm{g}_{k}^{t}-\bm{\hat{b}}_{k}^{t}\|$ or $\|\bm{g}_{k}^{t}+\bm{\hat{b}}_{k}^{t}\|$. However, bound both in each proof is cumbersome. Instead we can proceed as follows. Re-define $\bm{g}_{k}^{t}$ as

[TABLE]

and define the matrix

[TABLE]

With these new definitions, $\mathrm{dist}(\bm{g}_{k}^{t},\bm{\hat{b}}_{k}^{t})=\|\bm{g}_{k}^{t}-\bm{\hat{b}}_{k}^{t}\|$ and $\text{mat-dist}(\bm{G}^{t},\hat{\bm{B}}^{t})=\|\bm{G}^{t}-\hat{\bm{B}}_{t}\|_{F}$.

As an aside, we should point out that, even if some columns of a matrix change sign (are multiplied by $(-1)$), its singular values do not change. Thus, the minimum singular value of $\bm{G}^{t}$ remains the same with or without the above re-definition.

We need to do something similar to the above for $\bm{x}^{*}_{k}$’s as well. Define ${\tilde{\bm{x}}^{*}_{k}}^{t}=\bm{x}^{*}_{k}$ if $\|\bm{\hat{x}}_{k}^{t}-\bm{x}^{*}_{k}\|\leq\|\bm{\hat{x}}^{t}_{k}+\bm{x}^{*}_{k}\|$ and ${\tilde{\bm{x}}^{*}_{k}}^{t}=-\bm{x}^{*}_{k}$ otherwise. Define the corresponding matrix $\tilde{\bm{X}^{*}}^{t}$.

Clearly $\text{mat-dist}(\tilde{\bm{X}^{*}}^{t},{\bm{X}^{*}})=0$. So, in the rest of the writing in this section, to reduce notation, we will re-define

[TABLE]

With this, we can define the error/perturbation in $\bm{\hat{x}}_{k}$ as just

[TABLE]

and we have $\mathrm{dist}(\bm{\hat{x}}_{k}^{t},\bm{x}^{*}_{k})=\|\bm{\hat{x}}_{k}^{t}-\bm{x}^{*}_{k}\|=\|\bm{h}_{k}\|$.

A-D Proof of Lemmas 3.2 and 3.3

In this section, we let $\bm{a}_{ik}:=\bm{a}_{ik}^{(t)}$ and $\bm{y}_{ik}:=\bm{y}_{ik}^{(t)}$. Also, everywhere below, we remove the superscripts t for ease of notation. Recall that $\bm{x}^{*}_{k}={\bm{U}}\bm{g}_{k}+\bm{e}_{k}$ with $\bm{e}_{k}:=(\bm{I}-{\bm{U}}{\bm{U}}^{\prime})\bm{x}^{*}_{k}$.

Proof of Lemma 3.2.

To estimate $\bm{b}_{k}$, we first need to estimate $\bm{g}_{k}$ which requires measurements of the form $\bm{a}_{ik}{}^{\prime}{\bm{U}}\bm{g}_{k}$. Our measurements satisfy

[TABLE]

where $\nu_{ik}=\lvert\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}\rvert-\lvert\bm{a}_{ik}{}^{\prime}{\bm{U}}\bm{g}_{k}\rvert$ is the noise. We use these to obtain the estimate $\bm{\hat{b}}_{k}$ using RWF. By Theorem 2 of [10], if $m\geq Cr$, w.p. at least $1-\exp(-cm)$,

[TABLE]

where $c_{1}$ is a constant less than one. For our problem,

[TABLE]

Clearly, $\mathbb{E}[\|\bm{\nu}_{k}\|^{2}]=m\|\bm{e}_{k}\|^{2}$ and $\|\bm{e}_{k}\|^{2}\leq\sin\Theta({\bm{U}},{\bm{U}}^{*})\|\tilde{\bm{b}}^{*}_{k}\|$. Using Lemma A.1 with $t=m\delta_{b}\|\bm{e}_{k}\|^{2}$, $K_{X_{i}}=K_{Y_{i}}=\|\bm{e}_{k}\|$, and summing over $i=1,2,\dots,m$, we conclude that, w.p. at least $1-\exp\left(-c\delta_{b}^{2}m\right)$,

[TABLE]

where the last inequality used $\|\bm{e}_{k}\|\leq\delta_{t}\|\tilde{\bm{b}}^{*}_{k}\|$. Thus, using the above and $\|\bm{g}_{k}\|\leq\|\tilde{\bm{b}}^{*}_{k}\|$,

[TABLE]

By setting $T_{RWF,t}$ so that $(1-c)^{T_{RWF,t}}\leq\delta_{t}$, we get that $\mathrm{dist}\left(\bm{g}_{k},\bm{\hat{b}}_{k}\right)\leq C\delta_{t}\|\tilde{\bm{b}}_{k}^{*}\|=C\delta_{t}\|\bm{x}^{*}\|$ with $C=(\sqrt{1+\delta_{b}}+1)$. The above bound holds w.p. at least $1-\exp\left(-c\delta_{b}^{2}m\right)$ for a given $k$. By union bound, it holds for all $k=1,2,\dots,q$, w.p. at least $1-q\exp\left(-c\delta_{b}^{2}m\right)$. Hence, with this probability,

[TABLE]

For proving the third claim, recall that $\bm{\hat{x}}_{k}={\bm{U}}\bm{\hat{b}}_{k}$ and $\bm{x}^{*}_{k}={\bm{U}}\bm{g}_{k}+\bm{e}_{k}$. Let $\bm{h}_{k}:=\bm{x}^{*}_{k}-\bm{\hat{x}}_{k}$. We can rewrite $\bm{h}_{k}$ as $\bm{h}_{k}=\bm{x}^{*}_{k}-{\bm{U}}\bm{g}_{k}+{\bm{U}}\bm{g}_{k}-{\bm{U}}\bm{\hat{b}}_{k}$. Thus, by triangle inequality, and using $\|\bm{e}_{k}\|\leq\sin\Theta({\bm{U}},{\bm{U}}^{*})\|\bm{g}_{k}\|\leq\delta_{t}\|\tilde{\bm{b}}^{*}_{k}\|$,

[TABLE]

∎

Proof of Lemma 3.3.

Recall that $\hat{\bm{B}}\overset{\mathrm{QR}}{=}{\bm{R}}_{B}\bm{B}$ and so $\bm{b}_{k}={\bm{R}}_{B}^{-1}\bm{\hat{b}}_{k}$. Using Lemma 3.2, $\|\bm{g}_{k}\|\leq\|\tilde{\bm{b}}^{*}_{k}\|$, and right incoherence (which implies that $\|\tilde{\bm{b}}^{*}_{k}\|^{2}\leq{\sigma_{\max}^{*}}^{2}\mu^{2}r/q$),

[TABLE]

To lower bound $\sigma_{\min}({\bm{R}}_{B})$, observe that $\sigma_{\min}({\bm{R}}_{B})=\sigma_{\min}(\hat{\bm{B}})$. Using Lemma 3.2, the discussion of Sec. A-C, facts from Sec. A-A, and $\sin\Theta({\bm{U}},{\bm{U}}^{*})\leq\delta_{t}$,

[TABLE]

Using $\delta_{t}\leq c/\kappa\sqrt{r}$, $\sigma_{\min}({\bm{R}}_{B})=\sigma_{\min}\left(\hat{\bm{B}}\right)\geq 0.9{\sigma_{\min}^{*}}$. Thus,

[TABLE]

All of the above bounds used the bound from Lemma 3.2. Thus the above bounds hold w.p. at least $1-n^{-10}$ as long as $m\geq C\max(r,\log q,\log n)$. ∎

A-E Proof of the lemmas for Claim 3.4

In this section, we let $\bm{a}_{ik}:=\bm{a}_{ik}^{(T+t)}$ and $\bm{y}_{ik}:=\bm{y}_{ik}^{(T+t)}$. Also, at almost all places, we remove the superscript t.

All the proofs in this section use incoherence of $\bm{B}$ with parameter $\hat{\mu}=C\kappa\mu$ (by Lemma 3.3). This holds w.p. at least $1-n^{-10}$ as long as $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$.

Proof of Lemma 3.10.

Recall that $\mathrm{Term3}(\bm{W}):=\sum_{ik}(\bm{a}_{ik}{}^{\prime}\bm{W}\bm{b}_{k})^{2}$. We have

[TABLE]

Let $X_{ik}=|\bm{a}_{ik}{}^{\prime}\bm{W}\bm{b}_{k}|$. $X_{ik}$ is sub-Gaussian with sub-Gaussian norm $\|\bm{W}\bm{b}_{k}\|$. We use Lemma A.1 with $Y_{ik}=X_{ik}$ and $t=\epsilon_{3}m$, along with the following facts:

$\sum_{k}\|\bm{W}\bm{b}_{k}\|^{4}\leq\max_{k}\|\bm{W}\bm{b}_{k}\|^{2}\sum_{k}\|\bm{W}\bm{b}_{k}\|^{2}$, 2. 2.

$\sum_{k}\|\bm{W}\bm{b}_{k}\|^{2}=\|\bm{W}\bm{B}\|_{F}^{2}=\mathrm{trace}\left(\bm{W}\bm{B}\bm{B}{}^{\prime}\bm{W}^{\prime}\right)=\|\bm{W}\|_{F}^{2}=1$, and 3. 3.

$\max_{k}\|\bm{W}\bm{b}_{k}\|^{2}\leq\|\bm{W}\|^{2}\max_{k}\|\bm{b}_{k}\|^{2}\leq\|\bm{W}\|_{F}^{2}\max_{k}\|\bm{b}_{k}\|^{2}\leq\max_{k}\|\bm{b}_{k}\|^{2}\leq\hat{\mu}^{2}r/q$, w.p. at least $1-n^{-10}$ as long as $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$ by Lemma 3.3.

Using Lemma A.1 and the above facts, for a fixed $\bm{W}$,

[TABLE]

Now we develop an epsilon-net argument to complete the proof. This is inspired by similar arguments in [52]. Recall that $\mathcal{S}_{\bm{W}}=\left\{\bm{W}\in\mathbb{R}^{n\times r},\|\bm{W}\|_{F}=1\right\}$. By [54](Lemma 5.2), there is a set (called epsilon-net), $\bar{\mathcal{S}}_{W}\subseteq\mathcal{S}_{W}$ so that for any $\bm{W}$ in $\mathcal{S}_{\bm{W}}$, there is a $\bar{\bm{W}}\in\bar{\mathcal{S}}_{\bm{W}}$, such that

[TABLE]

and

[TABLE]

Pick $\epsilon_{net}=1/8$ so that $|\mathcal{S}_{W}|\leq 17^{nr}$. Also, define

[TABLE]

so that $\|\mathbf{\Delta}\bm{W}\|_{F}\leq\epsilon_{net}=1/8$.

Using a union bound over all entries in the finite set $\bar{\mathcal{S}}_{W}$,

[TABLE]

Next we extend the above to obtain lower and upper bounds over the entire hyper-sphere, $\mathcal{S}_{W}$. Define

[TABLE]

as the maximum of $\mathrm{Term3}(\bm{W})$ over $\mathcal{S}_{W}$. Since $\frac{\Delta\bm{W}}{\|\Delta\bm{W}\|_{F}}\in\mathcal{S}_{W}$,

[TABLE]

Using this, (14), and Cauchy-Schwarz, w.p. at least $1-2\exp{\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)}$

[TABLE]

The last equality just used $\epsilon_{net}=1/8$ and re-arranged the third term.

If $\theta_{W}/m<1$, we are done because then $\theta_{W}\leq m$. Otherwise, $\theta_{W}/m\geq 1$ and so $\sqrt{\theta_{W}/m}\leq\theta_{W}/m$. Using this and taking $\max_{\bm{W}\in\mathcal{S}_{W}}$ of (15),

[TABLE]

By assumption, $\epsilon_{3}<1/10$, and so the above implies that $\theta_{W}\leq 1.25(1+\epsilon_{3})m$.

Thus, w.p. $1-2\exp{\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)}$,

[TABLE]

We now obtain the lower bound on the minimum of $\mathrm{Term3}$ over the entire hyper-sphere. This uses (14), Cauchy-Schwarz, and the upper bound on $\theta_{W}$ from above. We have

[TABLE]

w.p. $1-2\exp{\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)}$. In the last line we substituted $\epsilon_{net}=1/8$ and used $\epsilon_{3}<1/10$.

All of the above bounds hold on the event in which $\bm{B}$ is $\hat{\mu}$ incoherent. This holds w.p. at least $1-n^{-10}$ as long as $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$ (this follows by Lemma 3.3).

Thus, if $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$, w.p. $1-2\exp{\left(nr(\log 17)-c\frac{\epsilon_{3}^{2}mq}{\hat{\mu}^{2}r}\right)}-n^{-10}$,

[TABLE]

∎

Proof of Lemma 3.11.

Recall that

[TABLE]

where

[TABLE]

First we will show that $\mathbb{E}[\mathrm{Term1}(\bm{W})]=0$.

[TABLE]

where we used $\bm{B}\bm{B}^{\prime}=\bm{I}_{r}$. Next, we will bound $\|\bm{p}_{k}\|$ and $\|\bm{P}\|_{F}$ where $\bm{P}:=[\bm{p}_{1},\bm{p}_{2},\dots,\bm{p}_{q}]$ and use these bounds to show that w.h.p. $\mathrm{Term1}(\bm{W})$ is of order $\epsilon_{1}\delta_{t}\|{\bm{X}^{*}}\|_{F}$ Using $\hat{\bm{X}}={\bm{U}}\hat{\bm{B}}$, $\hat{\bm{B}}\overset{\mathrm{QR}}{=}{\bm{R}}_{B}\bm{B}$ , $\bm{B}\bm{B}^{\prime}=I$,

[TABLE]

Thus,

[TABLE]

Thus, using Lemma 3.2,

[TABLE]

and, writing ${\bm{X}^{*}}-\hat{\bm{X}}=({\bm{U}}{\bm{U}}^{\prime}+(\bm{I}-{\bm{U}}U^{\prime}))({\bm{X}^{*}}-\hat{\bm{X}})$, using Lemma 3.2, $\hat{\bm{X}}={\bm{U}}\hat{\bm{B}}$, and $\bm{G}={\bm{U}}^{\prime}{\bm{X}^{*}}$,

[TABLE]

The last inequality used ${\sigma_{\max}^{*}}\leq\|{\bm{X}^{*}}\|_{F}$. Thus,

[TABLE]

Let

[TABLE]

To bound this, we add and subtract $\hat{\bm{X}}={\bm{U}}\hat{\bm{B}}={\bm{U}}{\bm{R}}_{B}\bm{B}$ from ${\bm{X}^{*}}$ and use the facts that $\bm{B}(\bm{B}^{\prime}\bm{B}-\bm{I})=0$ and $\|\bm{B}\bm{B}^{\prime}-\bm{I}\|\leq 2$ (by triangle inequality and $\|\bm{B}\|=1$). This gives

[TABLE]

We now use above bounds to apply Lemma A.1. Let $X_{ik}=\bm{a}_{ik}^{\prime}\bm{W}\bm{b}_{k}$ and $Y_{ik}=\bm{a}_{ik}^{\prime}\bm{p}_{k}$. Both are sub-Gaussian with $K_{X_{ik}}=\|\bm{W}\bm{b}_{k}\|\leq\|\bm{W}\|_{F}\|\bm{b}_{k}\|\leq\|\bm{b}_{k}\|$, and $K_{Y_{ik}}\leq\|\bm{p}_{k}\|$. Also, using Lemma 3.3, $\bm{b}_{k}$’s are incoherent w.p. at least $1-2q\exp(r\log(17)-\delta_{b}^{2}m)\geq 1-n^{-10}$ if $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$, i.e.,

[TABLE]

Thus, using this and incoherence of $\bm{b}^{*}_{k}$, with above probability,

[TABLE]

Applying Lemma A.1 with $t=m\epsilon_{1}\delta_{t}\|{\bm{X}^{*}}\|_{F}$, and using the above bounds on $\|\bm{p}_{k}\|$, $\|\bm{P}\|_{F}$ and $\max(\|\bm{b}^{*}_{k}\|^{2},\|\bm{b}_{k}\|^{2})$,

[TABLE]

Thus,

[TABLE]

Now we just need to extend our bound for all $\bm{W}\in\mathcal{S}_{W}$. We first extend it to all $\bm{W}$ in an epsilon-net of $\mathcal{S}_{W}$. By [54](Lemma 5.2), there is a net, $\bar{\mathcal{S}}_{W}$ so that for any $\bm{W}$ in $\mathcal{S}_{W}$, there is a $\bar{\bm{W}}$ in $\bar{\mathcal{S}}_{W}$, such that $\|\bar{\bm{W}}-\bm{W}\|_{F}\leq\epsilon_{net}$ and $|\bar{\mathcal{S}}_{W}|\leq\left(1+\frac{2}{\epsilon_{net}}\right)^{nr}$. Pick $\epsilon_{net}=1/8$. With this, $|\bar{\mathcal{S}}_{W}|\leq 17^{nr}$. Define $\mathbf{\Delta}\bm{W}:=\bar{\bm{W}}-\bm{W}$. We have $\|\mathbf{\Delta}\bm{W}\|_{F}\leq\epsilon_{net}=1/8$. Using union bound on the set $\bar{\mathcal{S}}_{W}$,

[TABLE]

To extend the claim to all $\bm{W}\in\mathcal{S}_{W}$, define

[TABLE]

Since $\frac{\Delta\bm{W}}{\|\Delta\bm{W}\|_{F}}\in\mathcal{S}_{W}$, $\sum_{ik}(\bm{a}_{ik}{}^{\prime}\Delta{\bm{W}}\bm{b}_{k})(\bm{a}_{ik}^{\prime}\bm{p}_{k})\leq\theta_{W}\|\Delta\bm{W}\|_{F}\leq\theta_{W}\epsilon_{net}$. Thus, using (18), for any $\bm{W}\in\mathcal{S}_{W}$,

[TABLE]

w.p. at least $1-2\exp\left(nr(\log 17)-c\frac{mq\epsilon_{1}^{2}}{\kappa^{2}\mu^{2}r}\right)$. Thus, taking the $\max_{\bm{W}\in\mathcal{S}_{W}}$ of the above equation and solving for $\theta_{W}$, with the above probability,

[TABLE]

All of the above bounds hold on the event in which $\bm{B}$ is $\hat{\mu}$ incoherent. Thus, if $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$, w.p. at least $1-2\exp\left(nr(\log 17)-c\frac{mq\epsilon_{1}^{2}}{\kappa^{2}\mu^{2}r}\right)-n^{-10}$, $\max_{\bm{W}\in\mathcal{S}_{W}}|\mathrm{Term1}(\bm{W})|\leq(8/7)m\epsilon_{1}^{2}\|{\bm{X}^{*}}\|_{F}$.

Finally, (17) also implies that

[TABLE]

We use the above bound in a later proof. ∎

Proof of Lemma 3.12.

Recall that $\bm{a}_{ik}:=\bm{a}_{ik}^{(T+t)}$ and same for $\bm{y}_{ik}$. Thus, these are independent of the current $\bm{\hat{x}}_{k}$’s.

By Cauchy-Schwarz,

[TABLE]

We can bound the first term using Lemma 3.10. Consider the second term. Since $\bm{c}_{ik}=sign(\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*})$ and $\hat{\bm{c}}_{ik}=sign(\bm{a}_{ik}{}^{\prime}\bm{\hat{x}}_{k})$, clearly $(\bm{c}_{ik}\hat{\bm{c}}_{ik}-1)^{2}=(4)\mathbbm{1}_{\left\{\bm{c}_{ik}\neq\hat{\bm{c}}_{ik}\right\}}$. To bound this term we use the following result.

Lemma A.2** (Lemma 1 of [10]).**

Let $\bm{a}_{i}$ be standard Gaussian random vectors. For any given $\bm{x}^{*}$, and $\bm{\hat{x}}$ independent from $\bm{a}_{i},i=1,\cdots,m$,

[TABLE]

for all $\bm{\hat{x}}$ that satisfy $\mathrm{dist}(\bm{\hat{x}},\bm{x}^{*})\leq 0.4$. Here $\mbox{erfc}(u):=\frac{2}{\sqrt{\pi}}\int_{u}^{\infty}\exp\left(-\tau^{2}\right)d\tau$ is the complementary error function.

Let $Q_{ik}:=\mathbbm{1}_{\left\{\bm{c}_{ik}\neq\hat{\bm{c}}_{ik}\right\}}\cdot\left(\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}\right)^{2}$ and let $Z_{ik}:=\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}$. Recall from Sec. A-C that $\mathrm{dist}(\bm{x}_{k},\bm{x}^{*}_{k})=\|\bm{h}_{k}\|$ with $\bm{h}_{k}=\bm{x}_{k}^{*}-\bm{\hat{x}}_{k}$.

We first upper bound $\mathbb{E}[Q_{ik}|\bm{\hat{x}}_{k}]$. For simplicity, we remove the subscripts $i$ and $k$ wherever these are not needed. Consider $\mathbb{E}[Q]=\mathbb{E}[\mathbbm{1}_{\bm{c}\neq\hat{\bm{c}}}\ Z^{2}]$. Observe that $\bm{c}=sign(Z)$ is a function of $Z$ and $Z$ depends on $\bm{a}$. Also $\hat{\bm{c}}$ is a function of $\bm{a}$. Thus both of $\bm{c},\hat{\bm{c}}$ are dependent on $Z$. Moreover $\hat{\bm{c}}$ also depends on $\bm{\hat{x}}$. We first bound $\mathbb{E}[Q|Z^{2},\bm{\hat{x}}]$ using the lemma stated above. For any $\bm{\hat{x}}$ that satisfies $\mathrm{dist}(\bm{\hat{x}},\bm{x}^{*})\leq 0.4$,

[TABLE]

The first inequality follows using Lemma A.2, the second is a standard upper bound on the erfc function [69]. Thus, for any $\bm{\hat{x}}$ that satisfies $\mathrm{dist}(\bm{\hat{x}},\bm{x}^{*})\leq 0.4$,

[TABLE]

Since $Z$ is zero mean Gaussian with variance $\|\bm{x}^{*}\|^{2}$, $Y:=Z^{2}/\|\bm{x}^{*}\|^{2}$ is standard chi-squared with one degree of freedom. Thus, $\mathbb{E}[Y]=1$. Using this and $\exp(-y/2)<1$, we get

[TABLE]

where $\Gamma(1/2)$ is the Gamma function evaluated at $1/2$ (can treat it as a constant). We will now use Lemma 3.2 to average over $\bm{\hat{x}}$. Let $E$ be the event that $\mathrm{dist}(\bm{x}^{*}_{k},\bm{\hat{x}}_{k})\leq\delta_{t}\|\bm{x}^{*}_{k}\|$ for all $k=1,2,\dots,q$. By Lemma 3.2, under the lower bound on $m$, this event occurs w.p. at least $1-n^{-10}$. On the complement event, we do not have any tight bounds on $Q$. However, trivially, $Q\leq\left(\bm{a}_{i}{}^{\prime}\bm{x}^{*}\right)^{2}$ always and so, for any $\bm{\hat{x}}$, $\mathbb{E}[Q|\bm{\hat{x}}]\leq\|\bm{x}^{*}\|^{2}$. Thus

[TABLE]

We will eventually set $\delta_{t}=c/(\kappa^{2}r)$ which is larger than $n^{-10}$. Thus,

[TABLE]

and so, if $m\geq C\max(r,\log n,\log q)$,

[TABLE]

As a side-note, we should clarify here that the dependence on $\bm{\hat{x}}$ matters in only the above expected value computation because this is the only term where we have upper bounded the expectation using Lemma 3.2. The expected value of the other two terms is the same for all values of $\bm{\hat{x}}$ and hence we ignore the dependence there. Everywhere else, we use this lemma only while obtaining the high probability error bounds, and of course we assume the bounds hold on the intersection of the desired event with $E$.

Next we show that, on the event $E$, whp, $\sum_{ik}Q_{ik}$ is of the same order. As shown in the proof of Theorem 1 of [10], $\bm{c}_{ik}\neq\hat{\bm{c}}_{ik}$ implies that888This follows from $(\bm{a}^{\prime}\bm{x}^{*})^{2}=(\bm{a}^{\prime}\bm{h})^{2}+(\bm{a}^{\prime}\bm{\hat{x}})^{2}+2(a^{\prime}\bm{h})(a^{\prime}\bm{\hat{x}})=(\bm{a}^{\prime}\bm{h})^{2}-(\bm{a}^{\prime}\bm{\hat{x}})^{2}+2(\bm{a}^{\prime}\bm{x}^{*})(\bm{a}^{\prime}\bm{\hat{x}})\leq(\bm{a}^{\prime}\bm{h})^{2}$. The inequality holds because $\bm{c}_{ik}\neq\hat{\bm{c}}_{ik}$ means that the last term is negative. $\left(\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}\right)^{2}\leq\left(\bm{a}_{ik}{}^{\prime}\bm{h}_{k}\right)^{2}$. Here $\bm{h}_{k}=\bm{x}_{k}^{*}-\bm{\hat{x}}_{k}$. Thus, $Q_{ik}=\mathbbm{1}_{\left\{\bm{c}_{ik}\neq\hat{\bm{c}}_{ik}\right\}}\left(\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}\right)^{2}\leq\mathbbm{1}_{\left\{\bm{c}_{ik}\neq\hat{\bm{c}}_{ik}\right\}}\left(\bm{a}_{ik}{}^{\prime}\bm{h}_{k}\right)^{2}\leq\left(\bm{a}_{ik}{}^{\prime}\bm{h}_{k}\right)^{2}$. Thus, it is a sub-exponential r.v., or equivalently it is a product of sub-Gaussian r.v.’s $\sqrt{Q_{ik}}$. Thus, we can apply Lemma A.1 with $K_{X_{ik}}=K_{Y_{ik}}=\|\bm{h}_{k}\|$. Pick $t=m\epsilon_{2}\delta_{t}^{2}\|{\bm{X}^{*}}\|_{F}^{2}$. Observe that

[TABLE]

The second inequality used $\|\bm{h}_{k}\|^{2}\leq\delta_{t}^{2}\|\bm{x}^{*}_{k}\|^{2}$ which holds on the event $E$. The third inequality used the fact that right incoherence implies $\|\bm{x}^{*}_{k}\|^{2}\leq\mu^{2}\kappa^{2}\|{\bm{X}^{*}}\|_{F}^{2}/q$. Thus, on the event $E$,

[TABLE]

Thus, using (20), and $\Pr(E)>1-n^{-10}$ (by Lemma 3.2), if $m\geq C\max(r,\log n,\log q)$, w.p. at least $1-n^{-10}-2\exp{\left(-c\epsilon_{2}^{2}mq/\kappa^{2}\mu^{2}\right)}$,

[TABLE]

Finally, combining (19), (21), and Lemma 3.10, if $m\geq C\max(r,\log n,\log q)/\delta_{b}^{2}$, w.p. at least $1-n^{-10}-2\exp\left(nr-c\epsilon_{3}^{2}\frac{mq}{\hat{\mu}^{2}r}\right)-2\exp{\left(-c\epsilon_{2}^{2}\frac{mq}{\kappa^{2}\mu^{2}}\right)}$,

[TABLE]

∎

Proof of Lemma 3.13.

Using facts from Sec. A-A, and using $\sigma_{i}({\bm{U}}^{*})=1$,

[TABLE]

To upper bound $\sin\Theta({\bm{B}^{*}}^{\prime},\bm{B}^{\prime})$, first notice that ${\bm{B}^{*}}^{\prime}$ and $\bm{B}^{\prime}$ are basis matrices. Thus, $\sin\Theta({\bm{B}^{*}}^{\prime},\bm{B}^{\prime})=\|\bm{B}^{*}\left(\bm{I}-{\bm{B}}^{\prime}\bm{B}\right)\|$. We have upper bounded $\|\tilde{\bm{B}}^{*}\left(\bm{I}-{\bm{B}}^{\prime}\bm{B}\right)\|$ in Lemma 3.11. Also recall that $\tilde{\bm{B}}^{*}=\bm{\Sigma}^{*}{\bm{B}^{*}}$ and $\|\tilde{\bm{B}}^{*}\|_{F}=\|{\bm{X}^{*}}\|_{F}\leq\sqrt{r}{\sigma_{\max}^{*}}$. Thus,

[TABLE]

Thus, using $\delta_{t}\leq 0.7/C\sqrt{r}\kappa$, we will have $\sigma_{\min}({\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}\bm{B}^{\prime})\geq 0.95{\sigma_{\min}^{*}}$. ∎

Appendix B Proof of Lemma 3.9

We begin by defining a few quantities needed for writing an expression of $\hat{\bm{U}}^{t+1}$ in closed form.

We use the subscript $\bm{M}_{vec}$ to refer to the vectorized version of matrix $\bm{M}$. Thus, for example, ${\bm{U}}^{*}_{vec}$ is an $nr\times 1$ vectorized version of the basis matrix ${\bm{U}}^{*}$. When updating its estimate by LS, we first obtain an expression for the vectorized version $\hat{\bm{U}}^{t+1}_{vec}$ and then rearrange it as an $n\times r$ matrix. 2. 2.

Scalars such as $b(q)$ refer to the $q$-th entry of vector $\bm{b}$. 3. 3.

Define the diagonal matrix ${\bm{C}}_{k}:=\text{diag}(\bm{c}_{ik})$, and recall from the algorithm that $\bm{\hat{C}}_{k}=\text{diag}(\hat{\bm{c}}_{ik})$. Here $\bm{c}_{ik}:=\text{sign}({\bm{a}_{ik}{}^{\prime}\bm{x}_{k}^{*}})$ and $\hat{\bm{c}}_{ik}:=\text{sign}({\bm{a}_{ik}{}^{\prime}\bm{\hat{x}}_{k}})$ as defined in Lemma 3.9. 4. 4.

For $k=1,\ldots,q$, define

[TABLE]

Observe that $\bm{B}_{k,mat}$ and ${\bm{M}}$ are $nr\times n$ and $nr\times nr$ matrices respectively while $\bm{d}$ is an $nr\times 1$ vector. Moreover, for any $\bm{w}\in\mathbb{R}^{nr\times 1}$ it is easy to see that $\bm{w}^{\prime}\bm{B}_{k,mat}=\bm{b}_{k}^{\prime}\bm{W}^{\prime}$ where $\bm{W}$ is matrix version of $\bm{w}$ with $\bm{W}\in\mathbb{R}^{n\times r}$.

Recall from the algorithm that

[TABLE]

This is an LS problem, it can be solved in closed form to give

[TABLE]

which is $nr$-length vector. We get the matrix $\hat{\bm{U}}^{t+1}$ by reshaping this vector into an $n\times r$ matrix.

To simplify the above expression, first recall that $\bm{y}_{k}=|\bm{A}_{k}{}^{\prime}\bm{x}^{*}_{k}|$ and $\bm{C}_{k}$ is the diagonal matrix containing the signs of $(\bm{A}_{k}{}^{\prime}\bm{x}^{*}_{k})_{i}$: its $i,i$-th entry is the sign of $(\bm{a}_{ik}{}^{\prime}\bm{x}^{*}_{k})$. Thus,

[TABLE]

Since $\bm{x}^{*}_{k}={\bm{U}}^{*}\bm{\Sigma}^{*}\bm{b}_{k}^{*}$, we can rewrite $\bm{d}$ as

[TABLE]

Before proceeding further, we define a few more quantities.

Define

[TABLE] 2. 2.

Define the “expanded” singular value matrix which is of size $nr\times nr$,

[TABLE]

where $\sigma_{i}^{*}$ are the singular values of ${\bm{X}^{*}}$. 3. 3.

In order to separate the contribution of phase error from the rest, split $\bm{d}$ as $\bm{d}=\bm{d}^{(1)}+\bm{d}^{(2)}$ where

[TABLE]

Thus,

[TABLE] 4. 4.

Define the $nr$-length vector ${\bm{F}}_{vec}$ as follows

[TABLE]

and let ${\bm{F}}\in\mathbb{R}^{n\times r}$ be the reshaped matrix formed from ${\bm{F}}_{vec}$.

We will now show that

[TABLE]

This will be useful because when we try to bound $\sin\Theta({\bm{U}}^{*},{\bm{U}}^{t+1})$, the first term will disappear. To do this, we add and subtract the vector ${\bm{M}}{\bm{S}}\bm{\Sigma}^{*}_{big}{\bm{U}}^{*}_{vec}$ from $\bm{d}^{(1)}$. This gives

[TABLE]

Next we explain why the $n\times r$ reshaped matrix version of the vector ${\bm{S}}\bm{\Sigma}_{big}^{*}{\bm{U}}^{*}_{vec}$ equals ${\bm{U}}^{*}\bm{\Sigma}^{*}{\bm{B}^{*}}\bm{B}^{\prime}$. We have

[TABLE]

Matrix version of the above vector has $p$-th column as ${\bm{U}}^{*}\bm{\Sigma}^{*}\left(\sum_{k}\bm{b}_{k}^{*}b_{k}(p)\right)$. This implies that matrix version of this vector is ${\bm{U}}^{*}\bm{\Sigma}^{*}{\bm{B}^{*}}\bm{B}^{\prime}$. Moreover, it is easy to see that $\bm{B}_{k,mat}^{\prime}{\bm{S}}\bm{\Sigma}_{big}^{*}{\bm{U}}^{*}_{vec}={\bm{U}}^{*}\bm{\Sigma}^{*}{\bm{B}^{*}}\bm{B}^{\prime}\bm{b}_{k}$.

In the rest of this section, we use (22) to obtain the desired bound on $\sin\Theta({\bm{U}}^{*},{\bm{U}}^{t+1})$. Recall that $\hat{\bm{U}}_{t+1}\overset{\mathrm{QR}}{=}{\bm{U}}^{t+1}{\bm{R}}_{U}$. Thus, ${\bm{U}}^{t+1}=\hat{\bm{U}}^{t+1}({\bm{R}}_{U})^{-1}=({\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}{\bm{B}}^{\prime}-{\bm{F}})({\bm{R}}_{U})^{-1}$ and so

[TABLE]

Since $\sigma_{\min}({\bm{R}}_{U})=\sigma_{\min}({\bm{U}}^{t+1})$, we have $\sigma_{\min}({\bm{R}}_{U})=\sigma_{\min}({\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}\bm{B}^{\prime}-{\bm{F}})\geq\sigma_{\min}({\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}\bm{B}^{\prime})-\|{\bm{F}}\|\geq\sigma_{\min}({\bm{U}}^{*}\bm{\Sigma}^{*}\bm{B}^{*}\bm{B}^{\prime})-\|{\bm{F}}_{vec}\|$. Thus,

[TABLE]

In the rest of this proof, we show that $\|{\bm{F}}_{vec}\|$ is upper bounded by $\mathrm{MainTerm}$. We have

[TABLE]

Consider the first term, ${\bm{M}}^{-1}$. Since ${\bm{M}}$ is a symmetric positive semidefinite matrix

[TABLE]

For all $\bm{w}\in\mathbb{{\bm{R}}}^{nr\times 1},\|\bm{w}\|^{2}=1$, we can write $\bm{w}^{\prime}{\bm{M}}\bm{w}$

[TABLE]

where $\bm{W}\in\mathbb{{\bm{R}}}^{n\times r}$ is the matrix version of $\bm{w}$ and $\bm{w}=\bm{W}_{vec}$. Recall that $\mathcal{S}_{W}=\{\bm{W}\in\mathbb{R}^{n\times r},\ \|\bm{W}\|_{F}=1\}=\{\bm{w}\in\mathbb{R}^{nr\times 1},\ \|\bm{w}\|=1\}$. Thus,

[TABLE]

Now consider the first term inside the parenthesis. Using the variational definition,

[TABLE]

It follows from definitions that

[TABLE]

Similarly

[TABLE]

and thus

[TABLE]

For the final term, $\|\bm{d}^{(2)}\|$, by variational definition we have,

[TABLE]

From definitions we know that

[TABLE]

and thus

[TABLE]

Combining (25) - (28),

[TABLE]

Combining the above bound with (24) proves the lemma.

Appendix C Proof Sketch of Corollary 5.1

Suppose first that the subspace change times $k_{j}$ were known. By our assumption, $k_{j+1}-k_{j}>T\alpha$. Then the proof is almost exactly the same as that for the static case. The only difference is that, in the current case, every $\alpha$ time instants, we are using measurements corresponding to a new set of $\alpha$ signals (columns of ${\bm{X}^{*}}_{\mathrm{full}}$) but we use the estimate of the subspace obtained from the measurements for the previous $\alpha$ time instants. As long as the subspace has not changed between the two intervals, Claims 3.1 and 3.4 apply without change. Combining them, we can again conclude that $\sin\Theta({\bm{U}}_{{\mathrm{sub}},(j)}^{0},{\bm{U}}^{*}_{{\mathrm{sub}},(j)})\leq\delta_{\mathrm{init}}$ at $k=k_{j}+\alpha$, and that the bound decreases $0.7$ times after each $\alpha$-length epoch so that $\sin\Theta({\bm{U}}_{{\mathrm{sub}},(j)}^{T},{\bm{U}}^{*}_{{\mathrm{sub}},(j)})\leq\epsilon$ at $k=k_{j}+\alpha T$. By our assumption, $k_{j+1}>k_{j}+\alpha T$ so this happens before the next change.

The proof in the unknown $k_{j}$ case follows if we can show that, whp, $k_{j}\leq\hat{k}_{j}\leq\hat{k}_{j}+2\alpha$. This can be done using Lemma C.1 given below along with the following argument borrowed from [55, 70]. Consider the $\alpha$-length interval in which $k_{j}$ lies. Assume that, before this interval, we have an $\epsilon$-accurate estimate of the previous subspace. In this interval, the first some data vectors satisfy $\bm{x}^{*}_{k}={\bm{U}}^{*}_{{\mathrm{sub}},(j-1)}\tilde{\bm{d}}^{*}_{k}$, while the rest satisfy $\bm{x}^{*}_{k}={\bm{U}}^{*}_{{\mathrm{sub}},(j)}\tilde{\bm{d}}^{*}_{k}$. By our assumption, this interval lies in the “detect phase”. We cannot guarantee whether the change will get detected in this interval, but it may. However, in the interval after this interval, all data vectors satisfy $\bm{x}^{*}_{k}={\bm{U}}^{*}_{{\mathrm{sub}},(j)}\tilde{\bm{d}}^{*}_{k}$. In this interval, Lemma C.1 given below can be used to show that the change gets detected whp. Thus, either the change is detected in the first interval itself (the one that contains $k_{j}$), or it is not. If it is not, then, by Lemma C.1, whp, it will get detected in the second interval (in which all signals are generated from the $j$-th subspace). Thus, $\hat{k}_{j}\leq k_{j}+2\alpha$. See Appendix A of [55] for a precise proof of this idea. The key point to note here is that we are never updating the subspace in the interval that contains $k_{j}$ and hence we do not have to prove a new descent lemma that deals with the interval in which the subspace changes.

We will replace $\alpha$ by $q$ in the following lemma and its proof, in order to able to use bounds from earlier proofs. Thus in this lemma, we are considering a $q$-frame epoch.

Lemma C.1**.**

Consider the $(n-r)\times(n-r)$ matrix

[TABLE]

Assume that $\|\bm{Y}_{U}-\mathbb{E}[\bm{Y}_{-}]\|\leq\frac{\delta_{\mathrm{init}}{\sigma_{\min}^{*}}^{2}}{q}$. This is true by (10). Assume that $\sin\Theta({\bm{U}}_{{\mathrm{sub}},(j-1)},{\bm{U}}^{*}_{{\mathrm{sub}},(j-1)})\leq\epsilon$. Then,

If $\mathcal{J}_{q}\subseteq[k_{j},k_{j+1})$ (change has occurred), then

[TABLE] 2. 2.

If $\mathcal{J}_{q}\subseteq[\hat{k}_{j-1}+Tq,k_{j})$ (change has not occurred), then

[TABLE]

Proof.

This proof uses the following fact: For basis matrices ${\bm{P}}_{1},{\bm{P}}_{2},{\bm{P}}_{3}$ of the same size, $\sin\Theta\left({\bm{P}}_{1},{\bm{P}}_{2}\right)-2\sin\Theta\left({\bm{P}}_{2},{\bm{P}}_{3}\right)\leq\sin\Theta\left({\bm{P}}_{1},{\bm{P}}_{3}\right)\leq\sin\Theta\left({\bm{P}}_{1},{\bm{P}}_{2}\right)+\sin\Theta\left({\bm{P}}_{2},{\bm{P}}_{3}\right)$.

Define the $(n-r)\times(n-r)$ matrix

[TABLE]

Proof of item 1

[TABLE]

Also we have

[TABLE]

Thus using the facts from Sec. A-A and $\min_{k}\beta_{1,k}^{-}\geq 1.5$ (proved while proving Claim 3.1 for initializing ${\bm{U}}^{*}$),

[TABLE]

Proof of item 2

[TABLE]

It is easy to see that

[TABLE]

∎

Author Biographies