Low Rank Phase Retrieval

I Introduction

In recent years there has been a large amount of work on the phase retrieval (PR) problem and on its generalization. The original PR problem involves recovering a length-$n$ signal $\bm{x}$ from the magnitudes of its discrete Fourier transform (DFT) coefficients. Generalized PR replaces the DFT by inner products with any set of measurement vectors, $\bm{a}_{i}$. Thus, the goal is to recover $\bm{x}$ from $|\bm{a}_{i}{}^{\prime}\bm{x}|^{2}$, $i=1,2,\dots,m$. These magnitude-only measurements are referred to as phaseless measurements. PR is a classical problem that occurs in many applications such as X-ray crystallography, astronomy, and ptychography because the phase information is either difficult or impossible to obtain . Algorithms for solving it have existed since the work of Gerchberg and Saxton and Fineup . In recent years, there has been much renewed interest in PR, e.g., and in sparse PR, e.g., .

In more recent works, non-convex methods, that do not lift the problem to higher dimensions, have been explored along with provable guarantees . An alternating minimization (AltMin) technique with spectral initialization, AltMinPhase, was developed and analyzed in . The AltMin step of this approach is essentially the same as the old Gerchberg-Saxton algorithm . A gradient descent method with spectral initialization, called Wirtinger Flow (WF), was studied in . In , truncated WF (TWF), which introduced a truncation technique to further improve WF performance, was developed. It was shown that TWF recovers $\bm{x}$ from only $cn$ iid Gaussian phaseless measurements, while the number of iterations needed for getting an error of order $\epsilon$ is $c\log(1/\epsilon)$ (converges geometrically). AltMinPhase and WF require more measurements, $cn\log^{3}n$ and $cn\log n$, respectively. WF also has a slower convergence rate. Two recent modifications of TWF have the same order complexities but improved empirical performance.

Problem Setting. In this work, instead of a single vector $\bm{x}$, we consider a set of $q$ vectors, $\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{q}$, such that the $n\times q$ matrix,

has rank $r\ll\min(n,q)$. For each column $\bm{x}_{k}$ of $\bm{X}$, we observe a set of $m$ measurements of the form

The measurement vectors, $\bm{a}_{i,k}$, are mutually independent. Our goal is to recover the matrix $\bm{X}$ from these $mq$ phaseless measurements $\bm{y}_{i,k}$. Since we have magnitude-only measurements of each column $\bm{x}_{k}$, we can only hope to recover each column $\bm{x}_{k}$ up to a global phase ambiguity. We refer to the above problem as low rank phase retrieval (LRPR).

A motivating application for LRPR is dynamic astronomical imaging such as solar imaging where the sun’s surface properties gradually change over time . The changes are usually due to a much smaller number of factors, $r$, than the size of the image, $n$, or the total number of images, $q$. If the images are arranged as 1D vectors $\bm{x}_{k}$, then the resulting matrix is approximately low rank. As another potential application, consider a Fourier ptychography imaging system that captures a dynamic scene exhibiting a temporal evolution; this is often the case when observing live biological specimens in vitro. 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, then the matrix formed by stacking the frames next to each other can be modeled as a rank-$r$ matrix, where $r\ll\min(q,n)$. Similar applications involving a sequence of gradually changing images also occur in X-ray and sub-diffraction imaging systems. Moreover, if we are only interested in identifying the principal directions of variation of the image sequences, then the problem becomes that of phaseless PCA. PCA is often the first step for classification, clustering, modeling, or other exploratory data analysis.

Contributions. This work has two contributions. We propose iterative algorithms for solving the LRPR problem described above. Our solution approach relies on the fact that a rank $r$ matrix $\bm{X}$ can be expressed (non-uniquely) as $\bm{X}={\bm{U}}\bm{B}$ where ${\bm{U}}$ is an $n\times r$ matrix with mutually orthonormal columns. Its first step consists of a spectral initialization step, motivated by TWF, for first initializing ${\bm{U}}$, and then, the columns of $\bm{B}$. The remainder of the algorithm is developed in one of two ways: using a projected gradient descent strategy to modify the TWF iterates (LRPR1); or an AltMin algorithm, motivated by AltMinPhase, that directly exploits the decomposition $\bm{X}={\bm{U}}\bm{B}$ (LRPR2). Via extensive experiments, we demonstrate that both LRPR1 and LRPR2 have better sample complexity than TWF; with LRPR2 being the best. Moreover, when enough measurements are available for TWF to work, we show that the LRPR initialization can also be used to speed up basic TWF for solving LRPR.

Our second, and most important, contribution is a sample complexity bound for the proposed initialization to get within an $\varepsilon$ ball of the true $\bm{X}$. Our results show that, if the goal is to only initialize ${\bm{U}}$ with subspace recovery error below a fixed level, say $\varepsilon=1/4$, then a total of $mq=cnr^{2}/\varepsilon^{2}=16cnr^{2}$ iid Gaussian measurements suffice with high probability (whp). When $r$ is small, $nr^{2}$ is only slightly larger than $nr$ which is the minimum required by any technique to recover the span of ${\bm{U}}$. If the goal is to also initialize the $\bm{x}_{k}$’s with normalized error below say $\varepsilon=1/4$, then we need more measurements, but still significantly fewer than TWF. For example, if $r\leq c\log n$ and $q\geq cn$, then, only $16c\sqrt{n}$ measurements per column are required. We note that our guarantees assume that a different set of measurements is used for initializing ${\bm{U}}$ and $\bm{B}$ (see Model 3.1).

As seen in many earlier works, e.g., AltMinPhase , resampled WF [10, Algorithm 2 and Theorem 5.1] or TWF , the sample complexity of the entire algorithm is equal to or smaller than that of the initialization step for a fixed error levelFor AltMinPhase, the initialization sample complexity (for achieving a given fixed error) is $cn\log^{3}n$ while it is only $cn\log n$ per iteration for the rest of the algorithm. For resampled WF, it is $cn\log^{2}n$ for initialization and $cn\log n$ for the rest of the algorithm, while for TWF, it is $cn$ both for the initialization and for the complete algorithm.. This is why initialization guarantees are important.

