Greedier is Better: Selecting Multiple Neighbors per Iteration for Sparse Subspace Clustering
Authors: Jwo-Yuh Wu
Liang-Chi Huang
Wen Hsuan Li
Chun-Hung Liu
Rung-Hung Gau
Venue: Accepted by TMLR ()

Abstract: Sparse subspace clustering (SSC) using greedy-based neighbor selection, such as orthogonal matching pursuit (OMP), has been known as a popular computationally-efficient alternative to the standard  $\ell_1$-minimization based methods. However, existing stopping rules of OMP to halt neighbor search needs additional offline work to estimate some ground truths, e.g., subspace dimension and/or noise strength. This paper proposes a new SSC scheme using generalized OMP (GOMP), a soup-up of OMP whereby multiple, say $p(\geq1)$, neighbors are identified per iteration to further speed up neighbor acquisition, along with a new stopping rule requiring nothing more than a knowledge of the ambient signal dimension and the number $p$ of identified neighbors in each iteration. Compared to conventional OMP (i.e., $p=1$), the proposed GOMP method involves fewer iterations, thereby enjoying lower algorithmic complexity. Under the semi-random model, analytic performance guarantees are provided. It is shown that, with a high probability, (i) GOMP can retrieve more true neighbors than OMP, consequently yielding higher data clustering accuracy, and (ii) the proposed stopping rule terminates neighbor search once the number of recovered neighbors is close to the subspace dimension. Issues about selecting $p$ for practical implementation are also discussed. Computer simulations using both synthetic and real data are provided to demonstrate the effectiveness of the proposed approach and validate our analytic study.

Published in Transactions on Machine Learning Research (10/2023)
Greedier is Better: Selecting Multiple Neighbors per Iteration for Sparse Subspace Clustering
J Yuh Wu
jywn edu. tw
Institute of Communications Engineering National Yang Ming Chiao Tung University
Liang-Cie Chi Huang
lchuang Erryyau edu. bw
Institute of Communications Engineering National Yang Ming Chiao Tung University
Wen-Hsuan Li
vincent @ryyu. edu.
Institute of Communiccaions Engineering National Yang Ming Chiao Tung University
Chun-Hung Liu
chliu Bece. www.ams.sttt..mmmm edu
Department Electrical and Computer Engineering Mississippi State Universily
Bung-Hung Gau Institute of Communications Engineering National Yang Ming Chino Tung University
gourunghhng edu.
Reviewed OpenReview: /ttps://openreet net/ forum? id= dyD8IbSugm
Abstract
Sparse subspace clustering (SSC) using greedy-basee neighbor selection, such as orthogonal matching pursuit (OMP), has been known as a popular computationalll--ffiie alternative to the standard 4-mminimization based methods. However, existing stopping rules of OMP to halt neighbor search needs additional offline work estimate some ground truths, e.g., subspace dimension and/or noise strength. This paper proposes a SSC scheme using generalized OMP (GOMP), soup-up of OMP whereby multiple, say p(2 1), neighbors are identified per itcration to further speed neighbor acquisition, along with a new stopping rule requiring nothing more than a knowledge of the ambient signal dimension and the number p of identified neighbors in each iteration. Compared to conventional OMP (i.e., p 1), the proposed GOMP method involves fewer iterations, thereby enjoying lower algorithmic complexity. Under the semi-random model, analytic performance guarantees are provided. It is shown that, with a high probability, (i) GOMP retrieve more true neighbors than OMP, consequently yielding higher data clustering accuracy, and (ii) the proposed stopping rule terminates neighbor search once the number of recovered neighbors is close to the subspace dimension. Issues about selecting for practical implementation are also discussed. Computer simulations using both synthetic and real data are provided to demonstrate the effectiveness of the proposed approach and validate our analytic study.
1 Introduction
1.1 Motivation
Subspace clustering (Vidal 2011 Yang et al. 2008 Goh & Vidal 2007) is a key enabling technique in modern unsupervised machine learning and its principles be recapitulated follows. Consider noisy dataset.

---

Published in Transactions on Machine Learning Research (10/2023)
Algorithm 1 SSC-LASSO algorithm
Input: Obscrved datasct y V2..,, data matrix Y YI..yyl, regularization parameter 0
for i = 1 to N do
1) c arg min A|je: || + 1 |ly: Yc, 12 s.1. Gi 0.
2) Normalize of and let C = [c ......... 11 0 C, E ER
end for
3) Set. C= [c,]] - vyl, and G = [sg] where Pijj = |e, + |cg,1-
4) Form an N-node similarity graph in which the edge between nodes i and j has edge weight Di.j-
5) Apply spectral clustering to the similarity graph.
Output:
Partition y = 5i U..
y = (21, yNJ C R" whose ground truth obeys a disjoint union as
y 3....
(1)
where cluster 3'k C IR consists of [R/l > 0 noisy data points coming from a de-dimensional subspace Sk and Dil+... ... + Ve N. A partition of y into the form (L) is widely known as the union-of-subspacees model (Vidal 2011 which underpins panoply of practical data clusters ranging from human face images, hand-written digits, to trajectories of moving objects in videos. Given y with unknown L and da, 1 S k < L, the task of subspace clustering is to uncover the partition (1). Among existing solutions to this problem, sparse subspace clustering (SSC) (Liu et al., 2013 Li et al.
2017 Lu et al. 2019 Elhamifar & Vidal catalyzed by the witnessed success of compressive sensing (CS) and sparse representation (Baraniuk,
2013).
2007 Candès & Wakin 2008 Davenport el al. 2011 Elad 2010), has gained much attention because of its compelling experimental performance and provable performance guarantees. A key ingredient of SSC is identify for cach data point y; neighbor group by using the sparse representation technique. Formally, we collect all data points in y to form the matrix Y = [y1 V2]..NN E RaxN and try to express y, as a linear combination of columns of Y except yı, say, y; = Ye,, where c; [Gi,1 C.2...G,,v]
R- subject to Ci O. SSC aims to find an optimal sparse solution c, so that the columns of Y innexxd by the support of c₁ are correct neighbors of y; (i.c., from the same cluster). Computing of is typically done by solving the following i1miniimzaaiion problem (a.k.a. the Lasso regressor (Hastie et al. 2015)
C = arg min A/|c:l1a + |ly: Ycal12 s.t. C₁2 = 0,
(2)
where 0 < x < 1 is a regularization factor. Once oj is obtained, SSC accordingly constructs a similarity graph with the edge connecting y, and y, with weight Ji, = |c,r] + |cyal (see step 3 of Algorithm 1), followed by spectral clustering (von Luxburg 2007) for final data segmentation (see Algorithms and for outlines of SSC-LASSO and spectral clustering algorithms, respectively). Solving problem (2) is computationally demanding, especially for large-scale high-dimensional datasets. Therefore, low-complexity alternatives using greedy-based neighbor selection, e.g., orthogonal matching pursuit (OMP) (Davenport et al. 2011 Elad 2010), were proposed to perform on par with f-mminiiiaation in many cases (see Algorithm 3 for an outline of OMP-based SSC algorithm). SSC-OMP conducts neighbor identification by computing sequence of orthogonal projections. Specifically, for each y, OMP iteratively identifies one neighbor each time as the data point yields the peak absolute inner product when paired with the residual vector (the initial residual (i) y:). As a new neighbor is identified, its ideutity is added to the "already-delected" neighbor vector ro
index subset Am, and the new residual vector rm is updated as the orthogonal projection (.(0)) = (I (YT YY, Y Y-m)--Y? Am )r/" in which YA.. consists of the columns of Y indexed by Am. Once the neighbor identification process ends, the sparse representation vector oj is then obtained by solving a least-squares problem (see step 4 of Algorithm 3). Notably, existing slopping rules of OMP halt neighbor search if either the number of iterations reaches a predesignated maximum M or the residual power is below a threshold T [Tschannen & Bölcskei, 2018).. While the number M is closoly related to the ground truth subspace dimension, the threshold 7 is determined by the background noise strength. Hence, offine estimating the subspace dimension and noise strength is necessary.
2

---

Published in Transactions on Machine Learning Research (10/2023)
Algorithm 2 Spectral clustering algorithm
Input: Weighted adjacency matrix G = [91.5] E RNXN
1) Let L = I A A-1//GA-1/ where A ding (ai, aN} and a; = Eÿ= Di.-
2) Estimate the number of clusters = mink AR+11 dul. where A is the kth smallest eigenvalue of L.
3) Find VI, V2, VD, the eigenvectors of associated with the I smallest eigenvalues of L, and set NN
V = vi/llvilla /e/lvvall...v///llyy
4) Segment the N rows of V into T. clusters using the K-means algorithm.
5) Declare y; E 51 if the ith row of V is assigned to the /th cluster.
Output: Partition y = 51 ....UU.U U.p
Algorithm SSC-OMP algorithm
Input: Observed dataset y' = (2,, .Ny: data matrix Y yy..yyy,, maximun number of ilerations M, residual vector norm threshold 7
for i = to do
(4) = Yı. Ao = Φ.
Let m = 0, ro
if (m < M) and ([ m [22 > then
(0)
1) m m + 1.
2) Am Am UJ*, where jr arg max (y,, (0) ISitiSN
m-1))-
3) rm = =(I -I-
YI-YAn Y.m (YT, Am YAm) Y 1
end if
4) When the above procodure terminates with M0) itcrations, compute of
arg min cupp CA
lly, Ycl12.
5) Normalize or and let ट. = [c.11- 000-110.,.11 GNNN E E
end for
6) Set C C= [c,]] = [cj ... oyj, and [at]] where Gú,j = |c,,gl + legal-
7) Form an N-point similarity graph in which the edge between nodes and j has edge weight gi.j
8) Apply spectral clustering to the similarity graph.
Output:
Partition y Di U...U
In addition to algorithm development, investigating the mathematical performance guarantees of SSC us- ing fruitful analytical tools from CS also received considerable atlention. The vast. majority of related works focused on investigating sufficient conditions ensuring the so-called subspace detection property (SDP) Soottanokootabi & Candès 2012). that is, neighbor identification correct in the way that the cocfficient c f only il y, and y, are in the same cluster; see (Soltanolkotabi & Candès 2012 Soltanolkotabi et al. 2014
2019) regarding the mmmimiizzation solutions, and (yerr et al. 2013 Heckel & Bölcskei 2015 You et al. 2016b Tschannen & Böleskei 2018 pertaining to greedy search. However, SDP is neither necessary nor sufficient for perfect data segmentation. Indeed, as reported in many studies Ngg et al.) 2001 Vershynin 2018), known type of similarity graph effectuating successful clustering one configured with many intra-cluster and few inter-cluster edges. Evidently, this arises when the sparse regression misidentifies few neighbors, hence violating SDP. The downside of few falsely directed edges from cluster to cluster can be effectively compensated by spectral clustering (see (von) Luxburg] 2007) for more discussions on this issue). Fulfillment of SDP does not guarantee correct cluster- ing, especially when accompanied by meager neighbors. The reason is that, despite no inter-cluster edges, the resultant similarity graph is cut into excessively many isolated pieces; poor graph connectivity in this way tends to cause over-estimation of the number L of clusters, leading to a large data clustering erro (Soltanolkotabi ot al., 2014 Wu et al., 2021). The above facts altogether shod further light on the study of sparse regression for SSC. On the aspect algorithm design, the efforts shall be particularly geared towards
Wu et al. 2021 Wang & Xu 2016; Wang et al.
'nne similarity graph is constructed, e.g., using either & mimiimmmzzaiion or greedy methods, one further employ pruning schemes such Yang et al. (2020]: Qin al. (2233) to obtain an updated graph with improved connectivity.
3

---

Published in Transactions on Machine Learning Research (10/2023)
fast acquisition of plentiful neighbors so that the similarity graph can be fleshed out in right configuration. As to neighbor recovery performance guarantees, much remains to be explored in search of new analysis criteria that are not so stringent as SDP, especially able to reflect neighbor recovery error. In the frame- work of two-step weighted é1-minimzzation (Wu et al.) 2021), the neighbor recovery rate was analyzed, and specifically the probability that the sparse regressor produces at least (₂ 0) correct and at most (2 0) incorrect neighbors was found. Such a probabilistic characterization is intuitive and quite flexible that it directly takes account of the general case with neighbor misidentification. For the special case of ky 0, i.e. error-free neighbor identification, it can reveal how much chance SDP stands with no less than k recovered neighbors.
1.2 Paper Contributions
This paper aims at tackling the aforesaid challenges by revamping the OMP, considering its up-to-par performance, reduced computational complexity, and, most importantly, the inherent flexibility to boost neighbor acquisition. Pivoted on the Generalized OMP (GOMP), prominent variant of OMP that is first introduced in the literature of CS (Wang et al. 2012) and allowed to identify multiple neighbors per iteration, we propose a new sparse regression scheme for SSC, logether with an in-depth analytic study of its neighbor recovery performance guarantee. Specific technical contributions of this paper are summarized as follows.
We propose to employ GOMP an effective alternative to OMP for fast neighbor identification. In particular, we first point out that the deviation (caused by noise) of the residual vector from the desired ground truth subspace, pinned down by the Angle of Deviation (AoD), plays a pivotal role in neighbor identification. According to this fact, we then argue that GOMP, while digging out more neighbors per iteration, enjoys a smaller AoD, which makes it more resilient to noise corruption and able to achieve higher neighbor identification accuracy than OMP.
Efficient stopping rules are crucial for greedy-based neighbor selection. If the algorithm stops carly, we would end up with scant correct neighbors, yet if late, with overly many false ones; either case is apt to cause poor graph connectivity and eventual erroneous data clustering. Alongside the proposed GOMP we devise a new stopping rule geared toward fulfilling the dimension-awarre property, that is, the algorithm is halted once the number of recovered neighbors is fairly close to the subspace dimension (in general this is the right. moment. to leave off neighbor search since the residual thereafter is typically dominated by noise). Mathematically, the proposed stopping rule judges the ratio of residual norms over consecutive two iterations against threshold dependent on the ambient space dimension n, which is known once the dataset is given, and the number of neighbors identified per itcration that is at the designer's disposal. Advantagcously, this dispenses with an extra oflline estimation of the subspace dimension or noise strength as required in the existing
solutions (Soltanolkotabi & Candès
2012
Soltanolkotabi et al. 2014 Wu et al. 2021
Wang & Xu.
2016 Wang et al. 2019 Dyer et al.
2013
Heckel & Bölcskei 2015 You et al., 2016b
Tschannen &
Böleskei 2018).
Capitalized on the work (w et al. 2021 h and under the semi-random model (Soltanolkotabi &
Candès 2012; (Soltanolkotabi et al. 2014), we conduct recovery rate analysis to bear out the claimed merits of the proposed GOMP scheme. Supposing that the ground truth subspaces are well-separated from each other, we first derive an analytic probability lower bound for the event that at least km neighbors (0 < kmm S p) are correct (i.e., from the ground truth subspace) in the mith ileration. Such a local iteration-wise recovery rate result is then exploited to obtain the global recovery rate, namely, the probability lower bound for the event that at least k, correct neighbors in total are identified throughout. The obtained analytic formula shows that, for large data size N and small noise power, GOMP enjoys higher correct neighbor recovery rate than the conventional OMP, thereby confirming GOMP can facilitate fast acquisition of many correct neighbors. Finally, we show that, with a high probability, the proposed stopping rule possesses the dimension-aware property.
To implement the proposed GOMP method, the number p of neighbors identified per iteration should be set beforehand. By further analyzing the obtained recovery rate formulae, the impact
4

---

Published in Transactions on Machine Learning Research (10/2023)
of p on the recovery rate is first discussed. For real-world datasets, oftentimes unable to meet the assumptions required by the semi-random model, recovery rate analysis is rather difficult to carry through. We, therefore, conduct numerical simulations to investigate the selection of aimed at fulfilling the dimension-aware property. Interestingly, both the recovery rate analysis for the semi- random model and our simulation study for real-world datasets indicate that a large p is preferred when (i) the ground truth subspaces are well-separated from each other, or (ii) the data size N is large. Some guidelines for selecting p are suggested accordingly.
1.3
Connection to Previous Works
Efficient: sparse regression for neighbor identification has played a pivotal role in the success of SSC. The standard 41-miniization based method was first introduced in the landmark paper (Elhamifar & Vidal 2013), and many related solutions have been proposed since then. In (Soltanolkotabi et al. 2014 Wu et al.
2021 , iterative re-weighted & 4-miiimizationn was adopted to further improve neighbor identification accuracy
at the expense of higher algorithmic complexity. In YYou et al. 2018 Peng et al., 2013 Matsushima & Brbic 2019), compputationally-effieei solutions under the framework of é E-minimization were then proposed; the
basic idea therein was to pre-process a small amount of data points to acquire side information about the ground truth subspaces, based on which a prumed dataset can then be used to reduce computations. Notably, (You et al., 2016a] Panagakis & Kotropoulos) 2014) utilized mixed-norm regularization to further enhance connectivity of the similarity graph. In contrast to Ei-mininization, greedy search such as OMP is one widely considered low-complexity neighbor identification scheme (Dyer et al. 2013 Heckel & Bolcskei 2015 You ct al. 2016b Tschannen & Bölcskei 2018). Recently in (Chen et al., 2018 Zhu et al. 2019). modified OMP algorithms aiming at improving network connectivity have also been proposed: the methods therein utilized the already established neighbor connections to narrow down the candidate neighbor list, overall promoting neighbor recovery and consequently better network connectivity. It is worth noting that all the existing OMP- based solutions identify one neighbor per iteration, and employ stopping rules calling for a knowledge of the subspace dimension or noise strength. Boosting neighbor recovery via multi-neighhor identiicattion per iteration as well as the development of efficient stopping rules free from aforementioned side information is not yet addressed in the literature of SSC.
Regarding the study of mathematical performance guarantees, sufficient conditions for SDP under the 4 minimization framework have boen investigated in, say, (Elhamifar & Vidal 2013 Soltanolkotabi & Candès 2012) for the noiseless case, and (Soltanolkotabi et al. 2014 Wang & Xu 2016 Wang et al. (019)) for the noisy case. Elhamifar & Vidal (2013) utilized convex geometry techniques to derive sufficient conditions tailored for specialized subspace orientations (e.g., disjoint or independent subspaces), while Solttanloooo & Candès (2012) dealt with the gencralization to subspaces with a non-trivial intersoction. Under noise corruption, Soltanolkotabi et al. (2014) leveraged certain approximation of the LASSO functional and the restricted isometry property of the noisy data matrix to estimate the probability that SDP holds; for LASSO sparse regression, Waan & Xu (2016) further investigated sufficient conditions that the regularization parameter must satisfy in order to guarantee SDP. Recently, Wu et al. (2021) extended the study in (Soltanolkotabi et al. 2014) to provide recovery rate analyses for general neighbor recovery events. For SSC employing greedy neighbor identification, [yer et al. (2013) and You et al. (2016b) considered the noiseless scenario and derived sufficient conditions for SDP using convex geometry analysis; [Tchannenn & Bölcskei (2018) 2013) and (You et al., 2016b) to the Gaussian-noise setting, and derived probability lower bounds for the event the SDP holds. As far as we can soe, all existing studics of performance guarantees for OMP-based SSC revolved around the fulfillment of SDP; the general case when neighbor misidentifiiation occurs is left unaddressed.
then extended the results in (Dyer et al.]
To sum up, while GOMP has been investigated in CS WWng et al., 2012): ) its application and poten- tial impacts on SSC remain yet to be explored. This paper is a first step loward this goal. Thanks to multi-neighhor identification per iteration, the proposed GOMP method boosts neighbor recovery at lower algorithmic complexity as compared to conventional OMP. In addition, our newly developed stopping rule enjoys the dimension-aware property, thereby free from lif-line subspace dimension estimation. Morcover, we leverage recovery rate analysis to derive mathematical performance guarantees for general neighbor re-
5

