Improved algorithm and bounds for successive projection

I have one comment about this review (I'm a reviewer, and not the original author)

This paper seems to be a incremental work based on Gillis & Vavasis

This is a completely incorrect assessment. The proofs are clearly not the same. Working on the same problem/algorithm is not incremental work no matter how you spin it.

Thank you for your reply.

First, I am not totally satisfied with the explanation on the matrix completion.

In the field of the matrix completion research, ppl do not discuss whether $V$ equals $\tilde{V}$ but whether they are the same under some rotations. This is because we always have the relation $VO O^T \Pi = V \Pi$, where $O$ is an orthonormal matrix.

I suspect that this problem can be solved in the framework of matrix completion. All you need is to align the matrix to the vertices. This is not very hard. For example, you can use the fact that the points lies in the vertices, i.e., the linear combination coefficients are sparse.

Also, I am not satisfied with the explanation that sample number $n$ does not play any role. Your paper considers a noisy sensing relation. Even under the duplicate measurements setting, you can improve the performance by averaging the duplicate measurement.

For example, suppose you have $n$ measurements $r_i = v + \epsilon_i$, where $\epsilon_i$ is zero mean and $\sigma^2$ variance. Averaging $n$ cases will reduce the sensing noise variance to $\sigma^2/n$. In this way, I am a little surprised that the measurement plays no role. Or this is an artifact of this line of algorithm?

Dear reviewer KsqF, I do not want to get into such an argument on the issue such whether I spin this work or not. Still, I believe I am allowed to make an independent decision and am qualified to make one.

I reviewed the paper carefully and acted in a responsible way. I have provided some unique angles to this paper. Even the authors themselves agree with me, as you can see in their own words.

First, I pointed out that this problem can be studied in the framework of matrix completion. Even the authors themselves agree it is a future research direction they want to explore. I quote their words here "there may be a way to borrow some ideas from matrix completion literature to further improve the current algorithms. We leave this interesting direction to future work."

Second, I asked about the role of the sample number, a question seldom asked. And the author's response is not satisfactory, as I have pointed out in my response. For example, even with duplicate measurements, more measurements can be used to reduce the measurement error. This does not justify the missing role of the sample number in the reconstruction performance. The authors' responses do not address this.

I am open to changing my ideas if some reasonable arguments can be provided. Thank you.

Thank you for your comments and for appreciating that "this paper proposes a novel algorithm" and "the analysis seems to be solid". We already posted a quick response earlier. This current response is to supplement the earlier one with more details. We also encourage you to read our Response-to-All-Reviewers for other clarifications.

Main Points:

1. Comparison with [GV13]: "This paper seems to be an incremental work based on Gillis & Vavasis (published in 2013)."

We are sorry to hear this comment, but we believe it is due to some misunderstanding. It is our fault that our initial submission did not explain why the proof is non-trivial. We hope that you will change your mind after seeing the clarification here.

Our paper has 3 major contributions:

• (a) A new algorithm
• (b) the non-asymptotic bound for SPA,
• (c) the error bound for pp-PSA.

Only (b) is related to Gillis & Vavasis (2013) [GV13]. Even without (b), the contributions in (a) and (c) are already significant. For example, our analysis in pp-SPA uses theory of extreme order statistics to analyze the effect of KNN denoise, which is subtle.

Regarding the contribution in (b), we guess you were actually asking: "Is the proof a simple modification of the proof in [GV13]?"

The proof is completely different. The proof in [GV13] is based on controlling matrix norms. Their proof does not use explicitly the geometry behind this problem. To see this, we recall that $\beta(X)$ is equal to the maximum of $\Vert\epsilon_i\Vert$ among all points, and $\beta_{\text{new}}(X)$ is roughly the maximum of $\Vert\epsilon_i\Vert$ among those points outside the true simplex. In the proof of [GV13], they frequently control the difference of replacing $X$ by $\mathbb{E}[X]$ in many intertermediate matrices. Such steps are feasible only if one has a bound for the maximum of all $\Vert\epsilon_i\Vert$'s - i.e., the quantity $\beta(X)$. However, it is hard to repeat these steps by using $\beta_{\text{new}}(X)$. Similarly, if one has a lower bound for $s_K(V)$, one can control the inverse of $V'V$, but a lower bound for $s_{K-1}(V)$ cannot do so. To tackle these difficulties, we use a different proof approach. Our proof relies more on the geometry behind this problem. We define a "simplicial neighborhood" for each true vertex and prove that each estimated vertex must fall into one of these simplicial neighborhoods. This is a new proof strategy.

In addition, we have added a numerical example in this revised draft, where we calculate the two bounds explicitly. In some cases, the new bound is only 0.01 times the old bound in [GV13].

In summary, we agree that [GV13] is an influential work, but our paper is not incremental of [GV13].

2. More simulations: "There is only limited numerical experiments to validate the performance of the proposed algorithms"

Thanks for your comment. We have significantly expanded our simulations:

• We have included settings for every $d=1,2,..., 48, 49, 59$.
• We have added the comparison with robust SPA in all simulation settings.

The updated Figure 3 contains nearly 90 different settings. We hope that the numerical experiments are sufficiently rich now.

Other questions:

1. The relation between sample number $n$ and the reconstruction error.

The reconstruction error is $O(\sqrt{d+\log(n)})$ for SPA and $O(\sqrt{\log(n)})$ for pp-SPA. In addition, in some special cases (e.g., Theorem 2), the error of pp-SPA can be further reduced to nearly $O(1)$.

At first glance, it is a bit counter-intuitive that increasing $n$ does not decrease the rate. The main reason is that the configuration of population points $r_1, r_2, \ldots, r_n$ can be arbitrary. Here is a toy example. Suppose

In this example, no matter how large $n$ is, there is only one point, $X_k$, that is useful for estimating $v_k$. This is why increasing $n$ does not help reduce the error bound.

This insight is further supported by simulations. On the right panel of Figure 3, we plot the error curve as $n$ increases. It shows that the error does not decrease with $n$, which is consistent with our theory.

2. How is the algorithm compared with the matrix completion solution?

This is an interesting comment. The main difference is: We are interested in recovering $V$, but the matrix completion literature is interested in recovering $V\Pi$. We can find $\widetilde{V}\in\mathbb{R}^{d\times K}$ and $\widetilde{\Pi}\in\mathbb{R}^{K\times n}$ such that $\widetilde{V}\widetilde{\Pi}=V\Pi$ but $\widetilde{V}\neq V$.

For example, let's consider a special case with $K=3$ and $d=2$. The true simplex is a triangle in $\mathbb{R}^2$. We now consider an arbitrary triangle that contains the true triangle in its interior. Let $\widetilde{V}$ contain vertices of this new triangle, and let $\widetilde{\Pi}$ contain the barycentric coordinates of $r_i$'s in this new triangle. Then, $V\Pi=R=\widetilde{V}\widetilde{\Pi}$. However, $V\neq \widetilde{V}$.

This example suggests that it is impossible to naively apply matrix completion. However, there may be a way to borrow some ideas from matrix completion literature to further improve the current algorithms. We leave this interesting direction to future work.

3. What if we have $\gamma(V)\ll s_{K-1}(V)$?

Thanks for the comment. This will never happen, as a result of Lemma D in the supplementary material.

For your convenience, we re-prove Lemma D here. By definition, $\gamma^2(V)$ is the maximum diagonal entry of $V'V$. Note that

In addition, let $s_k(V)$ be the $k$th singular value of $V$, so that $s_k^2(V)$ is the $k$th largest eigenvalue of $V'V$. It follows that

Combining the above, we have $\gamma(V) \geq \sqrt{\frac{K-1}{K}} s_{K-1}(V)\geq \frac{1}{\sqrt{2}}s_{K-1}(V).$

Thank you for your helpful comments. We are glad that you appreciate our contributions in both algorithm and theory. We feel especially encouraged that you agreed that "theoretical results are novel and non-trivial". As you have pointed out, the analysis exploits the tail bounds of Gaussian random vectors and is quite subtle (especially to get the tight constants).

Below is our point-to-point response. We also encourage you to read our Response-to-All-Reviewers for other clarifications.

Main comments:

1. Figure 1: "It is not clear to me what the illustration in Figure 1 shows? What does 'idea simplex' denote?"

Thanks for the comment. In this revision, we have combined two panels of Figure 1 into a single plot, to make it easier to compare the true vertices and the estimates from SPA and pp-SPA. The "ideal simple" means the ground-truth simplex formed by the true vertices $v_1,v_2,\ldots,v_K$. We have changed the terminology "ideal simplex" to "the simplex formed by the true vertices" to make it clear.

2. Citations to the biased outward bound: "At times a statement is given without citation,... pg 2: 'since the SPA is greedy algorithm, it tends to be biased outward bound.', but it is not clear where this can be seen."

We realized that since we put SPA and pp-SPA in two separate panels in the old Figure 1, the "biased outward bound" was not easy to see. In the updated Figure 1, the two simplexes are put together, so this phenomenon is more obvious to see.

This phenomenon is caused by the way SPA is designed: At each iteration, it selects $\hat{v}_k$ by finding the point furthest to the original. No matter the original is inside or outside the true simplex, this will yield "biased outward bound". We think this phenomenon was already recognized in the robust SPA literature, although not using the exact phrase. In this revision, we added a reference Gillis (2019) in this paragraph.

3. The bound in Theorem 1: "How it shows that the previous results of (Gillis & Vavasis, 2013) is non-satisfactory. It is very difficult to see the comparison with the bound in (5) for $\beta_{new}$."

Your first question is about why the bound in [GV13] is non-satisfactory. We now give a toy example where we fix $(K, d) = (3,3)$, $v_1=(20,20,a)$, $v_2=(20,30,a)$, and $v_3=(30,20,a)$, for a parameter $a\in\mathbb{R}$. Recall their bound and our bound are $g(V) \cdot \beta(X)$ and $g_{\text{new}}(V) \cdot \beta_{new}(X,V)$, respectively, where $\beta_{new}(X,V)\leq \beta(X)$ is easy to see. We now compare $g_{\text{new}}(V)$ and $g(V)$. By direct calculations, we get the following table. It demonstrates that the improvement is substantial.

The main reason is that $s_{K-1}(V)$ is driven by the volume of the simplex, but $s_K(V)$ is driven by the location (center) of the simplex. As $a$ varies, the location of the simplex can change significantly, making $g(V)$ unsatisfactory sometimes.

Your second question is about how to interpret $\beta_{\text{new}}(X,V)$. In this revision, we added another toy example and a new figure, Figure 2, to illustrate this quantity. Intuitively speaking,

• $\beta(X)$ is the maximum noise in all all data points, no matter whether the data points fall inside or outside the true simplex.
• $\beta_{\text{new}}(X, V)$ only cares about the maximum noise of data points outside the true simplex or near the boundary of the true simplex.

As a result, we can add a large noise to any point deeply in the interior of the true simplex without changing $\beta_{\text{new}}(X,V)$, but $\beta(V)$ will be heavily increased. This improvement of $\beta(V)$ to $\beta_{\text{new}}(X,V)$ is critical in the analysis of the new algorithm, pp-SPA. The KNN denoise idea only helps decrease the noise at those points outside the true simplex. If we stick to the original non-asymptotic bound in [GV13], then we will not be able to derive the sub-sequent bounds for pp-SPA.

4. Simulations for larger $d$: "The numerical experiments are done only for very small sizes. I am not sure how big is the effect of denoising using low-dimensional PCA projection in this case."

This is a great point. We have expanded Experiment 1 by considering every $d$ from $1$ to $50$. The effect of low-dimensional PCA projection can be seen by comparing pp-SPA and D-SPA (a simplified version with the projection step skipped). We see that the estimation error of D-SPA deteriorates as $d$ gets larger. However, the estimation error of pp-SPA is almost invariant with respect to the increase of $d$. This shows that the hyperplane projection is crucial for the performance, especially when $d$ is large.

Other questions:

• "How does the theory differ in the case of 3-simplex in 2D? Is there a difference in the $K$ v.s. $K-1$ singular value?"

Thanks for the interesting question. Let's consider two cases: (a) The triangle lives in $\mathbb{R}^2$. Then, $V$ is a $2\times 3$ matrix that has only $2$ singular values. The 3rd singular value is exactly zero. (b) The triangle lives in a hyperplane in $\mathbb{R}^3$. In this case, $V$ is a $3\times 3$ matrix, where 2nd singular value is nonzero as long as the triangle is non-degenerate (i.e., it does not collapse to a line segment). However, the 3rd singular value depends on the location of the hyperplane. When the hyperplane crosses the origin, the 3rd singular value can be zero.

The above conclusions can be deduced from Lemma 3 in the revised draft.

• "How does the algorithm perform for higher dimensional problems?"

Please see our response to your Main Point 4 (in Part 1 of our responses).

• "In the second sentence on the first page you state that Gaussian assumption can be relaxed. Can the theory by applied also when the errors are not Gaussian?"

Yes, the theory can be extended to nonGaussian noise. Take $d=1$ for example. Suppose $\epsilon_i$'s are i.i.d. from a univariate mean-zero distribution $F$. Let $y_k$ be the $k$-th largest value in $\{|\epsilon_1|, \ldots, |\epsilon_n|\}$. In the Gaussian case of $F = N(0,1)$, the key fact we need for our proof is that

To extend our proofs to a general $F$, we need to replace $\sqrt{2 \log(n/k)}$ by another quantity $c_{n,k}$. Such results for a general $F$ can be found in the literature of order statistics. For example, when $F$ is the Laplace distribution, we have

Using similar proofs, we can extend our results to Laplace noise. For $d>1$, we change $|\epsilon_k|$ to $\Vert \epsilon_k\Vert$ in the definition of $y_k$. For settings where either $d$ coordinates of $\epsilon_i$ are independent or $\epsilon_i$ follows an elliptically contour distribution, our proofs can be generalized.

Thank you for the helpful comments. We are glad that you think our problem is well-motivated and that our algorithm is intuitive. Your comment is mainly on the presentation of theorems. In this revision, we have made several changes to make the theoretical results easy to digest. (PS: We also encourage you to read our Response-to-All-Reviewers for other clarifications.)

First, we move Table 1 to the beginning of Section 4. This table provides an explicit comparison of bounds for SPA, P-SPA, D-SPA, and pp-SPA. In the paragraph below Table 1, we also explain carefully how to use this table to see (a) the advantage of hyperplane projection, (b) the advantage of KNN denoise, and (c) the advantage of combining both. Long story short: pp-SPA yields a strict improvement over SPA in all cases.

Second, we have followed your suggestion to explain the orders of $(\alpha_n, b_n, t_n)$ in Theorem 2. Take the case where $K=O(1)$ and $d=o(\sqrt{n})$ for example. The orders of $\alpha_n$ and $b_n$ are directly from Equation (11), and the order of $t_n$ is in the theorem body:

The order of $t_n$ is only used to determine if it is the first or second case of Theorem 2, but $t_n$ itself does not appear in the bound. In addition, it always holds that $b_n=O(\sigma\alpha_n)$ and $\alpha_n=o(1)$. Therefore, $\alpha_n$ and $b_n$ in the bound will be dominated by $\sqrt{\delta_n\log(n)}$, as long as we choose $\delta_n$ that tends to $0$ properly slow. By choosing $\delta_n=\log\log(n)/\log(n)$, we obtain

This is how we get simplified bounds, which are now stated in Equation (12). The above derivation is also in the new paragraph below Theorem 2. At the same time, the bound for SPA is in Equation (14):

It is now easy to compare pp-SPA and SPA directly from Equation (12) and Equation (14).

Finally, we want to explain why we did not present these simplified bounds in Theorem 2 but the more complicated ones. The main reason is that we want to cover as broad settings as possible, including finite or diverging $K$, finite or diverging $d$, and fixed or diverging or diminishing $\sigma$. With the general form in Theorem 2, we are able to get different simplified bounds in different cases. Another reason is that we seek for not only the sharp order but also the tight constant. The "constant" comes from properties of the extreme order statistics of chi-square random variables, which analysis is very subtle. We need to specify $(\alpha_n, b_n, t_n)$ precisely to guarantee that the statements are all rigorous.

Thank you for your positive assessment of our paper and constructive comments. We are very encouraged that you "like the pp-SPA algorithm a lot" and you think it is "an added bonus ... to get tighter theoretical bounds".

Below is a point-to-point response to your comments. We also encourage you to read our Response-to-All-Reviewers.

Main comments:

1. Robust SPA: "Please try to include at least one robust variant of SPA/similar algorithm in the experimental section."

This is a great point. In the revision we have added the comparison with robust-SPA (Gillis, 2019) in all simulations. We used the Matlab code downloaded from $`bit.ly/robustSPA`$. The results are included in the updated Figure 3. We observe that pp-SPA outperforms robust-SPA, especially when $d$ is large.

We also added references to recent variants of SPA, including robust SPA, pre-conditioned SPA, rank-K approximate SPA, etc.. The hyperplane projection and KNN denoise idea have never been used in these algorithms. Furthermore, none of existing works has got as tight bounds as in our work.

2. Writing: "If this is accepted, please spent a few iterations on improving the writing. It's extremely terse at time."

In this revision, we have already re-written some sections (a summary of changes is in our response-to-all-reviewers). We tried to use figures/examples/words to explain the math and notations. Below are a few examples.

• New figure 2 and text therein: They provide an intuitive explanation of $\beta_{\text{new}}(X,V)$ in the non-asymptotic bound.
• The numerical example below Theorem~1: We explicitly calculate the bound in Gillis and Vavasis (2013) and our new bound in a numerical example, to show why the improvement can be substantial in some cases.
• Table 1: This table is a summary of rate comparison and has been moved to the beginning of Section 4.

If the paper is accepted, we are willing to spend more time improving the writing.

3. Whether to include the proof of Lemma 1: "You don't need to prove everything from first principles! For example, Lemma 1 can be given as a HW question for an honors linear algebra class."

We agree that Lemma 1 is a standard result regarding subspace approximation via SVD, which appears in tutorial materials. However, we hope to prove this lemma in the supplementary materials for completion, and we think doing so will increase the readability of our paper. We did not intend to have any credit for such results.

In this revision, we added a sentence in the first line of proof saying that "this is a quite standard result, which can be found at tutorial materials (e.g., https://people.math.wisc.edu/~roch/mmids/roch-mmids-llssvd-6svd.pdf). We include a proof here only for convenience of readers."

4. Statements of Theorems 2-3: "I struggle to see how anyone would cite/use your result without a ton of heavy lifting in understanding your proofs."

We recognized where the confusion came from. In order to make the results as general as possible, we use $\alpha_n$, $b_n$, and $t_n$ in the theorem statements. To get the simple expressions of rates, one still needs to plug in the orders of $(\alpha_n, b_n, t_n)$ from the theorem body. We agree that this is not reader-friendly.

In the revision, we have explicitly explained how to plug in the order of $(\alpha_n, b_n, t_n)$ to get simplified bounds. For example, in the newly added paragraph below Theorem 2, we showed carefully how to use Theorem 2 to get a bound like

In addition, we have included a table in the beginning of Section 4 to summarize all the bounds.

Minor comments:

• Editing of text: We have made all the changes as you suggested.
• Lemma 1: Please see our response to your Main Point 3.
• An explanation of "it is inefficient if we directly apply SPA": We added more explanations on page 2 as follows: "However, since the true simplex vertices $v_1,\ldots,v_K$ lie on a hyperplane of $(K-1)$-dimension, if we directly apply SPA to $X_1, X_2, \ldots, X_n$, the resultant hyperplane formed by the estimated simplex vertices is very likely to deviate from the true hyperplane, due to noise corruption. This will cause inefficiency of SPA."
• If bias outward of SPA is a general phenomenon: This is indeed a general phenomenon (e.g., see Gillis (2019)), and it is caused by the way SPA is designed. At each iteration, the estimated vertex is (by definition) the point in the projected data cloud that is furthest away from the origin. Therefore, the output vertices are likely those points that lie on the furthest positions away from the center of the data cloud. This is why we always have the biased outward bound phenomenon. It is now seen more clearly in the updated Figure 1.