Our problem setting assumes a different (mutually independent) set of measurement vectors is used for imaging each column $\bm{x}_{k}$. This is critical for guaranteeing the improved sample complexity of our solution approach over single-vector PR methods because this is what ensures that the $mq$ matrices $\bm{y}_{i,k}\bm{a}_{i,k}\bm{a}_{i,k}{}^{\prime}$ are all mutually independent conditioned on $\bm{X}$. Hence, we can exploit averaging over $mq$ such matrices when estimating ${\bm{U}}$. If $\bm{a}_{i,k}=\bm{a}_{i,1}$ (same $\bm{a}_{i}$’s are used), then this benefit disappears since only $m$ of the above matrices are mutually independent. We demonstrate this in Table I (last column) in Sec. V. We discuss the practical implications of our setting in Sec. III-D.

Paper Organization. In Sec. II, we develop the proposed LRPR initialization approach (LRPR-init). We obtain sample complexity bounds for it in Sec. III. In Sec. IV, we explain how LRPR-init can be used to develop iterative algorithms for LRPR that are either faster than basic TWF (LRPR+TWF) or need a smaller $m$ to work (LRPR1 and LRPR2). Numerical experiments backing our claims are shown in Sec. V. We prove our results from Sec. III in Sec. VI, and conclude in Sec. VII.

The algorithms proposed in this work are applicable for both real and complex measurements. Experiments are shown for both cases too. Moreover, as shown in our experiments, our algorithms also apply to noisy measurements. However, for simplicity, we state and prove our guarantees only for the real Gaussian measurements’ case. Their extension to complex Gaussian measurements is straightforward.

II Low Rank PR (LRPR) Initialization

LRPR-init is a two step approach. We first initialize ${\bm{U}}$ using a truncated spectral initialization idea . For this, define

However, as explained in , because $\bm{y}_{i,k}\bm{a}_{i,k}\bm{a}_{i,k}{}^{\prime}$ can be written as $\bm{w}\bm{w}^{\prime}$ with $\bm{w}$ a heavy-tailed random vector, more samples will be needed for the law of large numbers to take effect than if $\bm{w}$ were not heavy-tailed. To remedy this situation, we use the truncation idea suggested in and compute ${\bm{\hat{U}}}$ as the top $r$ eigenvectors of

The idea of truncation is to average only over those $(i,k)$’s for which $\bm{y}_{i,k}$ is not too far from its empirical mean.

Next we consider initialization of the $\bm{b}_{k}$’s. Define the matrix

Suppose that ${\bm{\hat{U}}}$ is independent of the $\bm{M}_{k}$’s. Then, from (2), conditioned on ${\bm{\hat{U}}}$,

The complete approach, LRPR-init, is summarized in Algorithm 1. Note that this uses the same set of measurements to recover ${\bm{U}}$ and $\bm{b}_{k}$’s. But, as seen from our numerical experiments, it still works well in practice. For our analysis in Sec. III, we assume that a new set of measurements is available for computing $\bm{\hat{g}}_{k}$, and thus $\bm{M}_{k}$ is independent of ${\bm{\hat{U}}}$.

Algorithm 1 also estimates the rank $r$ automatically by looking for the maximum gap between consecutive eigenvalues of $\bm{Y}_{U}$. As we explain in Sec. III-C, under a simple assumption on the eigenvalues of $\bar{\bm{\Lambda}}$, this returns the correct rank whp.

II-B Projected-TWF initialization

Another way to obtain an initial estimate of the low rank matrix $\bm{X}$ would be to project the matrix formed by the TWF initialization for each column $\bm{x}_{k}$ onto the space of rank $r$ matrices. This is summarized in Algorithm 2. However, as we show in Sec. V, Tables I and II, this approach performs much worse than LRPR-init. The reason is that it does not simultaneously exploit averaging of the matrices $\bm{y}_{i,k}\bm{a}_{i,k}\bm{a}_{i,k}{}^{\prime}$ over both $i$ and $k$.

III Sample Complexity Bounds for LRPR-init

In this section, we obtain sample complexity bounds for getting a provably accurate initial estimate of both ${\bm{U}}$ and of the $\bm{x}_{k}$’s whp. For simplicity, our results assume iid real Gaussian measurement vectors, $\bm{a}_{i,k}$. As will be evident from the proofs, the extension to complex Gaussian vectors is straightforward. In Sec. III-A, we provide a guarantee for the case when $\bm{X}$ is a deterministic unknown matrix with known rank $r$. These hold whp over measurement vectors $\bm{a}_{i,k}$. In Sec. III-B, we give results for the case of $\bm{X}$ being random with known rank $r$. These hold whp both over matrices $\bm{X}$ generated from the assumed probability distribution and over measurement vectors $\bm{a}_{i,k}$. In Sec. III-C, we show how we can extend both sets of results to the unknown rank case.

With measurements taken as above, we study Algorithm 3.

and $\bm{X}={\bm{U}}\bm{B}$. Thus, $\bar{\bm{\Lambda}}=\frac{1}{q}\sum_{k}\bm{b}_{k}\bm{b}_{k}{}^{\prime}$. Let $\bar{\lambda}_{\max}$ and $\bar{\lambda}_{\min}$ denote the maximum and minimum eigenvalues of $\bar{\bm{\Lambda}}$. Define

Consider an unknown deterministic rank $r$ matrix $\bm{X}$. Assume that the measurements of its columns are generated according to Model 3.1. Consider the output of Algorithm 3 (known $r$ case). Suppose that $r\leq cn^{1/5}$. For an $\varepsilon<1$, if

then, with probability at least $1-4\exp(-cn)-\frac{32q}{n^{4}}$,

Furthermore, if $q\leq cn^{2}$, then the above event holds with probability at least $1-c/n^{2}$.