---

Published in Transactions on Machine Learning Research (10/2023)
covery events. In view of the above achievements, study of SSC with GOMP can contribute to more well-rounded literature on SSC under the framework of greedy neighbor selection.
The rest of this paper is organized as follows. Section 2 first explains why GOMP can outperform conventional OMP, and then introduces the foundations behind the proposed stopping rule, Afterwards, the algorithmic complexity of OMP and the proposed GOMP are analyzed. Section analyzes the recovery rate and discusses the issue of selecting the number p of recovered neighbors per iteration. Section 4 provides numerical simulations to verify our theoretical findings and discussions in Section 3. Section 5 presents the proofs of the main mathematical results. Finally, Section 6 concludes this paper. To ease reading, some detailed technical proofs are relegated to the appendix.
2 Proposed SSC-GOMP
This section introduces the proposed SSC-GOMP scheme. We first brief in Section 2.1 the reason why multiple neighbor recovery in each iteration is favored, in an attempt to motivate our GOMP proposal. Tn Section 2.2, we then encapsulate the foundations behind the proposed stopping rule. Finally, in Section 2.3 we provide algorithmic complexity analysis to justify the computational elficiency of the proposed GOMP as compared to OMP.
2.1 Why Multiple Neighbor Recovery per Iteration?
Recall that OMP iteratively identifies a neighbor each time the data point when paired with the residual vector yields peak absolute inner product (sce step 2 of Algorithm 3). Hence, the orientation of the residual vector in each iteration, in particular, the degree to which it deviates from the ground truth subspace, is important for accurate neighbor identification. To formalize this notion, assume that we are to build a neighbor list for the data point y; coming from the cluster Va whose ground truth subspace is Sy. Impaired by noise, the residual vector (c) computed in the mth iteration (m 1) is perturbed outwards Sy. If mm.ll E Sk and 11 E S (the orthogonal complement of Sk), such
we write !/ (f) perturbation can be pinned down by the angle of deviation (AoD)
in.]l + I mm where ((i)
(f)
mm tan -1 -11= m. 111111111101 m.
(f)
(
(3)
whereby a large on means ( severely deviates from S₁. As the OMP algorithm itcrates, perturbation of the residual would become increasingly severe. This is mainly because, owing to orthogonal projection (see step 3 of Algorithm 1), the current residual rm is obtained by removing from the previous r the component lying in the subspace spanned by the already-selected neighbors, most of which are likely correct. Consequently, the magnitude (6) |l2 of the component mm. S₁ diminishes from iteration to iteration; instead, the term |r/0]_ll which reflects the strength of projected misidentified neighbors (if any) plus noise onto St, is typically non-decreasing with m. Put together, the cascade effect is, therefore, an increase in of with m, rendering the residual ( more and more prone lo neighbor misidentiiication as the algorithm iterates. In this regard, simple remedy for securing enough correct neighbors in few iterations (a. "small m" is favored) is therefore to identify multiple neighbors per itcration, say, the p data points (p > 1) corresponding to the largest absolute inner producls. GOMP is therefore polential solution to meet this goal. Using synthetic data, Fig. (a) clearly demoustrates GOMP yields smaller average AoD than OMP thanks to fewer iteration? this accordingly brings about higher neighbor identification accuracy, as illustrated in Fig. (b) (this issue will be elaborated in Section 3). The reduction in the number of iterations can moreover reduce algorithmic complexity, which is potentially appealing in real-time applications. Indeed, the major computational bottleneck of the OMP algorithm is the orthogonal projection operation (step 3 in Algorithm 3). To recover p(> 1) neighbors, conventional OMP requires p iterations, hence p orthogonal orthogonal projection only, and therefore is
projections. Instead, GOMP calls for just one iteration, more computationally efficient (detailed algorithmic complexity comparison of OMP and GOMP is given at
2 An analytic study of average AoD for GOMP/OMP is rather challenging, and is one of our future works.
6

---

Published in Transactions on Machine Learning Research (10/2023)
80
4
OMM (p- -- 1)
trne neighhor rite
avergee Aoi)
60
--MP (p=3)
40
0.6
20
GGOMP (p
e OMMM
1)
0.2
3 5
7
2
3
5 6
7
8
index of selected neighhor
index of selectod neighbor
Comparison of GOMP and OMP in terms of average AoD and empirical recovery rate. We consider synthetic data set of 135 vectors drawn from L 3 orthogonal subspaces, each of dimension 9, in an ambient domain R100; 45 data points per cluster. The data vectors are sampled uniformly from the intersection of the unil-sphere in I R100 with the ground truth subspace and are corrupted by zero mean Gaussian noise with variance 0,04. For GOMP, p = 3 neighbors are picked per iteration as those when matched to the residual yielding the largest three absolute inner products. A total number of neighbors are recovered using both OMP and GOMP. (a) Left: plot of average AoD upon detection of the kth neighbor, k < 9. For GOMP, every three neighbors are detected in each iteration based on the same residual, leading to a staircase AoD curve. Clearly, GOMP results in smaller AoD thanks to fewer iterations. (b) Right: plot of the truc neighbor rate, i.c., the fraction of true neighbors recovered, versus the index k of the detected neighbor. Benefiting from smaller AoD, GOMP is seen to improve neighbor identification accuracy.
Figure 1:
80
- re-ineerent case, OMP (p 1)
O cohereet case, OMP (p
frue neighhor rate
60
* inherennn case, GOMP (p- 3)
coherent case, GOMP (p=3)
0.8
AoD
40
0.6
-8 -coherent ease, GSMe (p 3) 0
inherenn rase,GOMP (p 3)
20
-echherent easc. OMP (p 1)
0.4
- inherent case, OMP (p- 1)
4
7
3 index & of selected neighbor
6
Index l of selected neighbor
Figure 2: An illustration of AoD of GOMP and OMP in the presence of data coherence. We consider a synthetic dataset of 273 vectors drawn from L = 3 orthogonal subspaces, each of a dimension 9, in an ambient domain 10:00; 91 data points per cluster. To generate coherent data, in each cluster the first 46 data points are uniformly drawn from the intersection of the unit-sphere in R1000 and the ground truth subspace and are corrupted by zero-mean Gaussian noise with variance 0.04, while the remaining 47th~91th points are "small perturbations' of the 2nd~46th points, by additive zero mean Gaussian noise with variance 10-4. The inherent case, wherein all 91 data vectors per cluster are uniformly generated as above, is included as the baseline. We set = 3 for GOMP. (a) Plots of AoD when identifying neighbors for the 1st data point. Since the neighboring points are pair-wise strongly coherent, it is highly likely a coherent pair is picked in each iteration. For both OMP and GOMP, data coherence even reduces AoD as compared to the inherent case. (b) The true neighbor recovery rates are higher, thanks to smaller AoD.
the end of this section). Considering all the above facts, we thus propose to adopt GOMP in place of OMP for neighbor identification.
Finally, we would like to comment on the the case when coherent data points are present. Assume that two correct neighbors of y, are coherent (so that they are aligned toward the same direction). and both are selected by GOMP in the same itcration. As such, the dimension of the subspace identified throughout this iteration is p 1--1 1, rather than p as in the inherent case (i.c., all data points are sufficiently uncorrelated
Remark:
with each other). To obtain the new residual vector (i) the current residual r(() is projected onto the orthogonal complement of the span of the 'already-selected" neighbors, which is of a higher dimension (one more) than the inherent case. In this way, the resultant signal component r(()) m., is better retained, leading to Vann-((1r:)) [[1x=]]) ||1//11=/]] m.] [la) (since the strength [1r(") i lla of the projected misidentified neighbors (if any) plus noise onto S. is roughly the same in both cases). This is confirmed by our experimental study as illustrated in Fig. The above results indicate that GOMP can
a larger (6) la and consequently a smaller AoD = tan
(6)
work well even in the presence of coherent data points.
7

---

Published in Transactions on Machine Learning Research (10/2023)
2.2 Halt When There Are About as Many Neighbors as Subspace Dimension
Since the daia point y; comes from the dæ-dimensional ground truth subspace Sk, a group of around da true neighbors would reach a "critical mass" to well explain y; and, if so, the residual from then on is highly apt to be dominated by noiso, standing very little chance to uncover more true neighbors (this will be born out by our mathematical analysis in Section 3). Grounded on this fact, the algorithm is expected to be halted once dj neighbors or so are available. We shall recall the stopping rule widely adopted in the literature (You et al., 2016b [Tschannen & Bölcskei 2018), which terminates neighbor search when either the number of iterations reaches a pre-set maximum M, or the residual becomes so small that ||r l12 S T for some threshold т > 0. Though implicit, this assumes the availability of prior knowledge about M and т; the former is arguably all about the subspace dimension de and the latter is closely related to the background noise strength, both of which can only be acquired through extra off-line estimation process. Considering that the dedicated overhead of parameter estimation could be costly, below we develop a new stopping rule which is per se aware of the subspace dimension without the need of knowing M or T.
To introduce the proposed approach, let us write the data point under consideration as y; = X; + + where x; is the noiseless signal point and e, is the additive noise. The residual vector r(()) which can be obtained from y, through sequence of m orthogonal projections (step 3 of Algorithm 3), can be expressed as
( m
+ t = II P.y. = II +o) [=1
II P + II 1=1
t=1
where PA is the orthogonal projection onto the orthogonal complement of the subspace spanned by the already-selected neighbors up to the /th iteration. Assume that most of the recovered neighbors up to the m iterations are from the correct subspace S₃ so that the projected signal is very small and the residual (i)
is strongly dominated by the projected noise, i.e., (1) ~ I Pre, This is typically the case once about many neighbors as the subspace dimension are recovered. If the noise e, is Gaussian, so is the residual rm which, being nearly isotropic, tends to distribute its power evenly over all the dimensions (about n - dx) of the orthogonal complement of the subspace spanned by all the already-selected neighbors; that is to say, each dimension shares a factor 1/(n-d3 of the total power ||r! I2- During the (m+ + 1)th iteration, 1 is then obtained from r by removing from it the components along the newly selected neighbors, implying (1)
that ||r//] 112 is close to x ||rm (1 - d&) ~ p x ||r ll//n, /, in which the approximation makes sense since the subspace dimension dy is in general very small in comparison with the ambient space dimension n. Taking the square root and using the triangle inequality, we then obtain the following condition to halt the neighbor search
| ll2
|| [1 11-0-1ll S (leW////
(5)
or equivalently,
|| m++ill 1 Vp/n.
(6)
The proposed halting rule 6 is dimension-aware because, for most cases, it is triggered once the residual rm is dominated by noise owing to the recovery of de neighbors or thereabouts. Notably, the left-hand-side
of G is random variable and, thus, there is no way of ensuring (6) always holds. Instead, under Gaussian noise assumption it is shown in Section 5.3 that, with m = [dz/p], inequality (6 holds with a probability
higher than 1 - 2pe Va/p. As a result, as p/n is close to zero (thus, n/p is very large), the proposed stopping rule 1G can highly likely be triggered [ The proposed halting rule (6) is dimension-aware because, for most
cases, it is triggered once the residual (6) is dominated by noise owing to the recovery of de neighbors or thereabouts. We should moreover note that the
of (6) admits the form of a residual norm
ratio, which advantageously rids off the knowledge of noise strength. Indeed, since r I Pie,, we have
(6)
m [[ 1l2 = am le.lll, for some 0 < am < 1.
2
3Such brobabilistic interpretation of stopping criteria is not uncommon, e.g., the widely considered stopping rule [|r [|2 Tsschannnn & Bökcskej 2018) cannot be guarantced to meet deterministically, as || || is random variable.
8

---

Published in Transactions on Machine Learning Research (10/2023)
Algorithm 4 SSC-GOMP algorithm with the proposed data-dependent stopping criterion
Input: Observed data set y (yI. 12, yN), data matrix Y [yI...N]
for == to N do
Let 0, ro
(E) yi,
1- 2y,, Ao Φ.
(1)
if (1 [m / (1)
m- 2 2 /p/n) then
1) m - m + 1.
2) Am = 1m U Tmm where Tm is the set of cardinality p such that
Ky,, (0) m---l (( m-1)), Vj Tiu 9 Tm.
3) (0) = -I- Y... (Y) Am Y ) ) m-
end if
4) When the above procedure terminates with M6 iterations, compute cr
arg min I - Yell2. caupple))AA m0011
5) Normalize the column vector c be unit-norm. Let c = |c 1., C E R.N.
end for
6) Set C= [- - [oj ... -iv, and G [oc]l where Jü,j lca.9 + loz,l-
7) Form an N-point similarity graph in which the edge between nodes and j has edge weight Di.j.
8) Apply spectral clustering to the similarity graph.
Output: Partition y' 5i U...U
The condition 7 motivates the stopping rule of the form You et al. (2016b);
Tschannen Böleskei (0018)
whereby the threshold - depends on the noise strength |e|ll2- Using 7 the residual norm ratio in the proposed stopping rule 6 reads
m++11
R
0m++1 for some 0 < am+1 < am < 1.
(8)
Qlm
which is clearly independent of the noise strength ||c: The docision threshold in (6) assumes nothing more than a knowledge of the ambient dimension n, which always known in advance, and the parameter p, which is at the designer's disposal. As such, it is free from the need of an extra offline estimation of the subspace dimension or noise strength, in marked contrast with the existing solutions You et al. (2016b); Tschannen & Böloskei 2018). We summarize the proposed SSC-GOMP algorithm in Algorithm 4. of Algorithm 4 that the support of of is A -1; namely, the indexes of 'already-detectted" neighbors up to the (M - 1)th iteration, thereby precluding those identified in the M iteration, i.c., J E Tamo. To see the reason behind this, recall from step 3) of Algorithm 4 that the residual in the M iteration reads
We nole from step 4)
( MO
-M-- YA,
arto (Y Anto YA. YAMO 'YT 1000 --0-1
If the stopping condition
6 is triggered, we have ||r 1 lla Rs [|r(()
(1)
Mi Moyll ll2, which together with implies
YA (YT 200
YA YAMDD) -YT
YYT
( M R0 O.
(10)
- is nearly orthogonal to yj, Tar,) which are Since Agri -1 U Tast,, 100 implies that (1) therefore unlikely to be neighbors y:.
Remark: In case that the projected signal component I Pix; is non-zero, our arguments of deriving 16 remains valid. This is because, as long as sufficiently many true neighbors are recovered, the term Ih P,x, actually acts as a noise; thus, the residual (1) II Pix; + II""; Pie, is indeed "noise-only" irrespective of whether Ih Pix; is zero or not. To see this, let stack the "already-seloceed" neighbors up the mtth
iteration matrix (i) e g=xpm, which admits the form ( x E. where (1) and E the signal point and noise matrices, respectively. In case there are sufficiently many correct neighbors recovered,
have x; (0) for some c. In this way, the projected signal component reads
9

---

Published in Transactions on Machine Learning Research (10/2023)
Table 1: Comparison of GOMP and OMP in terms of running time (x10-5 second).
number of iteration
2
3
4
8
9
OMP (p 1)
3.2
3.4
3.4
3.5
3.6
3.7
3.8
3.9
4.0
GOMP (p 1p=38 3)
3.9
4.1
4.2
X
X
X
X
X
X
m
M
P
P
I
Y
(
(1))
(
X;
m
m
m
1=1
P
71
P
Y(i
(Y
(
(x; - -)
m
(11)
Pi
--
(i)
(Y
-
-
(x: (0) c -m
m
m
m
1=1
Pr
(i) Y
(
( c)
( c). P
[=1
The residual vector accordingly becomes Fm =
-I
P
Em!
c)
+0++1
Ple. +
II
P c e.). which
is again a projected noise. Hence, the proposed stopping rule makes sense and st.ill works even I Pie is non-zero.
2.3 Algorithmic Complexity
We end this section by analyzing the algorithmic complexity of the proposed GOMP and OMP. To ease discussion, assume that pM neighbors are to be recovered. Then OMP requires pM iterations, in which the mth iteration has complexity O(nN + nm + nm3 + m²), 1 S m S p&, respectively (Sturm & Christensen 2012). Instead, GOMP requires just M iterations, the mth itcration with complexity O(nN + npm + n²pm)² (pm))), S m & M, respectively. Notice that, for S m S M, the prath iteration of OMP has the same complexity the mth iteration of GOMP. Hence, it is clear that the OMP scheme calls for additional (p 1)M iterations, among which the complexity is OnnN + nm + nm² + m²) for 1 S m S pM, m # p. 2p.p..... Mp. To further justify the low-complexity advantage of the proposed GOMP, we consider the data set used in Fig. [ and compare the running time of the proposed GOMP (p = 3) and OMP (p = 1); the results are listed in Table [ From the table, we observe the following.
Rumning time increases with m, the number of iteration. This is mainly because, as m increases, the number of the identified neighbors, and hence the dimension the least squares problem, in step 4) increases, leading to longer running time.
Even though the mth iteration of GOMP and the pmth iteration of OMP have the same order of complexity, the running time of the former is slightly higher. The reason behind this is that, while OMP computes the maximal absolute inner product, GOMP seeks the first largest ones: this entails additional sorting efforts and therefore higher running time. Despite this, GOMP involves (p - 1)M less iterations, and overall less running time, than OMP. Based on Table the running time of GOMP is about 37.5% 47.5% of that of OMP.
3 Theoretical Results
In this section, we present the recovery rate analysis for the proposed SSC-GOMP. Under the semi-random model assumption, Section 3.1 first derives analytic recovery rate formulae. Section 3.2 then discusses the guidelines for selecting the number p of recovered neighbors per iteration.
10

---

Published in Transactions on Machine Learning Research (10/2023)
3.1 Recovery Rate Analysis
To formalize malters, the data vector is assumed to follow the standard additive noise model, that is,
y, = + e,, Iiis
(12)
where X; E R is the unit-norm noiseless signal vector and e; E R" is the noise. The analyses below are built on the popular semi-random model (Soltanolkotabi & Candès 2012 Soltanolkotabi et al. 2014 Wu et al.) 2021 Wang & Xu, 2016), i.e., the ground truth subspaces Si.... St in the partition (1) are fixed but otherwise unknown, whereas the data vectors and noise are random. Such a model is widely used in the theoretical study of SSC, thanks to its interpretability and amiability to analysis. Similar to the previous
works (Wu et al. 2021 Wang et al. 2019 Heckel & Böleskei 2015 Tschannen & Böleskei 2018), the following assumptions are made in the sequel.
Assumption 1. For each S i S N, the signal vector x; R is uniformly sampled from B5, where B& _ {x x E IR",, ||x||| 1} Sg is the intersection of the unit sphere with the subspace Sk. Assumption 2. For each S S N, the noise c; E RR, 1 N, are i.i.d. Gaussian random vectors with zero mean and covariance matrix (o² In),,,,, and are independent of the signal vectors x, 's.
O
To gauge the degree which two subspaces are separated away from each other, we recall the affinity between two distinct subspaces Sy and Si that defined to be (Soltanolkotabi & Candès 2012
off(S St)
(13)
/iin(d d:}
where columns of Us (U,, respectively) form an orthonormal basis for Sa (St, respectively). Assumption 3. The subspace affinity satisfies
9/3d,, + o
max k
(8 - 120) di log N S 4 4oo N
(14)
where < 1.
Notably, Assumptions guarantees that different subspaces are well separated from each other; affinity conditions akin to (14p are also needed in many existing studies of performance guarantees for SSC (Wul et al.] 2021; Wang et al. 2019 Heckel & Böleskei 2015 Tschannen & Bölcskei 2018). Without loss of generalily, we assume in the sequel XN St, therefore yN É Vr., and our goal is to identily neighbor group of yx. Under the above three assumptions it can be shown that GOMP stands a high chance to recover many true neighbors in each of the first [dt/p] iterations. More precisely, we have the following theorem.
Thcorem 1. (Iteration- Wise Recovery Rate): Let f2i,k2.....k k..... be sequence of integers satisfying 0 S km S p, for 1 S m & M S [dz//l- Under Assumplions) 1 to the proposed SSC-GOMP oblains at least kim true neighbors at the mth iteration, S m & M, with a probability exceeding
0 de (M
M
2e(N - D/d) km + 1) N//1L
p p-km+1
1 1 Ne-Ne- n/8
m=1
(15)
2 (v.
1)2 Villdde- km (c(12)111) 1) -1)
Rm-1
4 2c
1 > 0).
+
+
π
Km -
N2
where 1()) is the indicator function, c > is a constant, and di is the dimension of the subspace SL.
Proof: See Section 5.1.
Further scrutiny reveals the lower bound 15 high in most practical cases. Indood, with a large ambient dimension n and small noise level o, the second and third terms in (15) are kept small; regarding the last summation, the first and third terms scale like Nil--81og N/dd2/p and N-2, respectively, whereas the
11

---

Published in Transactions on Machine Learning Research (10/2023)
second term decays exponentially fast \// increases, thanks to 0 < < 1 (see Assumption 3 When
specialized to p and km = 1 for all I m < M, T.e., the case with conventional OMP subject to all
detected neighbors being true, the lower bound 15 then reads
de de-P(M11)
2e(N - 1VEE)) N8log N/dz
(v (v2- = π
Dill-d-11
4 2c + N2
N/5
n/8 -6
(16)
-
also reported in (Tschannen & Bölcskei 2018 which addressed sufficient conditions ensuring correct neighbor recovery. Based on the iteration- recovery rate result given in Theorem [1 the following corollary further establishes the global recovery rate, namely, the probability of the event that GOMP succeeds in recovering a specified total number of true neighbors throughout the M < [dz/p] iterations.
Notably, probability lower bound akin to 16
Corollary 1. (Global Recovery Rate): Let 0 S ho S pM be an integer and write k = Mge + Tt, where 0 rt & 1 Under the same selup as in Theorem 0 GOMP can recover at least ki brue neighbors in total throughout M(S [dz/pl) iterations with a probability higher than
de dc-P(M-1)
1 1 N/88/8
VT
Vel-dz dr 90-1
2e(NV - 112)
p-
2 πT
e((DV 4t
4+ 2c
N8 Nidi
+
N2
(17)
De dal-de-4: -
2e(N - Vr
p p-
2 π
[e ee(D'/) 1)
9
4+ 2c
(M rr)
+
1(4
0).
qt 1 N78 N/dx,
qt
N2
Proof: See Section 5.2.
By following our examination of (15 it is easy to check the probability lower bound [17 is high in most practical cases. Morcover, with the aid of (17). GOMP is to yield higher true neighbor recovery rate as compared with the conventional OMP. To better illustrate this, we assume without loss of generality that group of pM neighbors is to be found, while demanding at least k = kM true neighbors, for some 1 < k S p (that is to say, k/p of the total neighbors are true). Hence, GOMP requires M iterations, and OMP p times more. Under this sctting, the lower bound 17 for GOMP reads
de PM+p
2c(N (V2))
Ap-k+1
1 Ne -n/8
n/8
M
VT
(p - k + 18 N/dz
Di -
(18)
2
el ee(DV - 1) k-1
(4 2c)
M
-T
h: -
N2
whereas for OMP the bound becomes
de -PP--1 2kMe(N - (D2))
2
(4 2c)
N-Ne-ao a/8
kM
(19)
N8 log N/ds
N2
Note that 180 differs from 19 in the last four terms, among which it is clear that 6 (a/v\)) -pM-p
(a/ TT) dy. pll+1 for small noise level 0, and (1+26)M (1+2c) because k > 1. Since the 4th term NO
N²
M
2e( De
p-++1
in 18) and the 4th term 28Me(N- in (19) scale like N/dd NoloRN74
(p=K+1)N==xxNWA
and 1/11 N/dx)), respectively, both vanish whenever the data size N is very large. Also, since 7 < (
both the 5th term M (V²r)
Diul-de-R
k-1
in (18) and the 5th term kM (v²7) D'il-d-- in 19
decay exponentially fast to zero the cluster size [rr.] grows, therefore negligible. Accordingly, when the
12

---

Published in Transactions on Machine Learning Research (10/2023)
subspaces are well separated from each other and the data/chuster size is large enough, GOMP improves the true neighbor recovery rate.
Theorem and Corollary [ specify the neighbor recovery rate up to the [dz/p]th iteration, upon which GOMP is expected to be halted by the proposed stopping rule (6) thanks to the availability about du recovered neighbors. Such a conjecture provably true with a high probability, as established the next theorem.
Theorem 2. Under the same setup as in Theorem 1
GOMP accompanied by the proposed stopping rule halts neighbor search till the [dz/p] th iteration with a probability exceeding
de (n-n([de//)--
Ne N1/
6
o
2pe - vn/v
(20)
dr.
2e(NN- 1Drl))
+ (V 2 7) N81og N/d2
D1e 12 (e((ee(())))) - 1)
1))
p-1 4 2c
p p-1
+
N2
Proof: See Section 5.3.
Under Assumptions to Theorem 2 provides an analytic probability lower bound for the event that the proposed stopping rule is aware of the subspace dimension, i..., it terminates the GOMP algorithm when the number of recovered neighbors is close to the subspace dimension dr. We further remark that when the noise power increases slightly further, the dimension of the span of data points in V would exceed dy; accordingly, the algorithm would stop with more than [dz/p] iterations in order to recover more neighbors commensurate with the increased dimension. As the noise power grows higher, the residual is soon dominated by noise so as to trigger the stopping condition ( rendering the algorithm terminated fewer than [dr/p] iterations with scant neighbors. This phenomenon will be seen in our simulation study.
3.2 On Selection of p
The performance of the proposed GOMP algorithm depends on the parameter p, i.c., the number of identified neighbors per iteration. With the aid of the recovery rale analyses for the semi-random model, below we first discuss the selection of p aimed at improving the recovery rate. The results will shed light on the selection of p for practical datasets.
Recall the assertion of Theorem that, with a high probability, the proposed stopping rule ( halts the algorithm once dr neighbors or so are recovered and, so, the GOMP algorithm yields high global recovery rate throughout all iterations according Corollary [ Hence, a natural criterion for selecting p is to maximize both the probability lower bounds (17) and (20 in this way, the proposed algorithm stands a high chance of recovering many true neighbors and only few false neighbors. However, we should note that both probability lower bounds are rather complicated functions in and the ground truth parameters (such as the data size N, subspace dimension dr,, cluster size [)rr,,, noise power and the subspace affinity r). An explicit and tractable rule for choosing toward recovery rate enhancement is therefore very difficult to obtain. Even if such a solution can be found, necessarily depends on the ground truth paramcters, which are nonetheless unknown to the designers. Despite this, by analyzing the lower bounds (17 and (20) we can still summarize certain interesting properties of the recovery rate as the parameter p varies.
For this, let us examine cach individual term in the probability lower bounds 17 and (20 to see how they change with p. Suppose that total number M, = pM of neighbors are recovered, among which at least
<<(> 0) are true. To ease discussions, we assume that he 2 M so that ge = [M/M] > 0, where [+] is the floor function; hence, on average at least 4 true neighbors are found in each iteration. We can first observe from 17 the followings:
a) Clearly, the 3rd term - (T)e in the lower bound level is so small that o < vT.
17 increases with p when the noise
b) Following the discussions below 19.. the 4th and 7th terms of 17 scale like -111- (8 log N/dr.llp-s)
and [I (8 log [(2=)](r- respectively. Since
p- -- (M/ = Lo(ks/MM)),
(21)
13

---

Published in Transactions on Machine Learning Research (10/2023)
which increases with p. we conclude that the 4th and 7th terms of (77 also increases with p.
c) To check the 5th and 8th terms of (17), we first note that (i) the subspace affinity bound T < 1 (see (14p), thereby TV 2/2/ < 1, and (ii) since the number ki of recovered true neighbors never
exceeds the cluster size [V,],, we must have [Dr] > ke > at
/ leading to e((D'z) - 1)/9 >
Hence, it can be readily seen that the 5th term T
12 ---------
(=(>)t)-1)) q1
" and the 8th
(V²r)
D22 dd -a (
9-1
r+) ())) grows as p increases, we therefore conclude that the 5th and gth terms decrease with p.
term -MM-
decrease with gr- Notably, since q: Lp(k//M))
d) Finally, the sum of the 6th and 9th terms of 117 equals -M(422) p.
-M+(1+ which increases with pN2
Combining (a)~(d), we know that the 3rd 4th, 6th, 7th and 9th terms of (17 increase with p, whereas the 5th and 8th terms decrease with p. In particular, if the subspace affinity is very small, the 5th and 8th terms in 17 can be neglected so that the lower bound (17) increases with p. Hence, as long as the ground truth subspaces are well-separated from each other, large p is preferred. On the other hand,
consider the scenario that the data size N and cluster size [ are large; therefore the 4th 5th, 7th, and 8th terms in the lower bound (17 are vanishingly small (sce also the discussions below then reduced to
19 ). The lower bound 17 is
1 Ne 1/8 N- -6
σ VT
dr+
M1 2c)
(22)
p.N²2
Clearly, with Mr, N, dL, and a being fixed, the lower bound 22 increases with p. Hence, can also conclude that, as long as the data size N and cluster size [Drr] are large enough, a large is preferred. Next, we consider the lower bound (20). To ease discussion, we assume that p divides dt. Based on our analyses on the lower bound
17) (see the discussions (a) to (d) above) we note that:
e) The 3rd term -6(a/ P(dd/--(dd -6(o/VR) increases with p when the noise level is so small that < VT.
f) The 4th term 2pe v/p negligible since the ambient dimension n > p.
g) Clearly, the 5th term - du - increases with p. N7dz
h) Since 2/2/ < 1 and e((D' - 1)/p (see discussion (c)), (V++)
1D'cl 1D211-d2-2 )
(e(D>)-1)) p-1
increases with p. obviously [dz/p] decreases with p. Hence, it is hard to tell whether the 6th term - -Tdel (v DD (<(D%)-1) of (2 increases with or not.
i) Finally, the 7th term - de
4+2c N
increases with p.
By (e), (f), (g), and (i), we know that the 3rd 4th 5th, 7th terms of (0) increase with p. Still, if the subspace
affinity T is very small, the 6th term of 200 is negligible, leading to the conclusion that (20) increases with
Hence, a large p is favored for well-separated ground truth subspaces. On the other hand, when N and [
are large, the 5th and 6th terms are vanishingly small (again, see the discussions below (1)) so that (0)
increases with p. As a result, large p is preferred whenever the data size N and cluster size Vry are large. It is worth noting that the inferred tendency of (20) p varies totally agrees with that of (17). Finally, we
should note that the cluster size D/z is unknown beforehand. For balanced datasets, wherein all clusters are about cqually large, DLl is typically in direct proportion to the data size N; in this case, we then would like to increase p for large-scale datasets. We can therefore conclude that, for well-separated subspaces (small affinity) or large-scale balanced datasets (large N and [il)). a large p is preferred (the same tendency also seen in our simulation study using both synthetic and real-world datasets). This is expected since, in the former case, the data point under consideration is surrounedd by many true neighbors, and in the latter case tends to be strongly correlated with many neighbors. Fither way, the residual would quickly diminish in
14

---

Published in Transactions on Machine Learning Research (10/2023)
just few iterations, and one should therefore adopt large p to boost neighbor recovery before the residual is depleted of the component from the desired subspace.
We should note however that our recovery rate analyses are developed under the semi-random model, that is, the signal points are uniform (Assumption [] and noise is Gaussian (Assumption 2), which are seldom met in practice. Hence, for real-world datasets, performance guarantees similar to the semi-random model can hardly be obtained. Without such analytic metrics, a simple and intuitively reasonable way of judging good p would be the fulfilment of the dimension-aware property; this is because the availability of about di neighbors is likely to begat similarity graph with good connectivity, ie., sufficiently many intra-cluster edges and few inter-cluster edges. Based on our simulation study (see Section 4 for more details) it is observed that the dimension-aware property tends to hold for wide range of when datasets enjoy a small subspace affinity, such as the Coil-100 (Nene et al. 1996) and MNIST datasets (LeCun et al., 1998). but seems to fail for datasets with medium-to-large subspace affinities, such as the Devanagari (Acharya et al., 2015) and Extended Yale B (Lee et al. 2005) datasets. In particular, the Devanagari datasct is subject large noise corruption and prefers large p; contrarily, the Extended Yale B is under small noise corruption and prefers a small p. We would therefore conclude that the dimension- property holds when the subspace affinity is small; this is because there would then exist ample intra-cluster neighbors for the data point under consideration, making it highly likely to find about dy correct neighbors no matter how many are recovered in each iteration. The above simulation findings also suggest that the selection of p is highly dependent on the affinity and noise power, which are nonetheless known unless additional offline estimation is further conducted. In the absence of such side information, we would then propose to select a large p when the data size N is large; this is justified by our simulation results and also agrees with our previous claim in the scmi-random model casc. In particular, Extended Yale B, Coil-100, and NCKU human face (Chen & Lien] 2009) datasets (all with thousands of data points or less), while p- 5, 6 for large-scale datasets, e.g. MNIST, Devanagari, and Cifar-10 (Krizhevsky) 2009) datasets (with tens of thousands of data points or more).
would suggest choosing p 2, for small-scale datasets such
Remark: In the context of sparse signal recovery via COMP, the selection of p, in general, depends on the sparsity level, but the development of explicit rules for choosing p still remains a difficult open problem (Fu
ct al. 2022
4 Experimental Results
In this section, numerical simulations based on both synthetic data and real-world data are used to cor- roborate our thcoretical study and illustrate the performance of the proposed GOMP method. Synthetic data generated are similar to (Elhamifar & Vidal 2013 Wu et al. 2021 The ground truth is a union of three subspaces S1, S2 and 83 of R-50 with an equal subspace dimension d1 =d2 = d3 = 6; separation between two distinct subspaccs is gauged by the subspace affinity defined in
13 Noiscless signal points
are uniformly drawn at random from each subspace and are corrupted by Gaussian noise with zero mean
and standard deviation a. The sampling density is defined to be ths - 1D/x//d,, 1 I k < 3. Regarding the neighbor identification performance metrics, we consider the true neighbor rate (TNR) (Wu et al. 2021 that is,
EGNET (( + #0) EiSIJEN I(ct) 0) T 2 {(.i))y E M1 for some 1 151LL) L},
TNR A
(23)
where cf,, is the jth entry of the computed sparse representation vector of using COMP/OMP (see step 4 of Algorithm 3 and Algorithm 4), and the average number of recovered neighbors (ANRN)
ANRN
& Σ || Ioo/,,
(24)
which assesses the ability to boost neighbor acquisition. As in (Soltanolkotabi & Candès 2012
Soltanolkotabi
et al. 2014 Wu et al. 2021 You et al. 2018 Matsushima & Brbic 2019), the global clustering performance is evaluated based on the correct clustering rate (CCR), defined as
CCR A (# of correctly clustered data points)/N.
(25)
15

---

Published in Transactions on Machine Learning Research (10/2023)
0.1.1 5
0.3, d
pm0.5. 5, d
0.85
0.9
TNR
0.65
0.85
0.8
p=
0.8
0.75
0m (0m7
0m0.3 7/m7
0m0.5 (
0.85
0.85
0.9
TNE
0.9
TNE
0.85
0.85
0.8
0.75
3.2
0.4
0.8
10.8
0.2
0.4
0.6
3..8
1
1.2
0.2
0.4
0.8
0.8
1.2
0.1. (
0.3.3
0.5, d
0.85
0.85
0.9
YNIt
0.9
TNT
0.85
0.85
0.8
0.75
0.2
0.4
0.8
0.8
1.2
3.2
0.4
0.8
0.8
1
1.2
0.2
0.4
0.8
Figure 3: TNR versus noise standard deviation for nine different pairs of subspace affinity P and sample
density Φ; assume the subspace dimension is known beforehand and the number of iterations is preset to be M = [6/p] (p = 1,23).
1.1. 5
0.0.3
0,0,5 3
0.99
0.98
*
CCR
0.9
+
0.96
0.86
0.94
0.8
0,2
0.4
0.8
0.8
1,2
0.2
0.4
0.6
0.8
1.2
0.2
0.4
0.6
0.8
1.2
6=0
050.5
CCB
0.885
0.99
OCR
3
0.98
0.9
0.97
000
0.5
CCR
0.995
0.99
p
CCR
0.95
*
0.99
0.98
0.9
0.2
0.4
0.6
0.8
1
1,2
0.2
0.4
0.6
0.8
1.2
0.2
0.4
0.5
0.8
1.2
Figure 4: CCR versus noise standard deviation for nine different pairs of subspace affinity p and sample density Φ; assume the subspace dimension is known beforehand and the number of iterations is preset to be M [6/p] (p = 1,2.3).
4.1 Synthetic Data
We first use synthetic data to test the performance of the proposed method. To ease illustration, we consider the case that (f p for all < k # I < 3, i.e., the three subspaces are equally separated from each other, and that Φ1 d2 Ø3 o, with equal sample density. Associated with 9 different pairs (p. of subspacc affinity and sample density, Fig. 3 Jand g plot the simulated TNR and CCR, respectively, versus noise standard deviation σ; sub-figures on the same row (columm, respectively) correspond to an identical & (p, respectively), while p (Ø, respectively) is increased when going from left to right (top to bottom, respectively). In generating Fig. 3 and the total number iterations preset be M for OMP and M [6/p] for GOMP, respectively: this represents the ideal situation that subspace dimension (equal to 6) is known. It is first seen from Fig. that, as compared to the conventional OMP (p 1), the
16

---

Published in Transactions on Machine Learning Research (10/2023)
TNR
-
proposed 2. progoosed propoeed
TNN
propoeed 2. propooedd propoeed
TNNR
* propoeed p 2. proposed
0.9
0.9
0.9
p
propoeed
0.8
0.8
0.6
0.8
0.8
0.6
0.8
7=7
1=7
5. 1=7
TNN
P 3, proposed paogooed
TNNR
+- p 3, proppeed propueed
TNN
--P 3, propoeed p propueed
0.9
0.9
progosed
proposed
prroposed
0.8
0.8
0.6
0.2
0.4
0.6
0.6
0.2
0.4
0.6
0.8
02
0.4
0.6
0.8
0.3. km9
0.5, dmm
3. progored
proposed proposed
proponed 2, proposed
0.9
0.9
0.9
2. proposed
p
progoeee
progceed
p propoeed
0.7
0.
0.8
08
0.8
Figure 5: TNR versus noise standard deviation a for nine different pairs of subspace affinity p and sample 6 is used.
density o; the proposed stopping rulc
1, 5
0.3 5
05
5
0.8
0.8
0.8
COR
* p 3, peopused
pruposed
+ A-p peoposed
0.6
peoposed
0.8
propased peopssed
0.6
peoposed
0.4
p
peoposed
0.4
peoposed
0.2
0.2
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.0
02
0.4
0.8
0.8
0.
0,3, d=7
0.5. 4=7
0.8
0.8
0.8
0.6
- 3. proposed
0.8
*-P
3. proposed
0.6
- 3, peoposed
p propused
peupised
proposed
0.4
proposed
0.4
peoposed
0.4
peoposed
0.2
022
0.4
06
0.8
0.2
0.4
0.8
0.8
0.2
0.4
0.6
0.8
9.1
0,3
0.5, 4
0.8
0.8
0.8
06
+ 3, proposed
cCi
0.6
3, peopased proposed
0.8
* 3, peoposed
proposed
p= proposed
0.4
0.4
0.4
peoposed
proposed
peoposed
02
02
0.
0.4
0.6
0.8
0.2
0.4
0.8
0.8
02
0.4
0.6
0.8
Figure 6: CCR versus noise standard deviation a for nine different paits of subspace affinity p and sample
density Φ; the proposed stopping rule 5 is used.
proposed GOMP with multiple neighbor selection (p > 1) yields higher TNR. In particular, GOMP with = 3 achieves the highest TNR because it involves the fewest iterations (only two) and, thus, is subject to the least perturbation of the residual vector. Then, it is seen from Fig. 4 that improved TNR in turn leads to higher CCR. The above experiment conducted again by instead using the proposed stopping rule 6 to halt the algorithm, and the results are shown in Fig. Fand Fig. G It can be seen that the proposed GOMP still outperforms OMP. Fig. 7 further plots the ANRN of the two methods, both employing the proposed stopping rule (6) for different noise standard deviation o. We first observe from the figure that, when noise is so small that < 0.1, ANRN of all cases is about six, in support of the assertion in Theorem 2 that as many neighbors as the subspace dimension suffice to well explain the data point under consideration. When a grows up to 0.1, the proposed GOMP (with = 2, tends to retrieve more neighbors. As we have mentioned in the discussions at the end of Section 3, further noise corruption would enlarge the dimension of the span of the data cluster; as a result, the algorithm terminated after more than [dz/p] iterations so to recover more neighbors (than the ground truth subspace dimension) explain the data point. When
17

---

Published in Transactions on Machine Learning Research (10/2023)
0m dim5
0m..3 d
nm0.5 dim5
ANRN
6
JARN
2
20
0.2
0.4
0.6
0.8
0.2
0.4
0.8
0.8
0
0.2
0.4
0.6
0.8
0= f/=7
0=3.3 (1=7
0=0.5 df=7
* P
3
ANEWN
8
ANDA
6
P
ANRN
0.4
0.0.1 9lag
0.3, d
p=0.5, dlmg
ANIIN
ANNY
0
2
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
Figure 7: ANRN versus noise standard deviation for nine different pairs of subspace affinity p and sample
density Φ; the proposed stopping rule is used.
005. d=6
0.5. d=2
0=7.7 6=
TNR
0.9
TNR
0.8
TNR
0.7
p
0.7
0.6
P
0.8
0.5
0.2
0.4
0.8
0.8
1.2
0.2
D.4
0.6
0.8
1.2
0.2
0.4
0.8
0.8
1.2
0.7.
CCR
0.96
0.8
P 6
0.9
0.8
0.35
0.4
0.6
0.2
0.4
0.8
0.8
1 1.2
0.2
0.4
0.6
0.8 1
1.2
0.2
0.4
0.8
0.8
1.2
(a)
(b)
(c)
Figure 8: TNR and CCR versus noise standard deviation for three different pairs of subspace affinity p and sample density o; the subspace dimension (dg 6) is known. We set 6 for GOMP so that only one
iteration is conducted, whereas OMP conducts M = 6 iterations.
a gets even larger, the residual is then severely dominated by noise, only to halt neighbor search earlier, say, in less than dd./p] iterations; this ends up with fewer recovered neighbors and reduced ANRN. Finally, we consider the special case that subspace dimension dę is known, and set = dx(= 6) so that GOMP executes just one iteration. Specifically, we consider the following three affinity-density pairs (p. (0..5,,)))) (0.5,2), (0.7,6), and the total number of iterations for OMP set to be M 6. Fig. g plots the TNR and CCR curves. Fig. (a) clearly demonstrates that GOMP outperforms OMP when the subspace affinity p is small and the sample density is large. The results are expected since, with only one single iteration and p de. GOMP can perform well only when cach data point has at least p correct neighbors nearby, a situation fulfilled when subspaces are separated far away from each other (small p) and there are many neighbors around (large o). However, OMP performs better as the sample density & decreases to 2 so that correct neighbors are scant (soc Fig. (b)), or the subspace affinity p increases to 0.7, resulting in many incorrect neighbors (see Fig. (c)). Hence, while GOMP conducts just one iteration and enjoys a smaller AoD, lack of enough good neighbors outweighs the benefit of small AoD.
18

---

Published in Transactions on Machine Learning Research (10/2023)
(a)
(b)
D
2
S
5 6
7
8 9
9h L
J
(c)
(d)
Figure 9: Samples of four real datasets. (a) The Extended Yale B human face dataset. (b) The Coil-100 object image datasct. (c) The MNIST handwritten digit datasct. (d) The Devanagari handwritten character dataset.
Table 2: Running time (in minutes) of OMP (p 1) and GOMP (p 2, 4, when 12 neighbors are recovered.
Rumning time (min)
p= = 1 p =
=
Extended Yale B
1.51
0.85
0.61
0.48
0.3
Coil-100
2.83
1.45
0.98
0.75
0.51
MNIST
194
100
68
50
37
Devanagari
218
130
91
72
47
4.2 Real-Worke Data
We proceed to evaluate the performance of the proposed method by using the following four real-world datasots (an illustration of their test samples is shown Fig. 9
The Extended Yale B human face dataset (Lee et al. 2005 comprised of photos of 38 people each
with 65 images of 192 x 168 pixels. To reduce the computational cost, we use the dimensionality
reduction technique as (Flhamilar Vidal 2013
Peng al.] 2013) Matsushima & Brbic 2019)
for reducing the image size to 48x42.
The Coil-100 dataset (Nene et al. 1996.. composed of 7200 images (each with size 32x32) of 100 objects (72 images for each object).
The MNIST dataset. (TCCn et al. 28x28; at least 6000 images for each digit.
1998
containing 70000 images of handwritten digits with size
The Devanagari dataset. (Acharya et al. acters and digits. We use the character image part, which consists of 72000 images of 36 characters with size 32x32 (2000 images per character).
2015) containing images of handwritten Devanagari char-
In the first part of the simulations, we compare the proposed GOMP with OMP using the above four real datasets. Assuming that a total number of 12 neighbors are to be recovered, Table 2 Jists the running time for OMP (p 1) and GOMP (p 2,3.4, 6); the experiments are conducted using MATLAB 2016b on a desktop with an AMD 3700X CPU of 3.6GHz and 64 GB RAM. It can be seen that GOMP enjoys shorter running time, confirming the computational efficiency of GOMP. Table 3 then compares the resultant TNR, showing that GOMP yields higher TNR when the same number of neighbors are to be found. We note that, unlike the synthetic data case, the TNR of the two greedy methods is less than 0.5 for Yale-B, Coil-100, and Devanagari datasets, in disagreement with the high-TNR assertion of our recovery rate analyses in Section 3. This because real-world datasets seldom meet the assumptions underlying the semi-random model (Assumptions] and : neither the signal points are uniform nor the noise Gaussian (c.g., real human face data are typically corrupted by sparse outliers (Elhamifar & Vidal 2013) Soltanolkotabi & Candès 2012)) Hence, the developed recovery rate formulae in Section 3 may not fully explain the simulated results when
19

---

Published in Transactions on Machine Learning Research (10/2023)
Table 3: TNR of OMP 11=1) and GOMP (p = 2,3, 4, 6)
when 12 neighbors are recovered.
TNR
p = 1
p
p
p
4
p
Extended Yale B
0.37
0.41
0.46
0.49
0.58
Coil-100
0.20
0.21
0.26
0.33
0.45
MNIST
0.58
0.61
0.64
0.66
0.68
Devanagari
0.2
0.23
0.30
0.39
0.51
Table 4: FDR of OMP (p = 1) and GOMP (p= - 2,3, 4,6)
when 12 neighbors are recovered.
FDR
p 1
p
p
p 4
p
Extended Yale B
0.85
0.86
0.87
0.87
0.90
Coil-100
0.94
0.93
0.93
0.94
0.95
MNIST
0.80
0.81
0.81
0.83
0.84
Devanagari
0.51
0.56
0.60
0.65
0.72
Table 5: CCR of OMP (p = 1) and GOMP (p = 2,3, 4, 6)
when 12 neighbors are recovered.
CCR
p 1
p
p=
p 4
p
Extended Yale B
0.47
0.49
0.56
0.63
0.67
Coil-100
0.29
0.32
0.37
0.42
0.48
MNIST
0.43
0.51
0.63
0.64
0.66
Devanagari
0.26
0.31
0.35
0.37
0.38
real-world datasets are considered. Despite TNR may not be high in certain cases, we would like to remark that the existence of false neighbors does not largely degrade final data clustering accuracy. To see this, we further consider the fcature detection rate (FDR) (Hockel Böleskei 2015), defined to be
1
c c
FDR & N
(26)
where č is the vector containing the entries of cj supported true neighbors. If FDR is high, the weights on the false edges in the similarity graph are small, meaning that the existence of false neighbors is somehow not harmful. Table A presents the computed FDR; it can be observed that both OMP and GOMP achieve high FDR (above 0.7 for most cases). Table then lists the CCR, showing that GOMP achieves higher CCR. thanks to higher TNR.
Finally, we illustrate the ANRN of GOMP when employing the proposed stopping rule (6) and discuss the selection of the parameter p aimed fulfilling the dimension-aware property. For this purpose, knowledge about ground truth subspace dimensions and noise level is needed. To determine the subspace dimensions, we vectorize all data points associated with cach cluster, stack them into a matrix, and calculate the normalized singular values (with respect to the maximal one). For each datasel, we choose the first 9 clusters for singular value computation and the results are plotted in Fig. [0] respectively. Observe from the four figures that the normalized singular values of all datascts exhibit sudden slump from decrease slowly: the curves are seen bend to the floors around 7~9 for Yale-B, 5~77 for Coil-100, 10~15 for MNIST, and 15~20 for Devanagari, respectively. The above singular value index ranges then used as the
0.1, beyond which they
20

---

Published in Transactions on Machine Learning Research (10/2023)
normalized singslar value
Excended Yale B
normalized lized singular value
Coil-100
*- chuster
- chnstrr
chuuster /
clusser
+
-chhstee
0.8
cluster
0.8
cluster + #ccluster
*tchuster 8
0.6
chuster
- chaster
chuster 0
0.6
cluster 3 - *-cluster
+-claster 9
0.4
0.4
0.
0.2
0
0
10
20
30
40
50
10
20
30
40
50
Index (a)
Index
(b)
normalized singuiar valuue
MNIST
normalized singular value
Devanagari
claster chcusster&
*cchuster
+-chuster
0.8
claster chustee
+#-clsteer
0.8
+ cluster
chusee
#-cluster
0.6
claster 3 -# chuster
- +ccuster
0.6
cluster 3
chuster
+ cluster 9
0.4
0.4
0.2
10
20
40
50
Index (e)
Index (d)
Figure 10:
Normalized singular values of data matrices of the first 9 clusters of (a) Extended Yale B dataset, (b) Coil-100 datasct, (c) MNIST datasct, and (d) Devanagari datasct.
Table 6: Estimated subspace dimensions, noise powers, subspace affinity, data size, and ANRN of GOMP with proposed stopping rulc.
Coil-100
MNTST
Devanagari
Extended Yale B
Range of subspace dimension
[6,10]
[10,15]
[15,20]
[7.9]
Estimated noise powers
0.036
0.194
0.199
0.032
Estimated subspace affinity
0.447
0.553
0.628
0.706
Data size
7200
70000
72000
2470
ANRN (p 2)
4.7
8.6
6.2
7.5
ANRN (p 3)
6.3
11.6
7.4
11.2
ANRN (p 4)
8.4
13.5
9.4
14.3
ANRN (p - 6)
12.1
17.7
16.8
19.6
estimated subspace dimensions. Also, the noise power of each dataset is estimated as the averaged sum of the rest nondominant singular values. Table 6 then shows the ranges of estimated subspace dimensions, average subspace alfinity, estimated noise powers, and the computed ANRN of GOMP using the proposed stopping rule (6) with p = 3, 4, We can first observe from the table that, for the Coil-100 and MNIST datasets (both with relatively small affinity, though the latter subject to large noise), their ANRNs' are within the ranges 1.7~12.1 and 8.6~17.7, respectively, which are fairly close to the estimated subspace dimensions 6, 10) and [10, 15). This implies that, for these two datasets, the proposed stopping rule (6) is aware of the subspace dimension. However, for the Devanagari dataset (with medium-to-lane to-large affinity) and Extended Yale B dataset (with large affinity, despite small noise corruption), their ANRNs' deviate largely from the subspace dimensions, and hence the dimension-avaae property seems to fail. We particularly observe that, for the Devanagari dataset, the algorithm halted earlier (in less than [d/p] iterations) with a small p (p- 2,3, 4), resulting fewer recovered neighbors and reduced ANRN. This is mainly because, with large noise corruption, the residual is quickly dominated by noise (hence, less than [d/p] iterations); in this case, large p (say p 6) to soon recover many neighbors is
trigger the stopping rule (6) early
21

---

Published in Transactions on Machine Learning Research (10/2023)
desired. Such tendency is also seen in the synthetic data case (see Fig. 7 Regarding the Extended Yale B dataset, GOMP with a small p (p = 2,3) is instead preferred. This is because, with a large affinity, different clusters are close to each other; if is large, the algorithm is likely to pick neighbors from incorrect clusters and is then misled to find neighbors from some subspace (other than the ground truth onc) with a larger dimension, thereby ending up with higher ANRN. Based on the above simulation findings we conclude that, as long as the subspace affinity small, the proposed stopping rule G is aware of the subspace dimension, allowing for a wide range of p. If the affinity is medium-to-large, a good p highly dependent on the noise power: large (small, respectively) p for large (small, respectively) noise power. Overall, the selection of p is therefore closely related to the subspace affinity and noise power, which are nonetheless unknown unless additional off-line estimation is further conducted. Notably, we can also observe from Table O that the good increases with the data size N; this has also been inferred based on our recovery rate analysis for the saniom model (see the discussions in Section 3.2). Hence a very simple guideline, without the need of knowing the subspace dimensions and noise strength, to choose small p (say, p == 2,3) for small-size datasets (with thousands of data points or less) and large p (p 5,6) for large-scale datasets (with tens of thousands of data points or more).
In the second part of this simulation, we compare GOMP using the proposed stopping rule 'proposed') with the following neighbor identification schemes:
6 (dubbed by
GOMP terminated after [d/p] iterations, in which is the ground-truth subspace dimension de- lermined as the floor of the average of the subspace dimensions within the range in Table [ This serves as the ideal implementation of GOMP that employs the knowledge of the subspace dimension (hereafter coined 'ideal").
OMP (Tschannen & Böleskei
2018
terminated after a pre-determined number M of iterations.
Active OMP (AOMP) (Chen et al.
2018 terminated after a pre-determined number M of iterations; the neighbor dropping probability is set to be 0.8, as used (Chen et al., 2018
Restricted connection OMP (ROMP) [Zu et al. 2019); following the same setting as in Zhu et al, 2019), the number of iterations set to be 3, and the number restricted connections set to be
SSC with conventional fi-minimization (LASSO) technique [Elhamifar & Vidal 2013; the regular- ization parameter À is obtained via an exhaustive search over the interval (0, 0.5), as (Matsushima)
& BBbbi 2019).
Exemplar-based subspace clustoring Youu ct al. 2018); the number of exemplar data points is set as 160, as used in (You et al. 2018), and the regularization parameter À is obtained via an exhaustive search over the interval (0, 0.5), as in (Matsushina & Brbic 2019).
Scalable sparse subspace clustering (S3C) (Peng et al. 2013 following the same setting as in (Peng et al. 2013), we randomly choose 1000 data points as "in-sample data" to find coarse estimates of the ground-truth subspaces; the regularization parameter X is obtained via an exhaustive search over (0, 0.5], as in (Matsushima & Brbie 2019).
Selective sampling- based scalable sparse subspace clustering (S5C) (Matsushima & Brbic 2019 : to carry out representation learning to find representative data points, the batch size is set to be 1, used in (Matsushima & Brbic. 2019), and the size of representation set is chosen to be 10L. Again, the regularization parameter X is obtained via an exhaustive search over (0, 0.5] as in (Matsushima & Brbic 2019).
Elastic net subspace clustering with oracle guided elastic net (ORN) YYou et al. 2016a); the noise trade-off parameter y is set to be 1, and the regularization parameter À is obtained via an exhaustive search over (0, 0.5), as (Matsushima & Brbic 2019).
Note that the last four methods, namely, ESC, S3C, S5C, and ORN, are redueed-complexity variants of the E-minimization scheme (Flhamifar & Vidal 2013 For both OMP and AOMP, we conduct exhaustive
22

---

Published in Transactions on Machine Learning Research (10/2023)
Table 7: Simulation parameters for GOMP, OMP, and AOMP.
Parameters
Extended Yale B
Coil-100
MNIST
Devanagari
GOMP-ideal (p)
3
4
GOMP-proposed (p)
2
3
4
6
OMP (M)
4
4
5
AOMP (M)
12
4
6
6
Table 8: CCR of different clustering algorithms (the blank denotes that the time limit of days exceeded).
CCR
OMP
AOMP
ROMP
SSC
ESC
S3C
S5C
ORN
proposed
ideal
Yale
0.73
0.81
0.77
0.81
0.44
0.33
0.71
0.83
0.84
0.84
Coil 100
0.46
0.35
0.38
0.79
0.23
0.14
0.53
0.76
0.71
0.79
MNIST
0.61
0.64
0.63
0.27
0.22
0.55
0.72
0.64
0.66
Deva. 0.28
0.33
0.28
0.20
0.15
0.33
0.40
0.35
0.38
Table 9: Rumning time (in minutes) of different. clustering algorithms (the blank denotes that the time limit of 7 days is excooded).
(minute)
OMP
AOMP
ROMP
SSC
ESC
S3C
S5C
ORN
proposed
ideal
Yale B
0.96
1.55
0.87
577
140
11
13
0.55
0.43
Coil 100
0.95
2.7
0.82
988
6
138
42
31
0.55
0.5
MNIST
77
56
54
83
109
46
69
52
37
Deva.
96
74
69
X
96
142
54
80
59
47
search over the integer set {1,2, 18) to find the best number M of iterations that achieve the highest CCR; such approach has been adopted in the simulation study of (Matsushima & Brbic 2019); for both GOMP- proposed and GOMP-ideal, the parameter for cach dataset is selected according to the suggested guidelines (see Table 7 for a list of the parameters used for OMP and GOMP algorithms). Table g shows the CCR of all methods. Compared to the competing greedy algorithms OMP, AOMP, and ROMP, GOMP PrOP-proposed (proposed) achieves higher CCR. As against the five 4-minimization based solutions, GOMP-propee proposed (proposed) is next to ORN (You et al.] 2016a). which essentially a LASSO regression modified to further promote connectivity. It is also observed that GOMP-proposed (proposed) performs just slightly inferior to GOMP-ideal (ideal), which requires knowledge of the subspace dimension. This indicates that the proposed stopping rule (6) in conjunction with the suggested is aware of the subspace dimension. Table g then lists the running time of all methods. The results show that COMP-proposed (proposed) is very computationally efficient: it enjoys the lowest running time in most cases. For the MNIST and Devanagari datascts, S5C is faster but is inferior to GOMP-ideal (ideal).
5 Proofs
5.1 Proof of Theorem
We set about the proof by defining the following per-iteration neighbor recovery event:
Em A {at least km true neighbors are obtained in the mth iteration), 1 m S M.
(27)
23

---

Published in Transactions on Machine Learning Research (10/2023)
It suffices to show the event N N== Em holds with probability as high claimed in 15 The following lemmas are needed for deriving Theorem
Lemma 1. (Davenport et al. distributed over the wnil-sphere of R and b E R be random vector independent of a. Then, for 2 0 we have
2011 Ex. 25); (Heckel & Bölcskei, 2015, eq. 59)) Let R be uniformly
Pri b/ e||b|l2]] 2e 2e-me
(28)
Prila < c b 22 - 2/TT..
(29)
Lemma 2. Let U, E Raxdy be a matriz whose columns form an orthonormal basis for subspace St, 1 1 For x; & St, we have
(N)
Pr
(xss
(1000))) mill
m.ll
>/>
max
LL
4 (og((N)
UTUULIE ddidz
2 N88 log N)/dz
(30)
Proof: See Appendix A.
0
Lemma 3. For x; St, we have
(N) Pr (xi
m
log N dr. 1 2cN-3
Xis
(31)
||
(N)
l12
m,
Proof: See Appendix B.
Lemma 4. TTchannen & Bölcskei 2018 Lemma 10) The event n (llevll2 S 30/2] occurs with a prob- ability at least - N//8
If km it follows Pr [ because the number of true neighbors is never negative. We then consider the case that km 0. Since the already chosen data vectors in all the previous m- iterations are orthogonal for all j E Am-1, i the "yet-to-bo-selected" candidate neighbors are Yj's for 21, N 1 Am Given that Am-1 A U U where the disjoint subsets A and A consist of, respectively, the indexes true and false neighbors already selected up to the first - iterations, D2/ - | and N [D2l - -A/--1 data vectors remain in the same and different clusters as yN, respectively. Let the true neighbor candidates Yty, where X1, E S₁ and iú & A be sorted according to 16,,, ( - 2 2 1y ( F 1) and likewise the false neighbor candidates y -
the residual r m- i..., Ky,,
---1) (N)
where X $ Sr. and fa & Af in the way that (y (2 1 2.. Ky 1 vnl
(N)
m-11- Also, let the
noiseless signal vectors xi, E St, To & A be ordered so that [(x, ) 2 ...2. 1(x
11
while those x7, SL, Fr & A A/ in the way (x7, r (12 2 ... 1x7 N-122
(M)
The event Em defined in 27 can be expressed in terms the above ordered statistics. By definition, Em occurs when at least km true neighbors, and so at most p - km false neighbors, are recovered in the mth itcration. That is to say, ytam is chosen as a neighbor but Fpp-Am+1 not, justifying the incquality
(y fy-mmmt
(NN m=11) Ky (22 The converse is obviously true. Hence we rewrite E.m. as
Em (y (N) (22 KY te,, (N)
(32)
The expression 332 involves ordered statistics of the absolute inner products between the residual and noisy data points yj's. whose distribution however quite complicated. To case our lower bound derivation, shall instead seek sufficient conditions, specified by absolute inner products between the residual and
24

---

Published in Transactions on Machine Learning Research (10/2023)
noise-free signal points xy's, l,, for Ky,,
(N) m-11) I (yt - 1) this allows employ the assumed well the ground-truth subspace orientation (Assumption
uniform distribution of x, (Assumption
to facilitate analysis. For this we first note from 12 that, for all S s N -1, we have
(y = (x,, 1m1 + (x,; (M) ) + (ej, 1-)
(33)
Then for all km: it follows
My k. rm
(M2 )) = [x,, m + (x,_ (N) + (e, (M))
(N)
) + (7. -m--)) |
m-1,
[X,,
mm1111//
(N)
m-
1 11-11 1/
max
(ety 1
(34)
(N)
m-1.j
max
1
where (a) holds since Xt, E St. for X1, E St, and (b) follows from the ordering of xz, Hence, at least km absolute inner products Kyr, (N) are no smaller than the term the right-hand-side (RHS) of (34;; in particular, for (y r 1 the kmth largest (yr, 1221, must have
t 1 (X
(N) 1n-111/1
(et (222 1)
max
(35)
A
Similarly, by the ordering of 1x=7, r (N) - 21 for all u 2 km + 1, we have
(y (NN S (x7, p- r 1m-211>11
1/X/ r m
(N)
max
(36)
max 1545N -
Putting 35 and 36 together, the condition 1 em 11: the event Em occurs, once the following inequality is true
(N) S (ytuum mm-1) is guaranteed, and hence
(N) |x7, km 11 mm-111
(N) m-11))
max N - 1
(x fa:
(37)
+
max
Ke r 112
max
(Cr,, r (N) < (x7,.. -
111
Recall that our goal is to show that AM Em occurs with a high probability. Under the semi-random model, we go on estimate the probability of which the inequality (37) holds for all S m < M, and in turn a lower bound for Pr(n Em The basic idea is to "split and then lump": find an upper bound for each individual left- haft-han--side (LHS) term (37) and estimate, one by one, the probability about which the obtained incquality holds, and similarly a lower bound for the RHS term along with an estimated probability of its validity; putting the results altogether gives sufficiett condition for (37 along with the desired probability lower bound.
25

---

Published in Transactions on Machine Learning Research (10/2023)
a) An Upper Bound for the 1st LHS Term of 37
For fixed a, we first note that
Pr[1(x7
(N)
m- 1.11 S a}
= 1 - Pr[l(x [x Am mm-1,1 > a
(N)
(4)
2 1 Prf3I CIN |z1 IA
with
||] = km + 1 s.t. 16 (M) - a, Ak aa 13
(38)
(6)
p km + )
- |2
Pril(x - >>, vKE1))
1
max
1: -imm---Am- #1
SCIN 1911 Wh-1
(e)
ee(V - (ill))
p- km
1
max
Prfl(x a, vk [,
p km +
EppE--Pm Fm
where in (a) the notation [Q] 2 (1.2...., Q.Q in (b) (2) _ ক and (c) holds by Cormen
ee al. 2009. inequality (C.5)). By Lemma 1 and since data points are independent, for | p km
we have
Pr 11/1/ w
max (Ng 12L
|U/
(N)
lz, Vk E E
< (/8log Ni/d2 jo-kmtt Am+1 +1
(39)
Vdid..
Combining (38 and 39 and setting max (lox4 UZ 1
(2) llz /Vdide, it then follows
1 UZllE
1(x
r (N) m-1.1//1
11121 I S max log(NN
|| |
(40)
LF
holds with a probability at least
2e 2e(N-1D'll
km
(p- R,. N1 1n Ny/dr
b) An Upper Bound for the 2nd LHS Term of 37 : We first note that
Pr
(X m- 1, )1 s
(N)
ISKEN -DNNN--i-- 1i 1i
1 Pr
1 1(xje ((1) ( > a}
max
1
2 1 - Pr(3k
[N - 112/ 1A s.t. 1 (N) m-_,1)>> a}
N N-13'x11-4 11-12211-WA--11
(1(x)n) ( 11-11 $ {())) }
1 Pr
(41)
N- 1all-AA
1- - Pr
(N) 1 a} (1(x
k=
N-D-II-IA
2
Σ
(N) pr
المددسي > at.
k=1
By Lemma
3 we have
Pr (x x
|/>/)>>>////
6 Glog N
[ ll2
(N)
< 2cN-3, 1 141 - -IA 1..
(42)
dr.
26

---

Published in Transactions on Machine Learning Research (10/2023)
Combing (41) and [42) and with = 6log N
n-dz || (N) 1l,, we obtain
(x/ r (N)
6 6log (1 ll2
max ISASN 1 1
(43)
n2 dp
holds with probability at least 2cN-?
c) An Upper Bound for the 3rd LHS Term of 37 first reach
Using similar techniques as in deriving 44
N-1
Pr
(e/u 1-
S salla
P > a}.
(44)
ISkEN- IA -
Using 28 and with ////le//lll and b (1)) follows
6log N
Pr Kojs 11.0)))>>>
|| llalls/ll
(45)
Pr
K
ele
(N)
6 6log N
|r (N)
112
- 2N
ler 12"
Under Assumption we have
1x12 |PR yNlle S lyxll2 = || + ewlle 1 ||ew|l2.
(46)
in which R(Y A inequality. Combing and 46 yields
is the column space of YA YAm- and the last inequality holds by the triangle
6 log N (1 + |lw (22)2]]] sealla < 2N-3
Pr (e rm
(47)
Setting a = V bog N (1 (evll2)
max
1ef2ll2, (44 and (47) imply
ISAEN-1VAL
14
max
(N) ke
1 1
6 6log
(1 + l|ew|(2)
max 1
1 ll
(48)
holds with a probability least -N
d) An Upper Bound for the 4th LHS Term of (37 By following the same procedures as from we readily show the inequality
44
6 log N
Ket
(N)
(1 ||eNl22)
C 12-
(49)
JAt
--
c) A Lower Bound for the RHS Term of 37 Next, go to find lower bound for the RHS term
of (37).
Using (29 and with a = Xs, and mm 1.11
(N)
we have
Pr
(N) [Xx,, r mill
m
112
15 D -1 1
(50)
Vrr.
which together with the assumption that X 's are independent yields
2 -1 I1,, Vj 1 - (v
DD' 1D'l-114 1il1-1Ah_,1-R.
Pr (x0, r 11.7)) |
-
Vel
π
(51)
I Ic
[|Z]] - A m
[],
| = 122 - -1A km-
27

---

Published in Transactions on Machine Learning Research (10/2023)
It then follows
Pr (Xtum*
1 (100 1 2 T
Vdr. [n [12]
= 1 1 (xtxm 1 (A) 1 < -
-1 l12]
Vdl..
2 1 Pr 31 C [D/ -1 -1-14 with |2] = -D'z - 1 11 - km + 1
s.t. [(xe, ./1,1)1 < Vdl. ||, l|,, KjEl}
2 1
11'21 -
Aim -1
max m- 1mmm
Pr(1(xes: -m. | < v |
ll2,
Vjelj
(52)
14-11%
IC110%.
(a)
(1322) - 1
2 12 11 km
2 1
max
km
1-1 1z2.- 17.1-144, 144
km
1210"%.11-1144 --
[e([)r.) -
km -1
(V²-)
2
Vel Vell-da-ho
2
-
kim
where (a) holds thanks to 51 Hence, the inequality
(x
1.111
||
(N)
(53)
vdr.
Am-1
127.1- 11---
Km
holds with a probability at least
e((2)e111) Rm-11
With the bounds occurs, when the following inequality
40
43
48
(49 and holds:
53
the inequality in 37 is true, and hence the event Em
1 Un dddd
6 log N n dp
max 4 log(N) LL
| |mm
(N)
1l +
||| |1r l2
(54)
+ 2
6 log N
(1+ + [|e
max IsiSN
/le/lle S I
-11mm-1,112:
Vdz
or equivalently,
LTTFEE dd
/3dr
mm111 112
(N)
2 2v |ex 12)
LL
88n dz) log N
(N)
8n log N
1- ll2
(N)
ISiSN
max
/e:112 5
4 log N
(55)
which is obtained by multiplying the inequality (54
throughout by the factor vdz x l/( (N)
l12).
It can then be concluded that occurs once (40
43
48
49
53 and (55) hold for all 1 S m S M.
Hence, a lower bound for Pr( } can be obtained by finding a lower bound for the probability that holds for all S m S M.
(5)
Still, our approach seeks upper bound for the LHS term in (55), hence a sufficient condition guarantecing
55 and proves this bound holds with high probability. Towards this end, we derive an upper bound
for (N) 1n-,,,11, 1 S m S M, an upper bound for ||e 1|2, S i S N, and a lower bound for 1 lla < m S M; lumping these altogether yiclds the claimed result. To begin with, we write
(N) m- 1, ll2 = IPss (I Y A (YT.... Y YAm_)) -YT "YTT. S |Pss I ||e2|22-
(56)
According to Temma 4 the event ~ (/le:/1a s 30/22 occurs with a probability at least 1 - -n/8/8 Hence, we conclude that the following sets of inequalities
||r m 1.11111
(N)
ll2 S 30/22,15 m < M,
(57)
28

---

Published in Transactions on Machine Learning Research (10/2023)
|e:l12 - 30/ SisN,
(58)
hold at once with probability at least - Ne-n/e Next, we will find a lower bound for ||| m- 112 Ilass. For < m < M, we have
(N) m-
P.S. (I (Y YAm-1) YY -
1:11
(59)
= =Psc Ymm cm-1) = YN.II Y Am-.1.) Cm-1>
where YN.ll - PSLYN, YAm-1.1 A Pss Y Am-i and cm-1 - argmin lly Y Cm-11l
Then lower
Cm-1
bound for ||r (1) m la is obtained as
| (M)
1/--11111 ||y..II - Y A1- 1n-1l12
2 min IYN.I YAm-1,1
(60)
IYN.JI Yam 1,11 (YT
Am-1. Y 11) ) 'YT Am- Vx.jil2
Z IYN.JI - -YM- (X
Au YAM-LA) - Ay 1, VN.ulla = ||
is subspace of dimension dr. p(M 1) and both the distributions of xN and eN rotationally invariant in St, the normalized yx.1// v. is uniformly distributed over St ng"=1. Let Vg [v1 V2...Nd, (MM-)((-1)] [=x((=-P(M-1) be matrix whose columns form orthonormal basis for B; augment V8 by adding extra p(M- - 1) columns that the dy columns of Vg [V1 V2... -P(MM-)----Va.. E Raxds, form orthonormal basis for SL- Clearly, we have /N.I/ I N. la Va, where a = (ai a2...udz] obeys uniform distribution over the unit sphere of Rdd and, in particular, Pa(y www.://l....llll Vgä, where ä [a ag... -MMiddppp(MM))) obeys the distribution
where B- = St. R(Y R(YAM-IL -11 )
Since YN.I XN + ex.ll
specified in (Knokhlov
2006 eq. (7)). Noticing [[P (. = ||Vaãl|2 ||]|2, the following set
of inequalities hold
Pr
Ps
YN.Il
)11 = Pr[ S A}
|y
π (dl - p(M 11)11111/22 1)1/2T((( 1)/2) 1/121/1 /allasA
r(dz/2)
(1 (1-1 ( (1)-M)11)2//2 1)-2)/2 dd
(p(Ms
dp (dl pl 11)/2 / IaalasA
(1 ||||3) |
2n
(61)
di (dd (dL-p(M-1) p(M-1)
1 dā
2
||
di.
(dl ((M-1))) v(de - p(M 1))
2
dr. (dL p(M-1)))2 GAMc-D(M-1),
2n
where (a) holds by [Knokhlow] 2006 eq. (7)), (b) follows from (oucaart & Rauhut
2013 eq. (8.1)) in which
u(r) is the volume of unlt-ball in R', and (c) true since v(r) S for all r E N Smith & Vamanamurthy
1989). Using (60) and 61 with x o / 2/dr, we can obtain
(N) m 1l2
2 di
a flly 15 m mmmml M,
(62)
V
holds with a probability at least 1 6(a/ ((M-)) Since yN.ll = x + ex. and ||x ll = (sce Assumption Dh. triangle inequality gives Ily ||2 2 - ||e,.|ll,, which together with (62 implies
|| (2)
--11111
2 V
o (1-1 , (l2), 15 m < M,
(63)
di
29

---

Published in Transactions on Machine Learning Research (10/2023)
holds with probability high [59- condition ensuring (55 by finding upper bound for the LHH of 55
With the aid of
57 (58) and
663
go to find sufficient
Once
(7
58
and 63
hold, we
can obtain the following set of inequalities
UUTUL| vdi
m.
(N)
3/3d + 12)
max 12
88(n d1) log Mir l12
8n log M1r
1n'mm1111
(N)
ISVS
max
|
117FU1.11
3v3dr. +30)
(64)
S max
1L
vdl
(8 1207) ) V log IV
3d 30)
(b)
max aff(S., Sr.) 121) (80oynnoog log
1st.
log N
where (a) follows from the definition
133 and (b) is truc thanks to Assumption 3 As result, when 17
58 and (63) hold, the inequality in 64 for all 1 < m S. M is guaranteed. By employing the union bound
technique can be concluded that the inequality (55) for all 1 <.m m
holds with a probability at least
6(ar /Tr -PMM)1
Ne-n/8, Finally, by using (40), (43)
48).
453
and (55, the proof is thus
completed again by using the union bound.
5.2 Proof of Corollary
The lower bound 117 simply obtained as the probability of the event that the numbers km 's of recovered
true neighbors throughout all M itcrations yield the maximal lower bound (15). That is say, (17) is the maximum of (15) over all leasible < km p, m < M; more precisely, the minimum objective of the following optimization problem
M
min (kI.....R..)
Jo
J(km)
m=1
(65)
Km
his
2, km E (0, 1, p}, VI S S M,
where
(e - 12
p-hm
J(km)
kim + 11 log N/dd
(66)
(v 27) D/ P-ll-de-Re
(e(1)'/)- - kmm 1
Am-1 4 2c +
1 1 > 0)
N²
Below, we will show that the two-level sequence
(k1,..., km)
qt 1, qt + .... 9t + 1, qt 9t
(67)
fold (M-r)) fold
solves 665 our proof.
consequently leading to 17 We recall the following definition and lemmaa, which are needed
Definition 1. (Pečarié et al. 1992 Let x = ( I2 E R.M and [Vr V2-.UM) T E T be two real vectors whose entries are ordered in the way that z(1] z(2] 2 2 Mi and li1 2 U22 2 2 V(M) respeclively. Then x is said to be majorized by y is
[mm]]
Um
(68)
m=1
m=1
and
Fm
ym
(69)
m=1
m==1
We say f A R / is Sclur-coveer if f(x) fly whenever is majorized by y, x,y E A.
30

---

Published in Transactions on Machine Learning Research (10/2023)
Lemma 5. Pečariié et al. 1992, Theorem 12.25) Let f A RM be permatation invoriant function, that is, f(x) = [(x) for all E A and permulation matrices P E RMMM, whose first portial derivalives erist
in A. Then f is Schur-conver in A if and only if
0f Of dx
(at - x,)
2 0, vx= = (o1 M2 E A
(70)
dzz
holds for all si
Using vector-matris notation, all we have to do is show
T
k
A
+
1, Ge + 1,..,, It + 1, Cr. 21 .......9
(71)
fold
(M re) fold
fold
solves the following optimization problem:
min
Ju Σ J(kTvm)
k= ki...... Rml
(72)
s.l. k T
kf vm E 10, p}. LA m S M,
in which vm 0..00 0....0 m- fold
the mth standard unit vector. Tt suffices to prove
Jn(k") S Jo(k)
(73)
for any k E Dj, the feasible set of 72 based on Lemin 5 For this, we first note that the objective Jo is not differentiable, since the function in (66) involves the indicator function. To rid of this difficulty, we
consider the differentiable surrogate J R > IR for J, constructed according to
J(x) (zx,, Vr 2 0.9 and æ 0,
(74)
and
J((0) = JL).
(75)
Thanks to (74,, the function J thus obtained satisfies J(kim) = (kim) for all km E(O, (0.,,...,,1) p}. and J'(km) J'(mm) for all km 2(1,2,....P) p). The corresponding differentiable surrogate for Jo accordingly given by
M
Jo(k)
Σ T(kTvm) v m=
(76)
Clearly, (k) = Jokk) for all E Dj, due to condition (74;; all the better, Jo Schur-convex in D, (a proof is given in Appendix C). Hence, for k, q E DD such that k is majorized by we have Jo(k) Jo(k) S Jo(q) (g). The inequality 73 is guaranteed once k" is majorized by any feasible k, which is indeed true as shown below. For feasible k = [k kg...kr]" such that kı 2 b2 2 2 kM without loss of generality. By Definition it suffices to show
Σ kn km VI 51885 M.
(77)
Assume otherwise that there exists < S M such that E k Σ kms 1 159 s 1, whercas
ES ktt Em=1 kmm- Then have hn > k, which together with k E (1t, at + 1} implies de kim for
I S M, ending in the following contradiction:
k1 Σ km Σ
Rm
km
Σ
q
km
Σ
K
k+
Σ
k kr.
(78)
31

---

Published in Transactions on Machine Learning Research (10/2023)
5.3 Proof of Theorem 2
Below we first derive equivalent condition for the proposed stopping rule that is more amenable to analysis. Recall the residual r/a)) obtained in the mth iteration the orthogonal projection of the previous residual
(N) onto R(YAm)" ) Since || 1r = |r>> |l2 + ( 1 ||2. obtained from the Pythagorean theorem,
the proposed stopping rule can be rewritten as
(N) Fm
S
(79)
(NV
11/2
m
where
(N)
(N) R(Y
F(N)
m
T'
(80)
is the normalized difference of the residual vectors. Rearranging the inoquality in 79 yiclds the following equivalent stopping condition
| S V2las/2vii
2 p/n.
(81)
Below show that, with high chance, SI does not hold when S [dz/pl and is achieved (hence, GOMP stops) with m = [dr/p] + Formally, we prove the following inequalities
| 2 > V 2/p/n p/n, v1 S m S [dr/pl
(82)
and
M 1 S V 22pp/n p/n
(83)
hold at once with probability as high as claimed by Theorem -
To procoed, let y, be a selected data vector in the mth iteration, 1 m S [dr./p]- As long as
53
57
58 and
63 hold, we have
Ky,, r 12 .)l
2
(x1mm
r mm
(N)
max
|(e,, ( m
16 M
(N)
1.11 ll2
/6log N
|| (AV) 1-1--111
max 15ISN
||e||2|
VdL
(lx( m-1111 1, 12
-11- 1122
6 log N
max I8iS
||:/122
(84)
V dr. (1) mn- mn-1,1 1, 11 + ||r (N)
(N)
n
1,11
(N)
v2(1 ///2)
3 /looog. No 2/n
1|2
/9d2, /4 + 4(1 30 2o/22)dd
in which (a) is true since yj is selected neighbor, (b) follows form 35
(c) holds due to
49
and 53), and (d) is obtained by combining 57 PR(YA.. (2/ (1
58 and 63
By definition 80 we have x(
(N) [12), hence
(N)
1y (N) 2
N) 2 P span spanty, (
r mmm1 [
)|l
(85)
(r,--11)0 112)
Ily || (2(2)) 1|2
which together with
84 implies
(N)
7 (N)
(N 112 2
2
77 11xyl12
(86)
|y5lle|11 mm
(N)
||ysllal|r (N)
m-11l2
32

---

Published in Transactions on Machine Learning Research (10/2023)
With the aid of (86), inequality (82) is guaranteed once n/llys122 V2VD/ni p/n, which is typically true since the ambient dimension n is drastically large.
As for 83 again by definition 80 we write
dz/pl++ = PR(Y Maza (T 1 (dz/v) [dz/p] lla).
/
(N) (N)
(87)
Let {bi,.., b(Ide//1++pp be an orthonormal basis of R(Y R(Y A arranged in a way that (bptl,ng b(pp is an orthonormal basis of R(Y Aaa2/))) Since (N)
[di/p] is orthogonal to R(Y (
(N( (dc//1 = for all p+1s11 +iisisssi ((dd/pp 1)p. Thus, it follows
112 122 5 V2Vp/n p/p/ pp/n}
= Pr [IPRR (PPRCY (
(T [dc///
(N)
///l/111/// (N)
[dr./m] 12)ll2 - 22/p/n p/n
2 Pr n [|P span (r (N)
(Td./21//// (12) S 2/ pt//a-p/n//a -p/n////
(N)
(88)
= =1 - Pr (bu) (Fde/)/FT(///)ll))) 12))l2 > 2
(N)
(N)
k=1
21
Pr
PP span (b. (1 (M)
Ide/// /||
Td/// 112)]12>>>>>>))))))>>))) 22vp/n-ppn
1l))l2
(N)
ppp////
k=
= -1- - Pr (
/de//// |1 [d)//o] 112. bb)l
2Vppn p/n/
1=1
According to our proof of Theorem
the neighbors selected in all [d//p] iterations (i.e., all columns of
Y A4a.//)) are correct when 440).
43 (48).
149
53
57
58 and 63) hold for all < m S [dz/p]- If
so, there then exists č E Fpl pp) such that the "noiseless" signal vector Xy and neighbor matrix X (g/x)) satisfy
XN
XA
(89)
Noting that Yauuaal = Aa E /m where /mlag/ol is the noise matrix, the residual vector at the (dr/p/lh iteration can be written
[d./pp = (I Y Yagazm (Y Mazam Y YAagmm) ) YY A
(N)
Araeia (XN + +
=(1 YMTAAAM Mhrar YA Aranin )XN + Y YA Č+eN)
(90)
=(I (I-YA Khraa Ya Araymm (EATa eN): eN))
which is the projection of the Gaussian random vector Espay/o ē ex onto R(Y ((MM))) Being a linear
combination rotationally invariant vectors, EAja, غ + ex remains so; this implies the projection r [dz/r1
(N)
is also rotationally invariant on (((YArum)) Then using 288 with = -fan/// h/anm//li / / lll and b = b
(N)
(di/m)
< k < p, obtain
Pr
K /d//m)
(N))
/d.//) ll2. baYl v 22/pp/ htcp/n-pn//aaa// VD 20 Va/p
(91)
which together with (88 implics
83
holds with a probability at least 1 2pe V//p. In summary, once the
incqualitics 44
43
48
49
133
57
58
and 633 hold for all 1 < m V [dz/p]. then (82 and 83
hold with a probability least
2pe
n/p
The proof is then completed by employing the union bound.
33

---

Published in Transactions on Machine Learning Research (10/2023)
6 Conclusion
Fast acquisition of many true neighbors underlies the success of SSC in real-world applications. Under the framework of greedy selection, which is pretty suited for low-complexity implementation, we propose a GOMP sparse regression scheme, together with a novel subsppaceddimensionavaar stopping rule, to boost neighbor recovery. Thanks to multiple neighbor identification per itcration, the proposed GOMP involves lewer iterations, thereby enjoying even lower algorithmic complexity than conventional OMP. In addition, the residual vector better stands up to noise corruption, consequently bringing about higher neighbor iden- tification accuracy. Our proposed stopping criterion is appealing in that it depends entirely knowledge of the ambient space dimensionn in marked contrast with the existing solution (Tschannen & Böleskei 2018) that requires extra off-line estimation of either subspace dimension or noise power. Besides algorithm de- velopment, in-depth neighbor recovery rate analyses were conducted to justiff the merits of the proposed GOMP; the obtained analytic results are further validated by computations using both synthetic and real datasets. Overall, our study presents a new greedy-based SSC schemi,, with provable performance guar- antees, that can pave the way for practical applications. Future work will analyze AoD of both GOMP and OMP under the considered semi-random model, particularly show GOMP enjoys smaller average AoD; an analytic study of AoD would further offer certain guidelines determining the optimum number of neighbors to recover in each iteration, for the purpose of improving the current scheme, which identifies a constant number (p 2 1) of neighbors throughout all iterations, in more dynamic environment.
Acknowledgments
This work is sponsored by the National Council of Science and Technology of Taiwan under grants NSTC 111-2221-E-A 11--22--------------0 A49-067-MYY, NSTC 112-2811-E-449-555-MMY2 and MOST 110-221-E-A19-022-MY3 The work of C.-H. Liu was supported in part by the National Science Foundation under Awards CNS-2006453 and ECCS-2210106.
References
S. Acharya, A. K. Pant, and P. K. Gyawali. Deep learning based large scale handwritten devanagari charac- ter recognition. International Conference on Software, Knowledge, Information Management and Appli- calions, pp. 1 1-6, 2015.
R. G. Baraniuk. Compressive sensing. IEEE Signal Processing Magazine, 24(4):118-14) 124, 2007.
E. J. Candès and M. Wakin. An introduction to compressive sampling. IEEE Signal Processing Magazine, 25(2):21 30, 2008,
J. C. Chen and J. J. J. Lien. A view-based statistical system for multi-view face detection and pose estimation. Image and Vision Computing, 27(9):1252 1271, 2009.
Y. Chen, G. Li, and Y. Gu. Active orthogonal matching pursuit for sparse subspace clustering. IEEE Signal Processing Letters, 252)):161-168, 2018.
T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algorillmss. MIT Press, 2009.
M. Davenport, M. Duarte, Y. C. Eldar, and G. Kutyniok. Compressed Sensing: Theory and Applications. Cambridge University Press, 2011.
E. L. Dyer, A. C. Sankaranarayanan, and R. G. Baraniuk. Greedy feature selection for subspace clustering. Journal of Machine Learning Research, 14(1):248722517, 2013.
M. Elad. Sparse and Redundant Representations: From Theory to Applications in Signal and Image Pro- cessing. Springer, 2010.
E. Elhamifar and R. Vidal. Sparse subspace clustering: Algorithms, theory, and applications. IEFE Trans- actions on Pattern Analysis and Machine Intelligence, 35:2765-277 2781, 2013.
34

---

Published in Transactions on Machine Learning Research (10/2023)
S. Foucart and H. Rauhut. A Mathematical Introduction to Compressive Sensing. Springer, 2013.
T. Fu, Z. Zong, and X. Yin. Generalized orthogonal matching pursuit with singular value decomposition. IEEE Geoscience and Remote Sensing Letters, 19:1-5, 2022.
A. Goh and R. Vidal. Segmenting motions of different types by unsupervised manifold clustering. IEEE Conference Computer Vision and Pattern Recognition, 2007.
T. Hastic, R. Tibshirani, and M. Wainwright. Statistical Learning with Sparsity: The Lasso and Generaliza- tions. CRC Press, 2015.
R. Heckel and H. Böleskei. Robust subspace clustering via thresholding. IEEE Transactions Information Theory, 61(11):6320 -344222--44,, 2015.
V. I. Knokhlov. The uniform distribution on a sphere in R Theory of Probability and Its Applications, 50 (3):886-399 2006.
A. Krizhevsky. Learning multiple layers of features from tiny images, technical report, university of toronto. 2009.
Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11))2278-25 2324, 1998.
K. C. Lee, J. Ho, and D. Kriegman. Acquiring linear subspaces for face recognition under variable lighting. IEEE Transoctions on Pattern Analysis and Machine Intelligence, 27(5):684 098666888- 2005.
C. Li, C. You, and R. Vidal. Structured sparse subspace clustering: A joint affinity learning and subspace clustering framework. IEEE Transactions on Image Processing, 26(6):2988 3001, 2017.
G. Liu, S. Yan Z. Lin, J. Sun, Y. Yu, and Y. Ma. Robust recovery of subspace structures by low-rank representation. IEEE Transactions Pattern Analysis and Machine Intelligence, 35(1):171 184, 2013.
C. Lu, Feng, Z. Lin, T. Mci, and S. Yan. Subspace clustering by block diagonal reproscntation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 41(2):487-301, 2019.
J. E. Marsden and M. J. Hoffman. Elementary Classical Analysis. W. H. Freeman Co, 1993.
S. Matsushima and M. Brbic. Selective sampling-basee scalable sparse subspace clustering. Advonces in Newral Information Processing Systems, 2019.
S. A. Nene, S. K. Nayar, and H. Murase. Columbia object image library (coil- 100), technical report cues- 006-96. 1996.
A. Y. Ng, M. 1. Jordan, and Y. Weiss. On spectral clustering: Analysis and an algorithm. Advances in Newral Information Processing Systems, 2001.
Y. Panagakis and C. Kotropoulos. Elastic net subspace clustering applied pop/rock music structure analysis. Pallern Recognilion Lellers, 38:46-5 53, 2014.
X. Peng. L. Zhang, and Z. Yi. Scalable sparse subspace clustering. IEEE Conference on Computer Vision and Pattern Recognition, pp. 430-437, 2013.
J. E. Pečarić, F. Proschan, and Tong, Conver Functions, Partial Orderings, and Statistical Applications, volume 187 of Mathematics in Science and Engineering. Academic Press, 1992.
Y. Qin, X. Zhang, L. Shen, and G. Feng. Maximum block energy guided robust subspace clustering, 1EEE Transactions on Pallern Analysis and Machine Intellgence, 45(2):2652 2659, 2023.
D. J. Smith and M. K. Vamanamurthy. How small a unit ball. Mathematics Magazine, 62 1989.
35

---

Published in Transactions on Machine Learning Research (10/2023)
M. Soltanolkotabi and E. J. Candès. A geometric analysis of subspace clustering with outliers. The Annals of Slatistics,
2238, 2012.
M. Soltanolkotabi, E. Elhamifar, and E. J. Candès. Robust subspace clustering. The Annals of Statistics, 42(2):669-655) 191-699--99,99 2014.
B. L. Sturm and M. G. Christensen. Comparison of orthogonal matching pursuit implementations. European Signal Processing Conference, 2012.
M. Tschannen and H. Böleskei. Noisy subspace clustering via matching pursuits. IEEE Transactions on Information. Theory, 64(6):4081,4104, 2018.
R. Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge, 2018.
Vidal. Subspace clustering. IEEE Signal Processing Magazine, 28(2):52 68, 2011.
U. von Luxburg. tutorial on spectral clustering. Statistics and Computing, 17:395-416, 2007.
J. Wang, S. Kwon, and B. Shim. Generalized orthogonal matching pursuit. IEEE Transactions Signal Processing, 60(12): 6216, 2012.
Y. Wang, Y. X. Wang, and A. Singh. A theoretical analysis of noisy sparse subspace clustering dimensionally-reduee data. IEEE Transactions on Information Theory, 65(2):685-706, 2019.
Y. X. Wang and H. Xu. Noisy sparse subspace clustering, Journal of Machine Learning Research, 17(1): 1-41, 2016.
J. Y. Wu, L. C. Huang, M. H. Yang, and C. H. Liu. Sparse subspace clustering via two-step reweighted i-minimization:: Algorithun and provable neighbor recovery rates. IEEE Transactions on Information Theory, (7)2):12 1263,,2222222200001 2021.
A. Yang, J. Wright, Y. Ma, and S. Sastry. Unsupervised segmentation of natural images via lossy data compression. Computer Vision and Image Underslanding, 110(2):212-22)) 2222, 2008.
J. Yang, K. Wang J. Liang, P. L. Rosin, and M. H. Yang. Subspace clustering via good neighbors. IEEE Transactions on Pattern Analysis and Mochine Intelligence, 42(6))1537-1144 2020.
C. You, D. Robinson C. Li, and R. Vidal. Oracle based active set algorithun for scalable elastic net subspace clustering, IEEE Conference Computer Vision and Pottern Recognition, pp.
2016a.
C. You, D. Robinson, and R. Vidal. Scalable sparse subspace clustering by orthogonal matching pursuit. IEFE Conference Computer Vision and Pattern Recognåtion, pp. 3918-999999999 3927, 2016b.
C. You, Li, D.P. Robinson, and R. Vidal. scalable exemplar-based subspace clustering algorithm for class-im iassiimbalanced data. Proceedings of the European Conference on Computer Vision, pp. 67 83, 2018.
W. Zhu, S. Yang, and Y. Zhu. Restricted connection orthogonal matching pursuit for sparse subspace clustering. IEEE Signal Processing Letters, 26(12):1892-1888 1896, 2019.
A Appendix
The following lemma needed for deriving Lemma 2
Lemma 6. (Soltanolkotabi & Candès 2012
Extracted from the proof of Lemma 7.5) Let a E Rdd be |||| 1} and Rdaxd: Then for a we have
distributed uniformly on {x x E Rgd
| ddz
Pr Ball2
< 2c -2212
(92)
36

---

Published in Transactions on Machine Learning Research (10/2023)
Assume x, St, + L. Then have x; U/b, in which b is uniformly distributed over the unit sphere of implies
d by Assumption Setting C= U and a 4og N/vdz, Lemma 6
Pr |1Txxll > 4og |
2e 8(log /2 2 (6 log N/dz
2
(93)
Vdzdi
N N)/dz
Since m/) E EL, we have m// /|/__111112 = ULa with |z|2 = 1. The Cauchy-Schwariz inequality implies
(N) ( / |2. [2)l = 12"U/xx1 12 |102,
(94)
and therefore
(xx, (N)
/I/ mm.llla > (l)
c
{IUXxIll > 4log(N))
U
(95)
vdrdi
ddrdi
With
93
and 95 it follows immediately
U1U7U
Pr X4, (M))
/mm.://lzz ||r////0]]11 la) 4log(NV)
$ N8 N\/d
(96)
Vdide
Appendix B
Let
D(v)= & (Y1............ yn) : / m. 112111122 =v [v]
(97)
f: R" x... x R" + R be the probability density function of (y1. V2..., vN), and S--1 be the unit-sphere of R Then for z E St n we have
Pr (z, //.1/llb /11
(N)
(N) 112)l > =}=
[f
(8) n8a-1 (2))))) D(v)
/
flw..ve....yyddddddyy..... (98)
(8t )n (r.))))) KnD(()
dv.
kenD(v)
where the last inequality holds by defining
x (e: Sk, n gm-1, |/ella S 30/2] ... fx 1x:+ex Sky ns"-1 //ee||2 S 30//2),
(99)
in which we assume that y: E Vk. hence x, E Sk., for 1 s & N. Since the maximum function is continuous and the composition of continuous functions is also continuous, (N) [ |la is continuous function of (2.., y.). Consequently, since (S ng n fr: (r. z) > c] is closed, the inverse image
of D n 1= (r: (r,,)) 2 c)) also closed (Marsden Hoffman 1993 Notice that K. Kn D(S), n 5--1) n (r (r, =}) is closed and bounded in the Euclidean space, and therefore compact by the Heine-Borts theorem (Marsden & Hoffman 1993). Since ((y1.22,.... is continuous the compact set
IC n D(S ns=-1) n [r: (r, 2 tf), by the extreme value theorem Marsden & Hoffman 1993) there
exists constant 8 such that
82,1,22..., -.00---a yN) 5. Vey1.Y2...., e EEEEENDDDSS nsn-))nnn (r: (r. 2) 2 t)),
(100)
37

---

Published in Transactions on Machine Learning Research (10/2023)
which together with
98
implies
(N)
pr
> )
/
[iK]+ + / F(y1. (Wy
(N)
11
m.
(st (n(nknn))(m((
01/c -Ne/
( non )n/r:
(101)
(0) c/A(S) n8m 1 2ce (0-
(8) n ()(())
where (a) holds by lemma with [AC] here denoting the Lebesgue measure of K,, (b) is true as we define
1(S to be the surface area of n and the constant is defined as c 2 (o K] + NN n 1n-1), and (c) follows from 288 With (101), we reach
Pref((
X mm / m (12)
(N) (N)
> > = / /Prrp (x (N) 11 11ll))) > cx,; v] fx. (v 1
R"
(102)
(0)
Zen (m 2en- (n- (v)dv 2en (n-d2)=212
where (a) holds by 101 in which fx, is the probability density function of xg. Let /6log N/(ognN(nn- dz).
Then 102
gives
Pr Xi, (.7) /1 /1 112)
/ (N)
6log N >
n1 dr S 2cN-3.
(103)
Appendix C
In this appendix, we prove the Schur-convexity of Jn. Notably, straightforward manipulations show
(km - ku) OJo
(
OJo
(kim (J'm) J' J'(kg)) O, Whim hig (1,2....,P)
(104)
akim akg
It then follows
JJo
aJoy
(km - kig aki aku
= (km (T) (km)) J (kig)) (km (J' J'(Kg))
(105)
y km, kn E (1,,,...,,p p).
where (a) holds by
74 and (b) is due to 104 When evaluated around ki E {1, (1,2,...,,))) p} and kg = 0, the
resultant partial derivative reads
(ki - -ky) ( ak... akig!
OJo
oJo
km (Jm((()))) J'(o) (a) kim (J J(())
(106)
ki (a'((km))))))))))))) J'())) 20.
where (a) holds by 74 (b) follows from and (c) is true due to 104 Clearly, the function Jo is
75
permutation invariant, together with 105 and 1006) ensures that Jo satisfies the conditions of Lemma and thus is Schur-convex.
38


Reviews:
[{"review_id": "5V1fXp4a2M", "reviewer": "Reviewer_dwHk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\nThe paper provides an interesting direction of using GOMP for SSC. I think this thought is interesting and novel, and the provided recovery rate analysis is perhaps useful for future research (although I am not an expert in analysis to comment on the usefulness of those results). Further, \n1. The paper is relatively well written and easy to read.\n2. There is sufficient background provided that makes the context of the contributions very clear.\n. The experiments are quite elaborate and show promise against OMP; however the performances do not appear better than other approaches to clustering (although I am not too worried about that).\n\nWeaknesses:\n1. The motivation for using GOMP, while intuitive, is less convincing. The key idea to use GOMP is that the angle of deviation (AoD) between the residual strength in the ground truth subspace and its orthogonal complement increases as more subspaces are found incrementally using OMP. However, when using GOMP, won't this problem persist if one does not know the dimensionality d_k of the  subspace? And in case, if one knows the subspace dimensionality d_k, then why not use p=d_k?\n\n2. Further, when working in a sparse subspace clustering setting, there may be several data points y_i that belong to the subspace S_k. Let's say two points y_i and y_j \\in S_k are very similar that selecting them both in a GOMP iteration does not add any benefit to AoD as such. Won't you need a degree of incoherence between the p selected neighbors in GOMP to make them useful for the subspace selection  This latter point also affects the claimed algorithmic complexity benefits as, in such settings, one will expend more compute for finding the p nearest neighbors while in the worst case ending up with the same number of iterations as in OMP for finding the subspaces in the SSC setting. In this sense, I think there appears to be a disconnect in the paper from considering GOMP as a means of replacing OMP in signal recovery against using GOMP for subspace selection in SSC; in the latter, the matrix Y of data points may contain many redundant data points and using p selected neighbors destroys the intended \"sparsity\" desired in SSC? I would think it may be important to include additional assumptions on Y to account for such issues.\n\n3. A second contribution of this paper is the proposal of a halting rule based on the subspace dimensionality as against the know-how of the number of iterations M or the strength of the noise \\tau used in prior methods. Considering Eq. (6), it is not clear to me why this stopping rule is significant or how can it work in a practical setting? More specifically, I presume p << d_k, and is assumed n-d_k ~ n, then p/n is also close to zero? For example, in the real experiments reported in the paper, p is about 6 and n is in thousands. If so, ||r_{m+1}|| / ||r_m|| is very close to 1? if n is large, this stopping criteria (||r_{m+1}||/||r_m|| \\leq 1 - \\sqrt(p/n)) is very close to 1, and in that case, how can one even use this criteria? With this sensitivity on p/n, how is this criteria better than using a thresholded stopping rule as in \\tau used in previous works? A second issue with (6) is that it is derived by assuming that in (4) the signal component \\Pi P_l x_i is close to zero after a few iterations. In reality, we do not necessarily know when this component is close to zero and it is unclear to me if this stopping criteria makes sense when the signal component is non-zero.  \n\n5. The recovery rate analysis appears quite complicated, and there are very strong assumptions that need to be made to make it comprehensible, such as the noise variance \\sigma is small (14). The experiments also do not adhere to the insights in the recovery rates. While, the paper provides various situations when the recovery rates could be understood, it does not appear to be quite aligned with the experimental results. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper explores the use of generalized orthogonal matching pursuit (GOMP) as a replacement for OMP in computing the sparse coefficients in sparse subspace clustering (SSC) -- the main insight being that the strength of the noise subspaces increases as more subspace components are recovered by OMP incrementally such that selecting many components in each iteration may lead to better recovery. The paper also proposes a new stopping criteria for GOMP iterations that depend only on the ambient signal dimension and the number of neighbors selected at every GOMP iteration. The recovery rate analysis is provided under the semi-random noise model. Experiments are provided on synthetic and real datasets and demonstrate promising results. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see above.\n\nOverall, I think the idea is interesting, however, the claimed contributions appear less rigorous and unconvincing. The theoretical results are too complex to have any takeaway insights. The experiments show promise. Overall, I am slightly aligning towards acceptance, although the paper needs to be revised to account for the problems I listed. ", "broader_impact_concerns": "The broader impact is missing, and as I see it, it is unclear. I do not see any ethical issues with this paper. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "XEiY9u0319", "reviewer": "Reviewer_ZEAg", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "strengths: The proposed algorithm seems novel.\nweakness: The paper is hard to follow with some typos, English and notational issues. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a generalized orthogonal matching pursuit (GOMP) algorithm where multiple neighbors are identified at each iteration together with a new stopping criterion. Compared to the classical OMP, GOMP requires fewer iterations and then reduces computational complexity. However, it is not clear whether GOMP takes more running time for each iteration than OMP and thus the claimed comparison fairness needs to be addressed. Recovery rate under certain probability is also theoretically analyzed for the proposed GOMP. Numerical results on synthetic and realistic data sets are provided. ", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Substantial revision on language and careful check of notional consistency are desired.\n1. In (1), replace the comma by the set union sign.\n2. In p.2, it is not clear about, 'errorless neighbor identification alone ... pessimistic condition.\"\n3. In Algorithm 1, it should be \"m<- m+1\" in 1) and \"m-1\" -> \"m\" in 4). It is not clear how the output is related to step 8). \n4. In p.11, \"MNST\" -> \"MNIST\"?\n5. In p.12, \"gen-erated\" -> \"generated\", \"beginequation\"? \n...", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "DKibrUIRc7", "reviewer": "Reviewer_9QED", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. Compared with OMP, the proposed GOMP costs less time for fast neighbor identification.\n2. The proposed GOMP employs a new stopping rule geared toward fulfilling the dimension-aware property for greedy-based neighbor selection.\n3. The paper also provides analytic performance guarantees for the proposed method.\n\nWeaknesses:\n1. This paper is relatively difficult to understand, especially for those who are not very familiar with SSC. In the abstract, the authors do not introduce what shortcomings exist in existing work and what the motivation is for doing the current work. After introducing the motivation and related work, the author did not list specific objective functions and explain the problems that exist in existing work. Instead, they directly begin to introduce their work, which makes it difficult for readers to understand the advantages of the proposed method. The author should address the above issues to make the paper easier to read.\n2. Many mathematical symbols lack explanation, especially in Algorithm 1, such as $\\Lambda_m$, $\\mathbf{Y}$, $\\mathbf{Y}_{\\Lambda_m}$. The author should provide corresponding explanations when symbols first appear.\n3. It is very good to provide theoretical guarantees. However, in the proof process in the paper, the authors heavily relied on the proofs of others, which greatly increase the difficulty of reading the paper. Usually, it is best to only cite other people's theorems or lemmas, and to write the entire theorem in the paper for readers' convenience, rather than citing the proof process.\n\nOverall, I think that this paper needs significant adjustments in order to reach the level required for publication.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, a novel SSC scheme is presented that utilizes generalized OMP (GOMP), an enhanced version of OMP that identifies multiple neighbors in each iteration. The proposed method also introduces a new stopping rule that only requires knowledge of the signal dimension and the number of neighbors identified per iteration. Compared to conventional OMP, GOMP requires fewer iterations and has lower algorithmic complexity. Furthermore, the proposed stopping rule does not require an offline estimation of the subspace dimension or noise strength. Then, analytic performance guarantees are provided under the semi-random model. Finally, numerical experiments are conducted to verify the clustering performance of the proposed method.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please refer to the weaknesses enumerated above.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Based on the strength of reviewer dwHk's recommendation, I am pleased to accept the paper. The reviewer raised several issues, but the authors had good responses that satisfied both the reviewer and myself.  The paper has some theory and also shows some practical benefits, and also has good discussion of motivation for GOMP over OMP.  The paper and figures are professional.

The other two reviewers, while less extensive in their comments, did not raise any issues or objections.

Before accepting the final version of the paper, I request that the authors do a final copy-edit themselves. For example, the title of Section 3 is "Theorectical results" (rather than "Theoretical"), so please at least run a spell-check.
Decision: Accept (poster)