Notice that our lower bounds depend on $\kappa^{2}$ where $\kappa$ is the condition number of $\bm{X}\bm{X}^{\prime}$. This is pretty typical, e.g., it is also the case in and many other works. It may be possible to remove this dependence by borrowing ideas from . A second point to note is that the probability of the good event depends inversely on $q$. This dependence comes from needing to ensure that each of the $q$ vectors $\bm{x}_{k}$ are accurately recovered. However, the dependence is pretty weak: when $q<cn^{2}$, the probability can be further lower bounded by $1-c/n^{2}$.

When the goal is to only recover ${\bm{U}}$ with subspace error at most $\varepsilon$ (and not the $\bm{x}_{k}$’s), the required lower bounds can be relaxed further. In particular, we have the following corollary.

Recall that ${\bm{U}}$ is an $n\times r$ matrix and hence has $nr$ unknowns. From Corollary 3.3, for a fixed $\varepsilon$, $\rho$, and $\kappa$, one needs a total of only $mq=cnr^{2}$ measurements to recover ${\bm{U}}$. When $r$ is small, e.g., $r=c\log n$, this is only slightly more than the minimum required which would be $nr$.

III-B Main Results for Random 𝐗𝐗\bm{X} - Known rank case

First consider an independent zero mean Gaussian model on the $\bm{b}_{k}$’s.

let $\bar{\lambda}_{\min}$ be its minimum eigenvalue, $\bar{\lambda}_{\max}$ its maximum eigenvalue, and $\kappa:=\frac{\bar{\lambda}_{\max}}{\bar{\lambda}_{\min}}$ its condition number. Assume also that, for all $k=1,2,\dots,q$,

This is ensured, for example, if $\max_{k}\lambda_{k,\max}\leq c\min_{k}\lambda_{k,\max}$.

In using Model 3.4, there are two main changes. The first is that we need to apply a law of large numbers result to show that $\frac{1}{q}\sum_{k}\bm{b}_{k}\bm{b}_{k}{}^{\prime}$ is close to $\bar{\bm{\Lambda}}$ whp. This will hold only when $q$ is large enough, and, hence, our result will also need another lower bound on $q$. The second change is that we need to replace $r\rho\bar{\lambda}_{\max}$ by $r(10\log n)\bar{\lambda}_{\max}$ in the lower bound on $mq$. This is the high probability upper bound on $\|\bm{b}_{k}\|^{2}$ under Model 3.4. Moreover, because of these two changes, the probability of the good event reduces slightly.

In the setting of Theorem 3.2, suppose that the $\bm{x}_{k}$’s satisfy Model 3.4. For a $\varepsilon<1$, if

then, the conclusions of Theorem 3.2 hold with probability at least $1-2\exp(-cn)-\frac{36q}{n^{4}}-\frac{20}{n^{2}}$.

As will be evident from the proof of Theorem 3.5, any random model that ensures that (a) $\max_{k}\|\bm{b}_{k}\|^{2}$ is bounded whp, and (b) $\frac{1}{q}\sum_{k}\bm{b}_{k}\bm{b}_{k}{}^{\prime}$ is close to $\bar{\bm{\Lambda}}$ whp will suffice. For example, even if the $\bm{b}_{k}$’s in Model 3.4 have nonzero and different means, a similar result can be proved. More generally, as we state below, a sub-Gaussian assumption works as well. The independence assumption on $\bm{b}_{k}$’s may also be weakened to any other assumption that ensures that (b) holds, however we do not pursue it here.

With this model replacing Model 3.4 on $\bm{X}$, Theorem 3.5 holds with probability at least $1-2\exp(-cn)-\frac{cq}{n^{4}}-\frac{c}{n^{2}}.$

In the proof of Theorem 3.5, only the proofs of Lemmas 6.8, 6.9 change.∎

III-C Main Results - Unknown rank case

We now turn to the setting where the rank $r$ is unknown and show how Theorem 3.2 can be modified for this setting. Other results are modified similarly.

Consider the rank estimation approach given in Algorithm 3. We have the following corollary.

Consider Algorithm 3 (unknown $r$ case). Assume the setting of Theorem 3.2 with $\varepsilon\leq 0.001$. If, in addition, $\kappa\leq 10$ and if $\bar{\bm{\Lambda}}$ is such that $\bar{\lambda}_{j}-\bar{\lambda}_{j+1}\leq 0.9\bar{\lambda}_{\min}$, then, with the probability given in Theorem 3.2,

Another way to correctly estimate $r$ is via thresholding.

Consider Algorithm 3 with rank estimated as follows. Set $\hat{r}$ as the smallest index $j$ for which $\lambda_{j}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})\geq 0.25\bar{\lambda}_{\min}$. Assume the setting of Theorem 3.2 with $\varepsilon\leq 0.001$. Then, if $\kappa<124$, then, with the probability given in Theorem 3.2,

The rank estimation approach of Algorithm 3 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 only when consecutive eigenvalues of $\bar{\bm{\Lambda}}$ (consecutive nonzero singular values of $\bm{X}$) are not too far apart. On the other hand, the thresholding based approach of Corollary 3.8 does not require any extra assumptions beyond those in Theorem 3.2 and $\kappa<124$. However it necessitates knowledge of $\bar{\lambda}_{\min}$.

IV Low Rank PR (LRPR) - Complete algorithm

So far we developed an initialization procedure that directly exploited the low-rank property of $\bm{X}$. Here, we explain three possible ways to develop a complete LRPR algorithm. The first, given next, uses LRPR-init to only speed up TWF.

If $q\geq c\sqrt{n}$, then this means that $\varepsilon=c\frac{\rho\kappa r^{2}(\log n)}{n^{1/4}}$ works. Combining this with [11, Theorem 1], we have the following corollary.

IV-B LRPR1: Low Rank PR via projected gradient descent

The simplest way to develop a complete algorithm that exploits the low rank property of $\bm{X}$ is to use a projected gradient descent approach to modify TWF. This projects the TWF output at each iteration onto the space of rank $r$ matrices. We summarize the complete LRPR1 approach (projected-TWF initialized with LRPR-init) in Algorithm 5. When $m$ is small, this results in significantly improved performance over TWF because it exploits the low-rank structure of the matrix $\bm{X}$ at each step. For an example, see Fig. 2(b).

IV-C LRPR2: Low Rank PR via Alternating Minimization

The third and most powerful approach is to modify the entire algorithm to directly exploit the low-rank property of the matrix $\bm{X}$, i.e., to use its decomposition as $\bm{X}={\bm{U}}\bm{B}$. This idea can be used to modify TWF or AltMinPhase (Gerchberg-Saxton algorithm) or, in fact, many of the other PR methods from literature, e.g., . As noted by an anonymous reviewer, the last two are significantly faster than Gerchberg-Saxton. TWF is truncated gradient descent to minimize the negative data likelihood under a Poisson noise assumption, where as AltMinPhase is an AltMin approach to minimize the squared loss function (data likelihood under iid Gaussian noise). For noise-free measurements, this distinction is immaterial, and all methods apply.

Modifying TWF for the set of variables ${\bm{U}},\bm{B}$ needs to be done with care, and needs to include a step that ensures that one of $\|{\bm{U}}\|$ or $\|\bm{B}\|$ does not keep increasing. An early attempt along these lines is given in .

We show the power of both LRPR1 and LRPR2 for recovering a real video from coded diffraction pattern (CDP) measurements in Fig. 1. As can be seen, with as few as $m=3n$ CDP measurements, both these methods significantly outperform basic TWFproj (Algorithm 5 initialized using Algorithm 2) and basic TWF (Algorithm 4 initialized using Algorithm 7). This experiment is inspired by an analogous experiment for recovering a regular camera image from CDP measurements reported in [11, Fig. 2]. While this is not a real practical application since the video used is a regular camera video of a moving airplane, this example illustrates two points: (i) many real image sequences are indeed approximately low-rank; and (ii) our algorithm has significant advantage over single vector PR methods for jointly recovering this approximately low-rank video. For a detailed explanation of this and some more such experiments, please see Supplementary Material and http://www.ece.iastate.edu/~namrata/LRPR/.

V Numerical Experiments

We discuss here the results of three sets of experiments. All experiments were done on a single laptop which had these specifications: Intel(R) CPU E3-1240 v5 3.50 GHz, Installed memory: 32 GB, System type is 64 bit.

Experiment 1. The first experiment shows the power of the proposed initialization approach, LRPR-init (Algorithm 1), by comparing its initialization error with that of TWF initialization (TWF-init, Algorithm 7) and of TWFproj-init (Algorithm 2). TWF-init does not use knowledge of rank, TWFproj-init assumes $r$ is known, while LRPR-init estimates the rank automatically as explained earlier. For a fair comparison with TWFproj-init, we also show the error of LRPR-init with $\hat{r}=r$. Data was generated as follows. The matrix ${\bm{U}}$ is obtained by orthonormalizing an $n\times r$ matrix with iid Gaussian entries; $\bm{b}_{k}$’s were generated as being iid uniformly distributed between $-1$ and $1$; and we set $\bm{x}_{k}={\bm{U}}\bm{b}_{k}$. Measurements were generating using (1).

When the product $mq$ is large, the rank $\hat{r}$ is correctly estimated by LRPR-init (Algorithm 1) either always or most of the time. We display a Monte Carlo estimate of the probability of $\hat{r}=r$ in the 2nd column. In these cases, LRPR with $\hat{r}$ known versus $\hat{r}$ estimated both have similar errors (3rd and 4th columns). Inspired by a reviewer’s concern, we also evaluate LRPR-init with $\hat{r}$ deliberately set to a wrong value $2r$ in the 5th column. As can be seen, the error degradation is gradual even with a wrong rank estimate.

Finally, Table I shows errors of LRPR-Same in the last column. This refers to LRPR operating on measurements of the form $\bm{y}_{i,k}:=(\bm{a}_{i}{}^{\prime}\bm{x}_{k})^{2}$. Because it uses the same $\bm{a}_{i}$’s for all columns $\bm{x}_{k}$, there are only $m$ (and not $mq$) mutually independent matrices to average over. Hence its errors are almost as large as those of TWF.

In Fig. 2(a), we compare the speed of error decay of TWF when initialized with either TWF-init (TWF) or with the proposed initialization, LRPR-init (LRPR+TWF). We used $m=8n$ (large enough $m$ for TWF iterations to converge). For $t=0,1,2,\dots,100$, we plot the error at the end of iteration $t$ on the y-axis and the time taken till the end of iteration $t$ on the x-axis ($t=0$ corresponds to initialization). As can be seen, LRPR-init takes longer to finish than TWF-init (the first ‘triangle’ is to the right of the first circle). However, because LRPR-init results in much lower initialization error, LRPR+TWF needs much fewer iterations to “converge”, and, so the total time taken by it to “converge” is also smaller.

If $m$ is reduced to $m=0.8n$ measurements, as can be seen from Fig. 2(b), neither of TWF or LRPR+TWF converge. Basic TWFproj also does not converge and this is because its initialization error is larger (for reasons explained earlier). However, both LRPR1 and LRPR2 converge. It is also apparent that LRPR1 is significantly faster than LRPR2. This is because its per iteration cost is lower.

If $m$ is reduced further to $m=0.6n$ (Fig. 2(c)), then LRPR1 does not converge whereas LRPR2 still does. This is because LRPR2 iterates directly exploit the split-up $\bm{X}={\bm{U}}\bm{B}$ whereas LRPR1 iterates first implement a TWF iteration and then project the resulting matrix onto the space of rank $r$ matrices.

VI Proofs of Theorems 3.2 and 3.5

The derivations in this section use many useful results about sub-Gaussian and sub-exponential r.v.’s and the $\epsilon$-net taken from . These are summarized in Appendix A. The lemmas that are not proved here are proved in Appendix B.

We first state a simple corollary of the Davis-Kahan $\sin\theta$ theorem [28, Sec. 2] that follows from it using Weyl’s inequality (see for a proof).

Consider a Hermitian matrix $\bm{D}$ and its perturbed version $\hat{\bm{D}}$. Define $\bm{H}:=\hat{\bm{D}}-\bm{D}$. Let $\bm{E}$ be the matrix of top $r$ eigenvectors of $\bm{D}$, and let $\bm{F}$ be the matrix of top $r$ eigenvectorsMore generally, $\bm{E}$ and $\bm{F}$ can be any matrices whose columns span the space of top $r$ eigenvectors of $\bm{D}$ and $\hat{\bm{D}}$ respectively. of $\hat{\bm{D}}$. If $\lambda_{r}(\bm{D})-\lambda_{r+1}(\bm{D})-\|\bm{H}\|>0$, then

In Sec. VI-B, we will use the above result with $\hat{\bm{D}}=\bm{Y}_{U}$ and $\bm{D}$ being the expected value of a matrix that is close to it. In Sec. VI-C, we will use it similarly for $\bm{Y}_{b,k}$.

Theorem 6.2 below is a simple generalization of Theorem 5.39 of .

Suppose that $\bm{w}_{j}$, $j=1,2,\dots,N$, are $n$-length independent, sub-Gaussian random vectors with sub-Gaussian norms bounded by $K$.

For an $\varepsilon<1$ and a given vector $\bm{z}$, with probability (w.p.) $\geq 1-2\exp(-c\varepsilon^{2}N)$,

For an $\varepsilon<1$, w.p. $\geq 1-2\exp(n\log 9-c\epsilon^{2}N)$,

The proof follows that of Theorem 5.39 in . It is given in the Supplementary Material. ∎

The matrix $\bm{X}$ is a deterministic unknown.

This just makes it simpler to simultaneously obtain subspace error bounds under the assumptions of both Theorems 3.2 and 3.5. Notice that, if we write $\bm{x}_{k}={\bm{U}}\bm{b}_{k}$, then the definitions of $\bar{\bm{\Lambda}}$, $\rho$ and $\kappa$ given in (6) and (7) in Sec. III-A imply that, under Model 6.3, $\bar{\bm{\Lambda}}=\frac{1}{q}\sum_{k}\bm{b}_{k}\bm{b}_{k}{}^{\prime}$, $\kappa$ is its condition number, and $\max_{k}\|\bm{b}_{k}\|^{2}=\max_{k}\|\bm{x}_{k}\|^{2}\leq r\rho\bar{\lambda}_{\max}$.

Thus, all we need now is to specify $\bm{\Sigma}^{-}$ and find a high probability upper bound on $\|\bm{Y}_{U}-\bm{\Sigma}^{-}\|$.

To this end, as also done in [11, Appendix C], we first lower and upper bound $\bm{Y}_{U}$ in order to replace $\frac{1}{m}\sum_{i}\bm{y}_{i,k}$ in its indicator function expression by a constant. Recall that $\bm{Y}_{U}$ is defined in (3) and that $\frac{1}{m}\sum_{i}\bm{y}_{i,k}=\bm{x}_{k}{}^{\prime}(\frac{1}{m}\sum_{i}\bm{a}_{i,k}\bm{a}_{i,k}{}^{\prime})\bm{x}_{k}$. By Fact A.3, item 3, in Appendix A, $\bm{a}_{i,k}$’s are sub-Gaussian with sub-Gaussian norm bounded by $c$. Thus, using the first part of Theorem 6.2, conditioned on $\bm{x}_{k}$, $|\frac{1}{m}\sum_{i}\bm{y}_{i,k}-\|\bm{x}_{k}\|^{2}|\leq\epsilon_{1}\|\bm{x}_{k}\|^{2}$ w.p. $\geq 1-2\exp(-c\epsilon_{1}^{2}m)$. The constant multiplying $\epsilon_{1}$ is moved into the $c$ in the probability. This bound holds for all $k=1,2,\dots,q$ w.p. $\geq 1-2q\exp(-c\epsilon_{1}^{2}m)$. This implies that, with the same probability, conditioned on $\bm{X}$, $\bm{Y}^{-}\preceq\bm{Y}_{U}\preceq\bm{Y}^{+},\ \text{where}$

Notice that $\bm{w}_{i,k}^{+}$ is $\bm{w}_{i,k}^{-}$ with $9(1-\epsilon_{1})$ replaced by $9(1+\epsilon_{1})$ in the indicator function. The following claim is immediate.

Conditioned on $\bm{X}$, w.p. $\geq 1-2q\exp(-c\epsilon_{1}^{2}m)$, $\|\bm{Y}_{U}-\bm{Y}^{-}\|\leq\|\bm{Y}^{+}-\bm{Y}^{-}\|$.

We obtain expressions for these in the next lemma.

By the triangle inequality and Lemma 6.4, conditioned on $\bm{X}$, w.p. $\geq 1-2q\exp(-c\epsilon_{1}^{2}m)$,

Conditioned on $\bm{X}$, w.p. $\geq 1-2\exp(n\log 9-\epsilon_{2}^{2}mq)$,

Since the claim of Lemma 6.6 holds with the same probability lower bound for all $\bm{X}$, it also holds with the same probability lower bound if we average over $\bm{X}$. The same is true for Lemma 6.4.

The next lemma bounds $\max_{k}\|\bm{b}_{k}\|^{2}=\max_{k}\|\bm{x}_{k}\|^{2}$.

Under Model 6.3, $\max_{k}\|\bm{b}_{k}\|^{2}\leq r\rho\bar{\lambda}_{\max}$. Under Model 3.4, w.p. $\geq 1-2q/n^{4}$, $\max_{k}\|\bm{b}_{k}\|^{2}\leq r(10\log n)\bar{\lambda}_{\max}$.

Combining Lemmas 6.6, 6.8 and 6.9, we can bound $\|\bm{Y}^{-}-\bm{\Sigma}^{-}\|$ under both models. We get the same bound on $\|\bm{Y}^{+}-\bm{\Sigma}^{+}\|$ as well.

Finally, we bound $\|\bm{\Sigma}^{+}-\bm{\Sigma}^{-}\|$ using the fact that, for all $\xi>\xi_{0}$, $\xi^{d}e^{-\frac{\xi^{2}}{4}}\leq b$ for any $d>1$.

Combining the above bounds, we conclude the following.

under Model 6.3, w.p. $\geq 1-p_{U,1}$, $\|{\bm{Y}}_{U}-\bm{\Sigma}^{-}\|\leq\epsilon_{U,1}\bar{\lambda}_{\max}$; and

under Model 3.4, w.p. $\geq 1-p_{U,2}$, $\|{\bm{Y}}_{U}-\bm{\Sigma}^{-}\|\leq\epsilon_{U,2}\bar{\lambda}_{\max}.$

Finally, to bound the subspace error of ${\bm{\hat{U}}}$, we also require a lower bound on $\beta_{1}^{-}$. This follows easily using the fact that, for all $\xi>\xi_{0}$, $\xi^{d}e^{-\frac{\xi^{2}}{4}}\leq b$ for any $d>1$.

Applying the $\sin\theta$ theorem, Theorem 6.1, and the last two claims above, we get the following result.

Let Model 1 be Model 6.3 and Model 2 be Model 3.4. If $\kappa\epsilon_{U,d}<1/16$ then, under Model $d$, w.p. $\geq 1-p_{U,d}$,

Since $\epsilon_{1}\leq\epsilon_{U,d}\leq\kappa\epsilon_{U,d}\leq 1/16$, using Lemma 6.12, $\beta_{1}^{-}\geq 0.5$. Thus, under Model $d$, w.p. $\geq 1-p_{U,d}$,

For notational simplicity, we let $\bm{x}=\bm{x}_{k}$, $\bm{b}:=\bm{b}_{k}$, $\bm{y}_{i}:=\bm{y}_{i,k}^{new}$, $\bm{a}_{i}:=\bm{a}_{i}^{new}$, ${\bm{Y}}_{\bm{b}}:={\bm{Y}}_{\bm{b},k}$ defined in (5) in Algorithm 3. Since the different $\bm{b}_{k}$’s are recovered separately, but using the same technique, this notation does not cause any confusion at most places. Where it does, we clarify.

In this section, we state all results conditioned on $\bm{x}$ and ${\bm{\hat{U}}}$. Under this conditioning, in all our claims, the probability of the desired event is lower bounded by a value that does not depend on $\bm{x}$ or ${\bm{\hat{U}}}$. Thus, the same probability lower bound holds even when we average over $\bm{x}$ and ${\bm{\hat{U}}}$ (holds unconditionally).

Notice that $\bm{x}={\bm{U}}\bm{b}$ can be rewritten as

We can further split $\bm{g}$ as $\bm{g}=\bm{v}\nu$ where $\nu=\|\bm{g}\|$ and $\bm{v}=\bm{g}/\nu$. Recall that we estimate $\bm{x}$ as

The vector $\bm{e}$ and the scalar $\nu:=\|\bm{g}\|$ satisfy $\|\bm{e}\|\leq\delta_{U}\|\bm{b}\|$, andusing $\nu=\|{\bm{\hat{U}}}\bm{g}\|\geq\|\bm{x}\|-\|\bm{e}\|\geq(1-\delta_{U})\|\bm{b}\|$, $(1-\delta_{U})\|\bm{b}\|\leq\nu\leq\|\bm{b}\|$.

Clearly, the top eigenvector of this matrix is proportional to $\bm{g}$ and the desired eigen-gap is $2\|\bm{g}\|^{2}\geq 2(1-\delta_{U})^{2}\|\bm{b}\|^{2}$. So, we can use $\bm{D}=2\bm{g}\bm{g}^{\prime}+\|\bm{g}\|^{2}\bm{I}$ for applying Theorem 6.1. It remains to bound $\|\bm{Y}_{b}-(2\bm{g}\bm{g}^{\prime}+\|\bm{g}\|^{2}\bm{I})\|$.

To bound $\|\bm{Y}_{g}-(2\bm{g}\bm{g}^{\prime}+\|\bm{g}\|^{2}\bm{I})\|$, we use following modification of Theorem 4.1 of .

The value of $\|\bm{Y}_{e,2}\|$ can be bounded in a similar fashion:

Using Theorem 6.1 and Fact 6.15, with the same probability,

If the numerator is smaller than $1-2\delta_{U}$, then

As before, we can average over $\bm{X}$ and ${\bm{\hat{U}}}$ and still get all the events above to hold with the same probability. Using the above bound, (12), and Lemma 6.16, we get the following.

VI-D Proof of Theorem 3.2

Using the expressions for $p_{U,1}$ and $p_{g}$, to get the probability of the desired event below $1-2\exp(-cn)-32q/n^{4}$, we need

VI-E Proof of Theorem 3.5

We use the same approach as above. With Model 3.4, both $\epsilon_{U}$ and $p_{U}$ are larger. We have $\epsilon_{U}=\epsilon_{U,2}$ which is equal to $3\epsilon_{3}$ plus $\epsilon_{U,1}$ with $\rho$ replaced by $10\log n$. Also, $p_{U}=p_{U,2}=p_{U,1}+2q/n^{4}+2\exp(r\log 9-c\epsilon_{3}^{2}q)+18\exp(-c\epsilon_{3}^{2}\frac{q}{r})$.

VI-F Proof of Corollary 3.7 and Corollary 3.8

Using Corollary 6.11, w.p. $\geq 1-p_{U,1}$,

and $\beta_{1}^{-}(\epsilon_{1}),\beta_{2}^{-}(\epsilon_{1})$ are defined in Lemma 6.5. Let $\beta_{1}^{-}:=\beta_{1}^{-}(\epsilon_{1})$. Suppose that $\varepsilon\leq 0.001$. From Sec. VI-D, $\epsilon_{1}\leq\epsilon_{U,1}\leq\varepsilon\leq 0.001$. Using $\kappa\leq 10$, $2\kappa\epsilon_{U,1}\leq 0.02$.

Thus, from (17), Weyl’s inequality , and $\bar{\lambda}_{j}-\bar{\lambda}_{j+1}\leq 0.9\bar{\lambda}_{\min}$, we conclude the following: for a $j<r$ and a $j^{\prime}>r$, w.p. $\geq 1-p_{U,1}$,

By Lemma 6.12, $\beta_{1}^{-}\geq 0.5$. Using this, we conclude that

for a $j<r$ and a $j^{\prime}>r$. The second row used $0.1\beta_{1}^{-}>0.04$ and the last row used $\beta_{1}^{-}-0.02>0.02$.

Thus, if $\kappa\leq 10$, and $\bar{\lambda}_{j}-\bar{\lambda}_{j+1}\leq 0.9\bar{\lambda}_{\min}$, then, under the assumptions of Theorem 3.2, w.p. $\geq 1-p_{U,1}$, $\lambda_{j}(\bm{Y}_{U})-\lambda_{j+1}(\bm{Y}_{U})$ is largest for $j=r$, i.e., $\hat{r}=r$. Using this and then proceeding exactly as before, we obtain Corollary 3.7.

To get Corollary 3.8, use (17), Weyl’s inequality , and $\kappa\leq 124$, to argue that, w.p. $\geq 1-p_{U,1}$, for any $j>r$,

Thus, w.p. $\geq 1-p_{U,1}$, $j=r$ is the smallest index for which $\lambda_{r}(\bm{Y}_{U})-\lambda_{n}(\bm{Y}_{U})\geq 0.25\bar{\lambda}_{\min}$ and hence the rank estimation approach of Corollary 3.8 returns $\hat{r}=r$.

VII Conclusions and Future Work

We obtained sample complexity bounds for LRPR-init and argued that, when $q/r$ is large, these are much smaller than those for TWF or any other single-vector PR method. Via extensive experiments, we also showed that the same is true for both the complete algorithms - LRPR1 and LRPR2. Between the two, LRPR2 has better performance, but also higher per iteration computational cost, than LRPR1.

In future work we will analyze the complete LRPR2 algorithm. This should replace the dependence of sample complexity on $1/\varepsilon^{2}$ by a dependence on $-\log\varepsilon$.

References

Appendix A Preliminaries

As explained in , $\epsilon$-nets are a convenient means to discretize compact metric spaces. The following definition is [30, Definition 5.1] for the unit sphere.

By Lemma 5.4 of , for a symmetric matrix, $\bm{W}$, $\|\bm{W}\|=\max_{\bm{x}:\|\bm{x}\|=1}\|\bm{x}^{\prime}\bm{W}\bm{x}\|\leq\frac{1}{1-2\epsilon}\max_{\bm{x}\in\mathcal{N}_{\epsilon}}\|\bm{x}^{\prime}\bm{W}\bm{x}\|.$

By [30, Corollary 5.17], if $x_{i}$, $i=1,2,\dots N$, are a set of independent, centered, sub-exponential r.v.’s with sub-exponential norm bounded by $K_{e}$, then, for an $\varepsilon<1$,

If $\bm{x}\sim\mathcal{N}(0,\bar{\bm{\Lambda}})$ with $\bar{\bm{\Lambda}}$ diagonal, then $\bm{x}$ is sub-Gaussian with $\|\bm{x}\|_{\varphi_{2}}\leq c\sqrt{\bar{\lambda}}_{\max}$. Moreover, if $\bm{y}=x_{1}\bm{x}$ where $x_{1}$ is a zero mean bounded r.v. with bound $M$, then $\|\bm{y}\|_{\varphi_{2}}\leq cM\sqrt{\bar{\lambda}_{\max}}$.

If $\bm{x}_{i}\sim\mathcal{N}(0,\bar{\bm{\Lambda}})$, for $i=1,2,\dots,N$, are $n$-length random vectors and $\bar{\bm{\Lambda}}$ is diagonal, then

This is a direct consequence of eq. 5.5 of which says that if $x\sim\mathcal{N}(0,1)$, then $\Pr(|x_{i}|>t)\leq 2\exp(-t^{2}/2)$ for a $t>1$. Using this along with the union bound first for bounding $\|\bm{x}_{i}\|^{2}=\sum_{j=1}^{n}(\bm{x}_{i})_{j}^{2}$ for a given $i$ and then for bounding its $\max$ over $i$ gives the above result.

Using [30, Lemma 5.5], if $\bm{x}_{i}$’s are sub-Gaussian random vectors with sub-Gaussian norm bounded by $K$, then the following generalization of the above fact holds: $\Pr\left(\max_{i=1,2,\dots,N}\|\bm{x}_{i}\|^{2}\leq K^{2}\cdot n\cdot 2\nu\right)\geq 1-CnN\exp(-c\nu).$

The following is an easy corollary of Cauchy-Schwartz for sums of products of vectors.

Appendix B Proofs of lemmas from Section VI

We prove the lemmas that were not proved in Section VI.

Then $\bm{Y}^{-}=\frac{1}{mq}\sum_{k}\sum_{i}\bm{w}_{i,k}\bm{w}_{i,k}{}^{\prime}$.

Let $D=\max_{k}\|\bm{b}_{k}\|^{2}=\max_{k}\|\bm{x}_{k}\|^{2}$.

We first argue argue that, conditioned on $\bm{X}$, each $\bm{w}_{i,k}$ is sub-Gaussian with sub-Gaussian norm bounded by $c\|\bm{x}_{k}\|\leq c\sqrt{D}$.

To show this easily, we use the strategy of [11, Appendix C]. Since $\bm{a}_{i,k}$ is rotationally symmetric, without loss of generality, suppose that $\frac{\bm{x}_{k}}{\|\bm{x}_{k}\|}$ is the first column of the identity matrix. Then, $\bm{w}_{i,k}=\bm{a}_{i,k}(\bm{a}_{i,k})_{1}\mathds{1}_{(\bm{a}_{i,k})_{1}^{2}\leq 9(1-\epsilon_{1})}\ \|\bm{x}_{k}\|$. With this simplification, $\bm{w}_{i,k}$ is of the form $\bm{x}x_{1}$ where $x_{1}$ is a bounded r.v. with bound $\sqrt{9(1-\epsilon_{1})}\|\bm{x}_{k}\|$ and $\bm{x}$ is Gaussian with zero mean and covariance matrix $\bm{I}$. Thus, using Fact A.3, item 3, it is sub-Gaussian with sub-Gaussian norm bounded by $c\|\bm{x}_{k}\|\leq c\sqrt{D}$.

Conditioned on $\bm{X}$, all the $\bm{w}_{i,k}$’s are mutually independent. There are $N=mq$ of them.

Let $D=\max_{k}\|\bm{b}_{k}\|^{2}$. Under Model 6.3, $D\leq r\rho\bar{\lambda}_{\max}$ by definition. Under Model 3.4, we use Fact A.3, item 4 with $n\equiv r$, $N\equiv q$, $\nu\equiv 5\log n$ to get $D\leq 10\log n$ w.p. $\geq 1-2rqn^{-5}\geq 1-2qn^{-4}$ since $r\leq n$.∎

Using $\xi\exp(-\xi^{2}/2)<1$ for all $\xi^{2}>8$,

Similarly, using $\xi^{3}\exp(-\xi^{2}/2)<3.1$ for all $\xi^{2}>8$,

Using $(\xi^{4}-\xi^{2})\exp(-\xi^{2}/4)<7.58$ for all $\xi^{2}>8$, Fact A.3, item 4 and $\epsilon_{1}\leq 1/9$,

Appendix C Supplementary Document

Since $\mathcal{N}_{1/4}$ is a finite set of vectors, all we need to do now is to bound $|\bm{z}^{\prime}\bm{W}\bm{z}|$ for a given vector $\bm{z}$ followed by applying the union bound to bound its maximum over all ${\bm{z}\in\mathcal{N}_{1/4}}$. The former has already been done in the first part. By Fact A.2 (Lemma 5.2 of ), the cardinality of $\mathcal{N}_{1/4}$ is at most $9^{n}$. Thus, using the first part, $\Pr\left(\max_{\bm{z}\in\mathcal{N}_{1/4}}|\bm{z}^{\prime}\bm{W}\bm{z}|\geq\frac{4\varepsilon K^{2}}{2}\right)\leq 9^{n}\cdot 2\exp(-c\frac{\varepsilon^{2}}{4}N)=2\exp(n\log 9-c\varepsilon^{2}N).$ By (18), we get the result. ∎

The proof is a simplified and clearer version of the proof of Theorem 4.1 of . The few differences are as follows: we truncate differently (in a simpler fashion); and we use different constants to get a higher probability of the good event.

Next consider the $(j,j)$-th entry for $j>1$. This can be bounded using a similar trick.

Appendix D Experiment details for Fig. 1

We used real videos that are approximately low rank and CDP measurements of their images. Each image (arranged as a 1D vector) corresponds to one $\bm{x}_{k}$ and hence the entire video corresponds to the matrix $\bm{X}$. We show results on a moving mouse video and on a moving airplane video (shown in Fig. 1). We show two results with “low-rankified videos”The original video data matrix $\bm{X}_{orig}$ was made exactly low rank by projecting it onto the space of rank-$r$ matrices where $r$ was chosen to retain 90% of the singular values’ energy. and one result with the original airplane video. The airplane images were of size $n_{1}\times n_{2}$ with $n_{1}=240$, $n_{2}=320$; the mouse images had $n_{1}=180$, $n_{2}=319$. Thus, $n=n_{1}n_{2}=76800$ and $n=57420$ respectively. Mouse video had $q=90$ frames and airplane one had $q=105$ frames.

The CDP measurement model can be understood as follows . First, note that it allows $m$ to only be an integer multiple of $n$; so let $m=nL$ for an integer $L$. Let $\bm{y}_{k}$ denote the vector containing all measurements of $\bm{x}_{k}$. Then $\bm{y}_{k}=|\bm{A}_{k}{}^{\prime}\bm{x}_{k}|^{2}$ where $\bm{A}_{k}=[(\bm{F}\bm{M}_{k,1})^{\prime},(\bm{F}\bm{M}_{k,2})^{\prime},\dots,\bm{F}(\bm{M}_{k,L})^{\prime}]$; each $\bm{M}_{k,l}$ is a diagonal $n\times n$ mask matrix with diagonal entries chosen uniformly at random from the set $\{1,-1,\sqrt{-1},-\sqrt{-1}\}$, and $\bm{F}=\bm{F}_{1D,n_{1}}\otimes\bm{F}_{1D,n_{2}}$ where $\bm{F}_{1D,n}$ is the $n$-point discrete Fourier transform (DFT) matrix and $\otimes$ denotes Kronecker product. Thus, $(\bm{F}\bm{x}_{k})$ is the vectorized version of the 2D-DFT of the image corresponding to $\bm{x}_{k}$

In this experiment, $n$ and $m$ are very large and hence the memory complexity is very large. Thus, the algorithm cannot be implemented using matrix-vector multiplies. However, since the measurements are masked-Fourier, we can implement its “operator” version as was also done in the TWF code . All matrix-vector multiplies are replaced by “operators” that use 2D fast Fourier transform (2D-FFT) or 2D-inverse-FFT (2D-IFFT) functions, preceded or followed by applying the measurement masks. This is a much faster and memory efficient implementation. Only the masks need to be stored. The EVD in the initialization step is implemented by a block-power method that uses 2D-FFT. The LS step is implemented using the operator-version of conjugate gradient LS (CGLS) taken from http://web.stanford.edu/group/SOL/software/cgls/. TWFproj and LRPR1 are implemented similarly.