Paper deep dive
LoomSC: Scalable Deep Subspace Clustering with Projector Factorization and Exact Spectral Reduction
Nairouz Mrabah, Youssef Melki, Mohamed Bouguessa, Riadh Ksantini, Shakeeb Murtaza, Tehseen Zia
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 89%
Last extracted: 10/11/2026, 4:18:12 AM
Summary
The paper introduces LoomSC, a scalable deep subspace clustering framework that addresses the computational bottlenecks of dense self-expression matrices and full-affinity spectral clustering. LoomSC utilizes a Latent Orthogonal Optimization Model that jointly learns latent features and a projector self-representation via two thin factors. This approach preserves orthogonality through alternating Procrustes and least-squares updates, constructs a nonnegative quadratic affinity, and reduces the spectral problem to an eigenproblem dependent only on factor width, achieving linear time and memory complexity. The method outperforms state-of-the-art baselines on image clustering benchmarks and scales to 500,000 samples.
Entities (10)
Relation Signals (7)
LoomSC → uses → Projector Factorization
confidence 95% · LoomSC addresses both bottlenecks through projector factorization and exact spectral reduction.
LoomSC → uses → Exact Spectral Reduction
confidence 95% · LoomSC addresses both bottlenecks through projector factorization and exact spectral reduction.
LoomSC → achieves → Linear-time complexity
confidence 90% · For fixed dimensions and iteration budgets, the complete pipeline has linear time and memory complexity in the number of samples.
LoomSC → outperforms → State-of-the-art Baselines
confidence 90% · LoomSC ranks first or second in all 15 dataset-metric comparisons against 9 state-of-the-art baselines.
LoomSC → scalesto → 500,000 Samples
confidence 90% · Synthetic experiments scale to 500,000 samples while maintaining at least 99.8% accuracy.
SSC → uses → l1 Penalty
confidence 85% · Sparse Subspace Clustering (SSC) [8] promotes sparse coefficients through an l1 penalty.
LRR → uses → Nuclear-Norm Regularization
confidence 85% · Low-Rank Representation (LRR) [10] favors a low-rank coefficient matrix through nuclear-norm regularization.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Dense self-expression matrices and full-affinity spectral clustering limit the scalability of subspace clustering. We introduce the Latent Orthogonal Optimization Model for Subspace Clustering (LoomSC), a framework that addresses both bottlenecks through projector factorization and exact spectral reduction. Motivated by the spectral structure of least-squares regression, LoomSC jointly learns latent features and a projector self-representation through two thin factors. Alternating Procrustes and least-squares updates preserve the sample factor's orthogonality while keeping the coefficient matrix implicit. We construct a nonnegative quadratic affinity that preserves the projector's support. An exact feature map then reduces its normalized spectral problem to an eigenproblem whose dimension depends only on the factor width. Neither the full affinity nor the sample Laplacian needs to be formed. Our analysis quantifies the projector approximation and identifies conditions for subspace preservation and within-subspace connectivity. For fixed dimensions and iteration budgets, the complete pipeline has linear time and memory complexity in the number of samples. Across five image-clustering benchmarks, LoomSC ranks first or second in all 15 dataset-metric comparisons against 9 state-of-the-art baselines. Its mean accuracy exceeds the highest baseline mean by 6.66 percentage points. Synthetic experiments scale to 500,000 samples while maintaining at least 99.8% accuracy.
Tags
Links
- Source: https://arxiv.org/abs/2610.10266v1
- Canonical: https://arxiv.org/abs/2610.10266v1
Trouble viewing inline? Open PDF directly →
Full Text
105,951 characters extracted from source content.
Expand or collapse full text
LoomSC: Scalable Deep Subspace Clustering with Projector Factorization and Exact Spectral Reduction Nairouz Mrabah Youssef Melki Mohamed Bouguessa Riadh Ksantini Shakeeb Murtaza Tehseen Zia †thanks: N. Mrabah and S. Murtaza are with the École de technologie supérieure (ÉTS), Montréal, QC, Canada. Y. Melki and R. Ksantini are with the University of Bahrain (UOB), Kingdom of Bahrain. M. Bouguessa is with the Université du Québec à Montréal (UQAM), Montréal, QC, Canada. T. Zia is with COMSATS University Islamabad, Islamabad, Pakistan. Abstract Dense self-expression matrices and full-affinity spectral clustering limit the scalability of subspace clustering. We introduce the Latent Orthogonal Optimization Model for Subspace Clustering (LoomSC), a framework that addresses both bottlenecks through projector factorization and exact spectral reduction. Motivated by the spectral structure of least-squares regression, LoomSC jointly learns latent features and a projector self-representation through two thin factors. Alternating Procrustes and least-squares updates preserve the sample factor’s orthogonality while keeping the coefficient matrix implicit. We construct a nonnegative quadratic affinity that preserves the projector’s support. An exact feature map then reduces its normalized spectral problem to an eigenproblem whose dimension depends only on the factor width. Neither the full affinity nor the sample Laplacian needs to be formed. Our analysis quantifies the projector approximation and identifies conditions for subspace preservation and within-subspace connectivity. For fixed dimensions and iteration budgets, the complete pipeline has linear time and memory complexity in the number of samples. Across five image-clustering benchmarks, LoomSC ranks first or second in all 15 dataset–metric comparisons against 9 state-of-the-art baselines. Its mean accuracy exceeds the highest baseline mean by 6.66 percentage points. Synthetic experiments scale to 500,000 samples while maintaining at least 99.8% accuracy. I Introduction In several cases, high-dimensional data can be modeled by a generation process from the union of multiple low-dimensional linear subspaces. In this context, subspace clustering is a prominent unsupervised strategy that aims to partition the data into clusters by identifying these subspaces, where each cluster is defined by exactly one low-dimensional subspace. The subspace assumption is widely used in visual analysis and has been successfully exploited in several applications, such as image segmentation [1], motion segmentation [2, 3], and face clustering [3, 4]. For example, face images of a person can be represented by a low-dimensional subspace for a fixed pose and under Lambertian reflectance [5]. The most recent subspace clustering methods [6, 7, 8] follow a two-step strategy. First, these methods learn an affinity matrix that captures the similarity between points within the same subspace. Second, a graph clustering approach, such as spectral clustering, is used to identify the clusters. The affinity matrix is constructed according to the self-expression property [8, 9] of the subspaces. More precisely, each data point can be expressed as a linear combination of points in the same subspace. Given a dataset X∈ℝdx×nX ^d_x× n with n data points of dimensionality dxd_x, the optimization problem associated with this strategy is: minC12‖X−XC‖F2+λθ(C), __C 12 \|X-XC \|_F^2+λ\,θ(C), (1) where C=(cij)∈ℝn×nC=(c_ij) ^n× n is a matrix that captures the self-expression property based on the first term; ‖.‖F \|. \|_F denotes the Frobenius norm, which is defined for any matrix M=(mij)M=(m_ij) by ‖M‖F=(∑i=1∑j=1mij2)12 \|M \|_F=( _i=1 _j=1m_ij^2) 12; the second term, controlled by the function θ, serves as a regularization to avoid trivial solutions (e.g., C=IC=I) and to adjust the structure of the self-representation matrix C; λ∈ℝ+∗λ _+^* is a balancing hyperparameter between the two terms. Based on the self-representation matrix, the affinity graph is constructed by computing its adjacency matrix Aabs=(aijabs)∈ℝn×nA_abs=(a_ij^abs) ^n× n, such that Aabs=|C|+|CT|2A_abs= |C |+ |C^T |2. Different regularizers encode complementary structural preferences. Sparse Subspace Clustering (SSC) [8] promotes sparse coefficients through an ℓ1 _1 penalty. Low-Rank Representation (LRR) [10] favors a low-rank coefficient matrix through nuclear-norm regularization. Multi-Subspace Representation (MSR) [11] combines these penalties, whereas Block Diagonal Representation (BDR) [6] directly encourages a block structure. Least Squares Regression (LSR) [9] uses squared Frobenius regularization and admits a closed-form solution. Its grouping effect encourages correlated samples to receive similar coefficients. These formulations seek to suppress cross-subspace connections while retaining informative relations within each subspace. Both requirements matter because a sparse affinity can separate samples from different subspaces yet fragment samples from the same subspace [12]. The self-expression framework nevertheless introduces two distinct computational bottlenecks. Explicitly storing a dense coefficient matrix requires (n2)O(n^2) memory, even when its estimation admits an efficient solver. Constructing a full affinity retains this quadratic storage cost. The spectral stage adds a second bottleneck. A dense eigendecomposition requires (n3)O(n^3) operations, while partial eigensolvers still perform repeated products with the sample-level affinity. Consequently, accelerating coefficient estimation alone does not ensure a scalable clustering pipeline. Both self-expression and spectral assignment must avoid dense sample-by-sample computations. Several strategies reduce these costs. Sampling methods replace the full dictionary with a smaller set of representatives [13, 14, 15, 16]. Subspace-preserving selection seeks representatives that span the underlying subspaces [17]. Sparse solvers accelerate coefficient estimation through greedy pursuit [18], active-set selection [19], or stochastic dictionaries [20]. Structure-aware methods exploit compact graph models [21] or structured affinity estimation [22]. Neural coefficient predictors [23] provide another route by learning a mapping that replaces repeated regression. Representative-based reduction is particularly attractive because it replaces an all-sample dictionary with m≪nm n representatives. Each sample is then described by m coefficients instead of n. Anchor-graph methods can also obtain the spectral embedding from a compact factor [15, 16]. Their effectiveness depends on how well the representatives capture the subspaces in the supplied feature space. Learning the coefficients improves the graph, but does not by itself adapt the image representation. This motivates a compact model in which both the representation and the directions used for self-expression are learned from the data. The representation itself is important because complex appearance variation and data corruption can weaken the linear subspace structure of raw observations. Kernel methods address nonlinear structure through an implicit feature mapping [24]. Deep subspace clustering [25, 26, 27] instead couples an auto-encoder with self-expression in the latent space. The reconstruction objective preserves information about the inputs. The self-expression objective encourages latent features that conform to a subspace model. Joint optimization therefore adapts the features to clustering rather than imposing self-expression on a fixed representation. Joint learning also makes computational compactness essential. A dense self-expression layer introduces n2n^2 sample-dependent coefficients that must be optimized together with the encoder [25, 26]. Neural parameterizations and mini-batch training can reduce training costs [28], but do not by themselves eliminate full-affinity spectral inference. Our goal is to retain both pairwise self-expression and normalized spectral clustering while learning nonlinear representations, without constructing dense sample-level matrices at either stage. To this end, we introduce the Latent Orthogonal Optimization Model for Subspace Clustering (LoomSC). Motivated by the spectral structure of LSR, LoomSC represents self-expression through a projector C=PP⊤C=P . The sample factor P∈ℝn×mP ^n× m has orthonormal columns and a prescribed width m. An equivalent two-factor objective separates latent directions from sample coefficients and couples both to the auto-encoder. Alternating Procrustes and least-squares updates refine these factors while preserving orthogonality. Neither the coefficient matrix nor a sample Gram matrix needs to be formed. The latent directions are updated using all samples rather than remaining fixed to selected representatives. Our analysis quantifies the projector approximation to LSR and identifies conditions under which the full row-space projector preserves subspace memberships and within-subspace connectivity. LoomSC extends this compact representation to spectral assignment through a nonnegative quadratic affinity. Squaring the projector coefficients preserves their zero pattern and thus retains the connections encoded by self-expression. An exact feature map represents this affinity using q=m(m+1)/2q=m(m+1)/2 features per sample. Its normalized spectral problem then reduces to a q-dimensional eigenproblem. This reduction avoids forming either the full affinity or the sample Laplacian. For fixed dimensions and iteration budgets, the complete pipeline has linear time and memory complexity in n. The same compact representation therefore supports feature learning, self-expression, and global spectral assignment. Contributions. (1) We introduce a two-factor projector model for joint feature and self-expression learning with orthogonality-preserving updates. (2) We construct a support-preserving nonnegative affinity and derive an exact reduction of its normalized spectral problem without forming full pairwise matrices. (3) We quantify the projector approximation to LSR and establish conditions for subspace preservation and connectivity. We also prove linear time and memory bounds for the complete pipeline under fixed dimensions and iteration budgets. (4) Across five image-clustering benchmarks, LoomSC ranks first or second in all 1515 dataset–metric comparisons against nine baselines. Its mean accuracy exceeds the highest baseline mean by 6.666.66 percentage points. Experiments on synthetic data reach 500,000500,000 samples while maintaining at least 99.8%99.8\% accuracy throughout the tested range. I Related work In line with the focus of this work, we discuss the advantages and limitations of two existing subspace clustering categories: the scalable subspace clustering category and the deep subspace clustering category. In that respect, we highlight the pertinence of introducing a subspace clustering approach that combines the advantages of both categories. I-A Scalable subspace clustering: Least-squares regression (LSR) [9] provides an efficient alternative to sparse coefficient estimation. Its squared Frobenius regularization encourages a grouping effect among correlated samples and admits a closed-form solution. However, solving the regression analytically does not eliminate the dense n×n× n coefficient matrix, where n is the number of samples. Its explicit construction and storage remain quadratic in n. Thus, improving the regression solver alone does not resolve the full-affinity bottleneck. LoomSC instead adopts a prescribed-width projector model and learns its thin factors directly. This changes the representation of self-expression rather than merely accelerating a dense coefficient solve. Sparse methods reduce storage by retaining a small number of coefficients per sample. SSC-OMP [18] replaces ℓ1 _1 optimization with orthogonal matching pursuit and establishes subspace-preserving guarantees under suitable conditions. EnSC [19] combines ℓ1 _1 and squared ℓ2 _2 regularization to balance subspace preservation and graph connectivity. Its oracle-guided active-set solver restricts optimization to candidate dictionary elements. Both approaches substantially accelerate coefficient estimation. Nevertheless, their dictionary searches still compare target residuals with the full dataset. Repeating these searches across samples retains a quadratic sample-size dependence for fixed sparsity and search budgets. Sparse spectral solvers reduce the subsequent graph-processing cost, but do not remove this coefficient-estimation cost. LoomSC avoids these separate full-dictionary searches by optimizing a shared factorization through thin matrix products and an orthogonal Procrustes update. Stochastic sparse subspace clustering (SSSC) [20] addresses weak within-subspace connectivity through dictionary dropout. For normalized samples, the expected dropout objective induces an additional squared ℓ2 _2 penalty, which encourages denser representations. Its scalable implementation combines coefficients from multiple reduced dictionaries through a consensus formulation. This improves connectivity while preserving a sparse affinity. However, with a fixed retained fraction, each dictionary still grows with n, and the reduced regression must be repeated for every target. Parallel execution accelerates these subproblems without changing their aggregate arithmetic complexity. LoomSC obtains compactness through a common learned latent factor rather than repeated stochastic dictionary reduction. The resulting sample relations are represented jointly without an ensemble of coefficient estimates. Representative-based approaches reduce the dictionary size more directly. LMVSC [15] learns nonnegative sample-to-anchor coefficients and computes its spectral embedding from a compact normalized factor. SGL [16] further couples anchor-graph learning with a spectral constraint that targets k connected components. These methods achieve linear sample-size scaling for fixed anchor and optimization budgets. Their graph weights are learned, but their reconstruction is performed using supplied features and an anchor dictionary. Consequently, they do not jointly adapt the image representation to the self-expression objective. LoomSC preserves a compact representation while learning both the encoder and the latent factor. After initialization, the columns of the latent factor are updated from all samples rather than remaining fixed anchor representatives. Unlike the nonnegative anchor weights in LMVSC and SGL, its sample-factor coefficients can be signed. This permits linear combinations beyond convex anchor assignments. Selective sampling-based scalable sparse subspace clustering (S5C) [29] selects representatives using stochastic estimates of subgradient violations. It then solves restricted sparse regressions for all samples. Its spectral stage exploits a graph with only O(nm)O(nm) nonzero entries, where m is the representative budget. Thus, S5C reduces both regression and spectral costs without requiring a dense affinity. However, its representation remains tied to selected samples in the feature space. LoomSC replaces this exemplar dictionary with adaptive latent directions and jointly refines the features used to construct them. A-DSSC [22] addresses a different source of difficulty: the compatibility between self-expression coefficients and spectral normalization. It first estimates coefficients and then obtains a doubly stochastic affinity through quadratically regularized optimal transport. A support-restricted solver exploits sparsity to accelerate this second stage. This gives a principled alternative to heuristic coefficient normalization, but retains separate coefficient-estimation and transport problems. LoomSC instead derives the affinity and its normalization from the same learned projector. For its squared affinity, each degree equals the squared norm of the corresponding sample-factor row. Thus, degree computation requires only the thin factor and no additional transport optimization. Compact spectral computation also appears in scalable attributed-graph subspace clustering (SAGSC) [30]. SAGSC constructs a projector from graph-smoothed features and uses a quadratic kernel to obtain a nonnegative factorized affinity. This permits spectral clustering without explicitly assembling pairwise weights. LoomSC applies this spectral principle to jointly learned image representations rather than fixed graph-smoothed features. Its two-factor objective allows the projector to evolve with the encoder, and its analysis connects the learned factors to dominant latent directions and subspace structure. The distinction is therefore the integration of adaptive representation learning with a compact spectral construction, rather than factorization alone. I-B Deep subspace clustering: DSCNet [25] integrates a differentiable self-expression layer into a convolutional auto-encoder. This couples feature learning with the reconstruction of each latent sample from the others. S2ConvSCN [26] additionally uses spectral assignments to supervise feature learning and coefficient estimation. These formulations make the representation responsive to clustering, but retain a sample-to-sample self-expression matrix. Additional supervision does not remove its quadratic parameter storage. LoomSC retains joint feature and self-expression learning while replacing this layer with two thin factors. Its factor updates admit closed-form solutions, and their dimensions are controlled by the prescribed factor width. Neural coefficient prediction offers another route to scalability. SENet [23] parameterizes pairwise coefficients with a neural network that can generalize to unseen samples. Its objective learns a coefficient-prediction rule for the supplied representation rather than jointly learning a new feature space for self-expression. This separates the number of network parameters from the number of coefficient entries. However, a compact predictor does not itself provide a compact spectral representation of the resulting affinity. LoomSC instead learns a sample factor whose algebraic structure supports both self-expression and the subsequent spectral reduction. PRO-DSC [28] jointly learns latent representations and self-expression coefficients. A log-determinant regularizer controls representation collapse under the conditions established in its analysis. Its implementation predicts coefficients through neural mappings and Sinkhorn normalization, which supports mini-batch training. The distinction from LoomSC arises at full-set clustering. PRO-DSC constructs the coefficient matrix over the clustering set and applies spectral clustering to the resulting affinity. For n samples, this inference route retains quadratic matrix storage despite efficient training. LoomSC does not reconstruct a full coefficient matrix after learning. The final sample factor directly defines a nonnegative affinity whose normalized spectral problem is solved in a feature space determined by the factor width. Thus, the computational reduction extends through spectral assignment rather than ending with network training. Other deep methods reduce cost by changing how the final partition is obtained. The large-scale extension of pseudo-supervised deep subspace clustering [27] trains on a subset and assigns the remaining latent samples through nearest-neighbor classification. LoomSC instead retains an explicit model of sample relations while avoiding their dense storage. All samples participate in factor learning, and the final partition is obtained by spectral clustering of the learned affinity rather than by extending labels from a training subset. Its efficiency follows from representing that affinity compactly while preserving the self-expression formulation. I LoomSC The proposed method develops a compact self-representation by exploiting the spectral structure of least-squares regression (LSR). We first analyze the LSR solution and motivate a projector approximation. We then express this approximation through two thin factors that are learned jointly with a deep auto-encoder. Finally, we construct a nonnegative affinity from the learned factors and obtain the spectral embedding through a reduced eigenproblem. Both phases operate directly on the factors, which makes their computational cost linear. Let the matrix X=[x1,…,xn]∈ℝdx×nX=[x_1,…,x_n] ^d_x× n contain the input samples xi∈ℝdxx_i ^d_x, and let k denote the number of clusters. I-A A Projector Approximation Motivated by LSR LSR [9] models each sample through a linear combination of the data collection and penalizes the energy of the coefficients. For the input X, its diagonal-unconstrained formulation is minC∈ℝn×n12‖X−XC‖F2+λ‖C‖F2,λ>0. _C ^n× n\ 12\|X-XC\|_F^2+λ\|C\|_F^2, λ>0. (2) The reconstruction term encourages self-expression. The Frobenius penalty controls coefficient magnitudes while allowing correlated samples to contribute jointly. Differentiation with respect to C gives the closed-form solution of this problem CLSR⋆=(X⊤X+2λIn)−1X⊤X,C _LSR=(X X+2λ I_n)^-1X X, (3) where InI_n is the n×n× n identity matrix. A direct evaluation of Eq. (3) requires forming the Gram matrix X⊤X X, at a cost of O(dxn2)O(d_xn^2), followed by the solution of an n×n× n linear system, which requires O(n3)O(n^3) operations using a dense factorization. Moreover, the resulting coefficient matrix contains n2n^2 entries and therefore requires O(n2)O(n^2) memory. Although the inversion cost can be reduced in particular dimensional regimes, explicitly constructing the dense self-representation remains quadratic in the number of samples. Thereofre, we examine the spectral structure of CLSR⋆C _LSR to derive a compact representation. Let r=rank(X)r=rank(X). Since X⊤X X is symmetric positive semidefinite with rank r, its compact eigendecomposition is X⊤X=Vxdiag(β1,…,βr)Vx⊤,X X=V_xdiag( _1,…, _r)V_x , (4) where Vx∈ℝn×rV_x ^n× r has orthonormal columns given by the eigenvectors of X⊤X X associated with its positive eigenvalues, and β1≥⋯≥βr>0 _1≥·s≥ _r>0 are the corresponding eigenvalues. Proposition 1 (Spectral structure of LSR). For λ>0λ>0, the LSR solution is symmetric and positive semidefinite, with rank r. Its nonzero eigenvalues are τj=βjβj+2λ,j=1,…,r, _j= _j _j+2λ, j= \1,…,r \, (5) and CLSR⋆=Vxdiag(τ1,…,τr)Vx⊤C _LSR=V_xdiag( _1,…, _r)V_x . The eigenvectors identify the sample relationships represented by LSR. Eq. (5) also shows that directions with βj _j large relative to 2λ2λ have eigenvalues close to one. This observation motivates replacing their spectral weights by one and retaining the corresponding subspace. For the full row space, define Crow=VxVx⊤∈ℝn×nC_row=V_xV_x ^n× n. The approximation error is ‖CLSR⋆−Crow‖F2=∑j=1r(2λβj+2λ)2.\|C _LSR-C_row\|_F^2= _j=1^r ( 2λ _j+2λ )^2. (6) Thus, the projector approaches the LSR representation as the regularization becomes small relative to the positive data spectrum. This gives a direct spectral justification for an idempotent self-representation. To obtain a compact model, choose a factor width m with k≤m<nk≤ m<n. When m≤rm≤ r, let Vx,m∈ℝn×mV_x,m ^n× m contain the first m columns of VxV_x. The corresponding approximation satisfies ‖CLSR⋆−Vx,mVx,m⊤‖F2=∑j=1m(1−τj)2+∑j=m+1rτj2.\|C _LSR-V_x,mV_x,m \|_F^2= _j=1^m(1- _j)^2+ _j=m+1^r _j^2. (7) For the retained directions j≤mj≤ m, the approximation replaces the LSR eigenvalue τj _j by the projector eigenvalue 11. The first term in Eq. (7) measures the resulting error. The remaining directions j>mj>m are discarded, and the second term measures the corresponding loss of the LSR spectral components. Since τj _j is monotonically increasing with βj _j, the ordering of the τj _j’s is identical to that of the positive eigenvalues of X⊤X X. Hence, for a fixed width m, retaining the first m eigenvectors gives the closest projector among the rank-m projectors formed from the eigenvectors of X⊤X X. Setting m=rm=r recovers the full row-space projector Crow=VxVx⊤C_row=V_xV_x . However, constructing this projector requires determining r and computing all r nonzero right singular vectors of X. For a dense X∈ℝdx×nX ^d_x× n, computing the full singular subspace costs O(dxnmindx,n)O(d_xn \d_x,n\), while storing VxV_x requires O(nr)O(nr) memory. These costs become prohibitive when r grows with the data size. Thus, LoomSC does not attempt to recover the complete rank-r projector. Instead, m is treated as a prescribed model width, and the corresponding thin factors are learned directly in the subsequent formulation. Accordingly, we parameterize the self-representation as C=PP⊤,P∈ℝn×m,P⊤P=Im.C=P , P ^n× m, P P=I_m. (8) The factor P stores the relationships between all samples through nmnm coefficients. Its orthogonality makes C symmetric, positive semidefinite, and idempotent, with rank m. It also fixes the coefficient energy: ‖PP⊤‖F2=tr(P⊤P)=m.\|P \|_F^2=tr(P P)=m. (9) Substituting Eq. (8) into the LSR objective in Eq. (2) yields the following optimization over fixed-width projectors: minP⊤P=Im‖X−XPP⊤‖F2. _P P=I_m\ \|X-XPP \|_F^2. (10) Since PP⊤P is an orthogonal projector, Eq. (10) is equivalent to maximizing tr(P⊤X⊤XP)tr(P X XP). Hence, for m≤rm≤ r, an optimal solution is P⋆=Vx,mRP =V_x,mR, where R∈ℝm×mR ^m× m is any orthogonal matrix. All such solutions span the same dominant m-dimensional subspace and yield P⋆P⋆⊤=Vx,mVx,m⊤.P P =V_x,mV_x,m . (11) The minimum reconstruction error is ∑j=m+1rβj _j=m+1^r _j. Thus, the exact solution requires computing the dominant right singular subspace of the current representation. This is practical when the representation is fixed, but it becomes less attractive when the features are learned jointly, since the dominant subspace must be recomputed as the representation evolves. To address this problem, LoomSC replaces the repeated spectral solution with a two-factor formulation. The resulting updates involve only thin matrix products and an n×mn× m Procrustes problem, while the full coefficient matrix C remains implicit. I-B Learning a Latent Representation We apply the compact projector model to a learned latent space. This lets the self-expression structure guide feature extraction. We use an encoder ℰE and a decoder D with parameter collections W and W W, respectively. They define Z Z =ℰ(X,W)∈ℝdz×n, =E(X;W) ^d_z× n, (12) X X =(Z,W^)∈ℝdx×n, =D(Z; W) ^d_x× n, where dzd_z is the latent dimension and zi∈ℝdzz_i ^d_z is the representation of xix_i. The vector x^i∈ℝdx x_i ^d_x is the ithi^th column of the matrix X X. For image inputs, convolutional layers extract spatial features, which are flattened to form the columns of the latent representation matrix Z. The decoder reconstructs the input from these features. Its reconstruction objective is Jrec(W,W^)=1n‖X−X^‖F2.J_rec(W, W)= 1n\|X- X\|_F^2. (13) This objective trains the representation to retain information about the observed samples. The projector objective in Eq. (10), applied to Z, imposes a complementary subspace structure. I-C Self-Expression Through Two Thin Factors To optimize projector self-expression in the latent space, we introduce a latent factor L∈ℝdz×mL ^d_z× m. Its columns represent learned directions in the latent space. The following decomposition connects the resulting two-factor model to projector self-expression: ‖Z−LP⊤‖F2=‖Z−ZPP⊤‖F2+‖L−ZP‖F2.\|Z-LP \|_F^2=\|Z-ZPP \|_F^2+\|L-ZP\|_F^2. (14) The identity in Eq. (14) holds for every P satisfying P⊤P=ImP P=I_m. The proof is provided in Appendix B. Minimizing its left-hand side with respect to L gives L=ZPL=ZP. Hence, we obtain: minL,P⊤P=Im‖Z−LP⊤‖F2=minP⊤P=Im‖Z−ZPP⊤‖F2. _L,\,P P=I_m\|Z-LP \|_F^2= _P P=I_m\|Z-ZPP \|_F^2. (15) The two-factor formulation separates the latent directions represented by L from the sample coefficients encoded by P. This separation gives simple block updates and requires only products involving thin matrices. For fixed Z, let rz=rank(Z)r_z=rank(Z) and write Z=UzΣzVz⊤Z=U_z _zV_z . The matrices Uz∈ℝdz×rzU_z ^d_z× r_z and Vz∈ℝn×rzV_z ^n× r_z have orthonormal columns. The diagonal matrix Σz∈ℝrz×rz _z ^r_z× r_z contains the positive singular values σ1≥⋯≥σrz>0 _1≥·s≥ _r_z>0. Proposition 2 (Reconstruction through dominant directions). For a fixed latent matrix Z, the optimal factor residual is minL,P⊤P=Im‖Z−LP⊤‖F2=∑j=m+1rzσj2, _L,\,P P=I_m\|Z-LP \|_F^2= _j=m+1^r_z _j^2, (16) where an empty sum is zero. For m≤rzm≤ r_z, a minimizing P spans an m-dimensional dominant right singular subspace of Z. For m≥rzm≥ r_z, an optimal projector contains the row space of Z and satisfies ZPP⊤=ZZPP =Z. In particular, m=rzm=r_z gives C=VzVz⊤C=V_zV_z . Proposition 2 connects the learned factors to the dominant spectral directions of Z. Thus, the spectral approximation developed above applies directly to the learned representation. It also quantifies how the width controls the amount of latent information represented by the model. The proofs of (14) and Proposition 2 are given in Appendix B. We couple the factor objective with the auto-encoder to learn the latent features and their self-representation jointly: minW,W^,L,PJ(W,W^,L,P) _W, W,L,P J(W, W,L,P) =Jrec(W,W^)+‖Z−LP⊤‖F2, =J_rec(W, W)+\|Z-LP \|_F^2, (17) subject toP⊤P to P P =Im. =I_m. The reconstruction term retains information about the input samples. The factor term encourages the latent codes to admit a common compact self-representation. Joint optimization allows the encoder to adjust this representation as the factors are refined. For sample i, the row coefficients pi=Pi:⊤∈ℝmp_i=P_i: ^m reconstruct its latent vector as LpiLp_i. I-D Alternating Optimization We first pretrain the auto-encoder with (13). We then initialize the columns of L using m latent samples selected by k-means++ seeding [31]. This spreads the initial directions across the observed latent features. Appendix E specifies the sampling rule. The factor P is initialized by its update below. Each joint iteration updates the network parameters, refreshes Z, and then updates P and L in sequence. Network parameters With L and P fixed, an Adam step [32] updates the encoder and decoder: (W,W^)←Adamt((W,W^),∇W,W^J),(W, W) _t\! ((W, W), _W, WJ ), (18) where t is the joint iteration index. The factor term contributes the latent gradient 2(Z−LP⊤)2(Z-LP ), which encourages the encoder to produce representations consistent with the current factors. To control activation memory, the samples can be processed in chunks ℬb\B_b\ while keeping L and P fixed: J=∑b∑i∈ℬb(1n‖xi−x^i‖22+‖zi−Lpi‖22).J= _b _i _b ( 1n\|x_i- x_i\|_2^2+\|z_i-Lp_i\|_2^2 ). (19) The gradients are accumulated over all chunks before a single update of (W,W^)(W, W). Hence, chunking reduces memory usage without changing the full-data objective. After the update, the encoder recomputes Z before updating P and L. Sample factor For fixed Z and L, orthogonality gives ‖Z−LP⊤‖F2=‖Z‖F2+‖L‖F2−2tr(P⊤Z⊤L).\|Z-LP \|_F^2=\|Z\|_F^2+\|L\|_F^2-2tr(P Z L). (20) Thus, updating P amounts to maximizing the trace term under P⊤P=ImP P=I_m. This is an orthogonal Procrustes problem. To solve it, we form the thin product Z⊤L∈ℝn×mZ L ^n× m and compute Z⊤L=UPΣPVP⊤,P←UPVP⊤,Z L=U_P _PV_P , P← U_PV_P , (21) where UP∈ℝn×mU_P ^n× m and VP∈ℝm×mV_P ^m× m have orthonormal columns. The diagonal matrix ΣP∈ℝm×m _P ^m× m contains the singular values. This update attains the maximum trace and preserves the orthogonality constraint. Latent factor For fixed Z and P, differentiation gives ∇L‖Z−LP⊤‖F2=2(L−ZP) _L\|Z-LP \|_F^2=2(L-ZP). The minimizing update is L←ZP.L← ZP. (22) Each column of L is refined using the current latent samples and their coefficients. Together, the two updates adapt the representation to the full dataset using an n×mn× m decomposition and thin matrix products. Proposition 3 (Descent of the factor updates). For fixed Z, each update in (21) and (22) minimizes the factor objective with respect to its block. Repeating these updates yields a nonincreasing sequence of objective values that converges to a finite limit. Appendix B proves this result. Algorithm 1 combines the factor updates with representation learning for a prescribed number of joint iterations. Algorithm 1 LoomSC 1: X∈ℝdx×nX ^d_x× n, k≤m<nk≤ m<n, Tpre,T_pre,T. 2: Cluster labels γ1,…,γn _1,…, _n. 3: Pretrain (W,W^)(W, W) for TpreT_pre updates using (13). 4: Compute Z using (12). 5: Initialize L with m sample columns selected by applying k-means++ to Z. 6: Initialize P using (21). 7: for t=1,…,Tt=1,…,T do 8: Fix L and P and accumulate the gradient using (19). 9: Update (W,W^)(W, W) using (18). 10: Recompute Z using (12). 11: Update P using the thin-SVD step in (21). 12: Update L←ZPL← ZP using (22). 13: end for 14: Obtain the cluster labels from Algorithm 2 using the final P and Z. 15: return γ1,…,γn _1,…, _n. I-E Subspace Structure of the Learned Projector The projector interpretation also explains how self-expression can preserve subspace memberships. Consider latent samples spanning subspaces S1,…,Sk⊆ℝdzS_1,…,S_k ^d_z with dimensions r1,…,rkr_1,…,r_k. For independent subspaces, their sum is direct and dim(S1+⋯+Sk)=∑a=1kra=rz. (S_1+·s+S_k)= _a=1^kr_a=r_z. (23) Proposition 4 (Subspace preservation and connectivity). Suppose all latent samples are nonzero and span the direct sum in (23). If m=rz<nm=r_z<n and the fixed-Z factor objective attains its minimum, then C=PP⊤C=P is block diagonal after grouping samples by subspace. Within a subspace, its nonzero off-diagonal coefficients define a connected graph whenever that subspace’s samples admit no partition into two nonempty sets whose spans form a direct sum. The result identifies conditions under which the projector both separates the generating subspaces and connects the samples within each subspace. Its proof uses the block structure of the row-space projector and appears in Appendix C. The next phase converts these coefficients into nonnegative graph weights while preserving their nonzero connections. I-F Scalable Spectral Clustering from the Learned Factors The learned factor P provides the pairwise coefficients through C=PP⊤C=P . We now use this structure to obtain the final partition directly from compact matrices. The construction has two steps. First, we convert the coefficients into nonnegative affinities while preserving their support. We then express the normalized spectral problem through a quadratic feature map whose dimension depends only on m. A nonnegative affinity with an explicit feature map The magnitude of Cij=pi⊤pjC_ij=p_i p_j measures the strength of the relation between samples i and j. Squaring this coefficient gives the symmetric affinity A=C∘C∈ℝn×n,Aij=(pi⊤pj)2,A=C C ^n× n, A_ij=(p_i p_j)^2, (24) where ∘ denotes entrywise multiplication. This construction assigns nonnegative weights to both positive and negative coefficients. It also preserves their zero pattern. Thus, any separation between subspaces encoded by C is retained in A. Moreover, the affinity is invariant to an orthogonal change of basis in P. The quadratic form in (24) admits an explicit feature representation. Let q=m(m+1)/2q=m(m+1)/2. For a symmetric matrix S∈ℝm×mS ^m× m, svec(S)∈ℝqsvec(S) ^q contains its m diagonal entries followed by its upper-triangular entries. Each off-diagonal entry is multiplied by 2 2. This scaling preserves the Frobenius inner product between symmetric matrices. We define ϕ(p) φ(p) =svec(pp⊤)∈ℝq, =svec(p ) ^q, (25) Fi: F_i: =ϕ(pi)⊤,F∈ℝn×q. =φ(p_i) , F ^n× q. Consequently, ϕ(pi)⊤ϕ(pj)=(pi⊤pj)2φ(p_i) φ(p_j)=(p_i p_j)^2, which gives A=FF⊤.A=F . (26) Hence, the n×n× n affinity matrix is represented exactly through q features per sample and does not need to be formed explicitly. We retain the diagonal entries of A as self-loop weights. The projector identity C2=C^2=C then gives a particularly simple degree formula, δi=∑j=1nAij=(C2)ii=Cii=‖pi‖22,∑i=1nδi=m. _i= _j=1^nA_ij=(C^2)_i=C_i=\|p_i\|_2^2, _i=1^n _i=m. (27) Each degree therefore requires only a row norm of P. Computing all degrees takes O(nm)O(nm) time. Let ℐ+=i:δi>0I_+=\i: _i>0\ denote the indices with positive degrees and let n+=|ℐ+|n_+=|I_+|. Spectral normalization is performed on these samples. Zero-degree samples receive labels through their latent features after clustering. An exact reduced spectral problem Let F+∈ℝn+×qF_+ ^n_+× q contain the rows indexed by ℐ+I_+ in increasing order. Define Δ =diag(δi:i∈ℐ+)∈ℝn+×n+, =diag( _i:i _+) ^n_+× n_+, (28) B B =Δ−1/2F+∈ℝn+×q, = ^-1/2F_+ ^n_+× q, G G =B⊤B∈ℝq×q. =B B ^q× q. The normalized affinity is BB⊤B , and its normalized Laplacian is In+−BB⊤I_n_+-B . Thus, the smallest Laplacian eigenvalues correspond to the largest eigenvalues of BB⊤B . The following result transfers their computation to G. Proposition 5 (Exact spectral reduction). For P⊤P=ImP P=I_m, the normalized feature matrix satisfies m≤rank(B)≤minn+,q.m (B)≤ \n_+,q\. (29) The positive eigenvalues of G and BB⊤B coincide, including multiplicities. If Gvj=νjvjGv_j= _jv_j with νj>0 _j>0 and orthonormal vj∈ℝqv_j ^q, then Bvj/νjBv_j/ _j are orthonormal eigenvectors of BB⊤B with eigenvalues νj _j. The proof is provided in Appendix D. Since k≤mk≤ m, the reduced problem supplies at least k positive eigenvalues. The embedding dimension is k, whereas m controls the factor width and q determines the reduced spectral dimension. The degree identity also gives one dominant eigenpair analytically, v1=svec(Im)m,ν1=1,Bv1=(δi)i∈ℐ+m.v_1= svec(I_m) m, _1=1, Bv_1= ( _i)_i _+ m. (30) We retain this vector and compute the largest k−1k-1 eigenpairs on v1⟂v_1 . This gives k dominant modes while providing a strictly positive coordinate for every active sample. The sample embedding and its normalized rows are H H =B[v1,…,vk]diag(ν1,…,νk)−1/2∈ℝn+×k, =B[v_1,…,v_k]diag( _1,…, _k)^-1/2 ^n_+× k, (31) H^i: H_i: =Hi:/∥Hi:∥2,i=1,…,n+. =H_i:/\|H_i:\|_2, i=1,…,n_+. Equation (30) ensures that every row norm is positive. Retaining the eigenvalue-one modes also preserves the component indicators when the affinity graph has k connected components. In that case, row normalization maps each component to a distinct embedding point, as established in Appendix D. Following normalized spectral clustering [33], we apply k-means to the rows of H H. Let γi∈1,…,k _i∈\1,…,k\ denote the resulting label of active sample i. For a zero-degree sample, we use its distance to the latent centers of these clusters, a _a =i∈ℐ+:γi=a, =\i _+: _i=a\, (32) ζa _a =1|a|∑i∈azi∈ℝdz, = 1|J_a| _i _az_i ^d_z, γi _i ∈argmina∈1,…,k∥zi−ζa∥22,i∉ℐ+. ∈ _a∈\1,…,k\\|z_i- _a\|_2^2, i _+. The clustering conventions in Appendix E ensure nonempty active clusters and specify initialization and tie-breaking. Together, these steps assign every sample to a cluster. The corresponding membership matrix is Γ∈0,1n×k ∈\0,1\^n× k with Γia=γi=a _ia=1\ _i=a\. Direct and implicit implementations When q is modest, we accumulate G from normalized feature rows and solve its symmetric eigenproblem directly. For larger q, the same operator can be applied by streaming the rows of P. Let smat:ℝq→ℝm×msmat:R^q ^m× m be the inverse of svecsvec on symmetric matrices. For any v∈ℝqv ^q, Gv=svec(∑i∈ℐ+pi⊤smat(v)piδipipi⊤).Gv=svec\! ( _i _+ p_i smat(v)p_i _ip_ip_i ). (33) Each operator evaluation takes O(nm2)O(nm^2) time. A symmetric partial eigensolver uses this action on v1⟂v_1 to recover the remaining dominant modes. Both implementations solve the spectral problem defined by (24). They use compact factors throughout, with C, A, and the sample Laplacian represented implicitly. Algorithm 2 Spectral clustering through the quadratic feature map 1: P∈ℝn×mP ^n× m with P⊤P=ImP P=I_m, Z∈ℝdz×nZ ^d_z× n, k≤mk≤ m, Rkm,TkmR_km,T_km. 2: Cluster labels γ1,…,γn _1,…, _n. 3: Compute δi _i by (27) and retain the indices ℐ+I_+ with δi>0 _i>0. 4: Set q=m(m+1)/2q=m(m+1)/2 and the feature ordering in (25). 5: Set v1,ν1v_1, _1 using (30). 6: if k>1k>1 then 7: if direct eigensolver then 8: Accumulate G=∑i∈ℐ+ϕ(pi)ϕ(pi)⊤/δiG= _i _+φ(p_i)φ(p_i) / _i using (28). 9: Compute its largest k−1k-1 eigenpairs on v1⟂v_1 . 10: else 11: Compute the largest k−1k-1 eigenpairs on v1⟂v_1 with a symmetric partial eigensolver using (33). 12: end if 13: end if 14: Recover H and normalize its rows by (31). 15: Run RkmR_km seeded k-means restarts using (54), with at most TkmT_km iterations each. 16: Retain the restart with minimum distortion and restore the original sample indices. 17: Assign zero-degree samples by (32). 18: return γ1,…,γn _1,…, _n. I-G Complete Algorithm and Computational Complexity Algorithm 1 combines representation learning with spectral clustering. Let TpreT_pre and T be the numbers of pretraining updates and joint outer iterations, respectively. The spectral routine is given in Algorithm 2. It uses RkmR_km clustering restarts with at most TkmT_km Lloyd iterations per restart. Representation learning Let cnetc_net bound the per-sample cost of network evaluation and backpropagation in one update, including the refreshed encoder evaluation. Pretraining takes O(Tprencnet)O(T_prenc_net) time. Seeding requires O(ndzm)O(nd_zm) time, and the initial Procrustes update requires O(ndzm+nm2+m3)O(nd_zm+nm^2+m^3). Each joint iteration takes O(ncnet+ndzm+nm2+m3).O\! (nc_net+nd_zm+nm^2+m^3 ). (34) These terms account for network training, the products Z⊤LZ L and ZPZP, and the thin SVD. The representation cost is the sum of pretraining and initialization costs plus T times (34). Gradient accumulation allows the activation workspace to be controlled through the chunk size. Spectral computation Computing the degrees requires O(nm)O(nm) time. In the direct implementation, constructing the reduced matrix G∈ℝq×qG ^q× q costs O(nq2)O(nq^2), and its eigendecomposition costs O(q3)O(q^3). Recovering the k-dimensional sample embedding requires O(nqk)O(nqk) operations. The total complexity is therefore O(nm+nq2+q3+nqk),q=m(m+1)/2.O(nm+nq^2+q^3+nqk), q=m(m+1)/2. (35) The dependence on n appears only in the construction of the reduced matrix and the recovery of the sample embedding. The eigendecomposition is performed in the q-dimensional space. For the implicit implementation, the matrix G is not formed explicitly. Let TeigT_eig denote the total number of matrix–vector products required by the iterative eigensolver, seigs_eig the maximum Krylov basis size, and ReigR_eig the number of restarts. Including operator evaluations, orthogonalization, and the projected eigenproblems, the total complexity is O(nm+Teig(nq+qseig)+Reig(qseig2+seig3)+nqk).O\! (nm+T_eig(nq+qs_eig)+R_eig(qs_eig^2+s_eig^3)+nqk ). (36) Each vector processed by a block eigensolver contributes one matrix–vector product to TeigT_eig. The quantities TeigT_eig, seigs_eig, and ReigR_eig depend on the convergence of the eigensolver and the required numerical accuracy. Since the implicit implementation evaluates products with G only when needed, it avoids storing the dense q×q× q matrix and is preferable when q is large. Clustering and memory The k-dimensional embedding gives a clustering cost of O(Rkm(Tkm+1)nk2)O(R_km(T_km+1)nk^2). Computing latent centers and assigning zero-degree samples adds at most O(ndzk)O(nd_zk). Let MnetM_net include network parameters, optimizer state and the activation workspace. A resident-data memory bound is O(n(dx+dz+m+k)+dzm+Mnet),O\! (n(d_x+d_z+m+k)+d_zm+M_net ), (37) plus O(q2+qk)O(q^2+qk) for direct reduction or O(qseig+qk)O(qs_eig+qk) for the implicit eigensolver. Quadratic feature rows can be evaluated as needed, so storing F is optional. The output labels require O(n)O(n) memory. Together, these bounds establish linear time and memory in n for fixed dimensions, architecture, factor width, cluster count and iteration budgets. The factor width m is a chosen model parameter, so both phases operate without estimating the data rank. Their shared use of compact factors makes scalability extend from representation learning to the final partition. IV Results and Discussion Protocol We evaluate LoomSC on five image-clustering benchmarks. The face datasets are Extended YaleB (EYaleB), ORL, and UMIST, with (n,k)=(2,432,38)(n,k)=(2,432,38), (400,40)(400,40), and (575,20)(575,20), respectively. We also use COIL-100 with 7,2007,200 images from 100100 object classes and a balanced subset of Fashion-MNIST with 5,0005,000 images from 1010 classes. Each sample is normalized to unit ℓ2 _2 norm. Tables I and IV specify the sample counts used in the comparisons. We compare with nine baselines that distinguish deep representation learning from computational scalability. The scalable non-deep methods comprise two groups. LMVSC [15], SGL [16], and S5C [29] use sampling or anchor-based reduction. EnSC [19], SSC-OMP [18], A-DSSC [22], and SSSC [20] accelerate self-expression estimation without learning deep representations. These methods target scalability, but do not all achieve linear complexity for the complete clustering pipeline. The non-linear-time full-affinity references are LSR [9] and PRO-DSC [28]. LSR is non-deep and computes a dense self-expression matrix. PRO-DSC is deep and supports mini-batch training. Its published spectral inference nevertheless constructs pairwise coefficients over the clustering set and processes the resulting affinity. Thus, scalable training does not eliminate its full-affinity spectral bottleneck. LoomSC combines deep representation learning with factorized self-expression and spectral assignment. Its complete pipeline has linear time and memory complexity in n for fixed architecture, dimensions, factor width, and iteration budgets. A detailed comparison of computational and matrix-storage complexities is provided in Appendix F. For LoomSC, we use the same architecture and optimization hyperparameters across all datasets. The encoder contains L=3L=3 convolutional layers with 1010, 2020, and 3030 output channels and kernel sizes 5×55× 5, 3×33× 3, and 3×33× 3, respectively. Each layer uses stride 22 and a ReLU activation, with no batch normalization. The decoder mirrors the encoder using transposed convolutions. The resulting convolutional representation has size 30×6×630× 6× 6, which gives a flattened latent dimension of dz=1080d_z=1080. We set the factor width to m=380m=380 for all datasets. The auto-encoder is first pretrained for 5,0005,000 epochs using full-batch Adam with a learning rate of 10−310^-3. No data augmentation, weight decay, or learning-rate scheduling is used. After pretraining, k-means with k-means++ initialization is applied to the latent representations, and its m centroids initialize the columns of L. The factor P is then initialized by the orthogonal Procrustes solution. Joint optimization is performed for another 5,0005,000 epochs using full-batch Adam with a learning rate of 10−310^-3. The reconstruction and self-expression terms have equal weights. At each joint-training epoch, the network parameters are first updated by Adam. The factor P is then updated by the exact Procrustes solution, followed by the closed-form update L=ZPL=ZP, as described in Algorithm 1. After joint optimization, spectral assignment is performed following Algorithm 2. The spectral embedding dimension is set to the target number of clusters k, after which k-means produces the final partition. We report clustering accuracy (ACC), normalized mutual information (NMI), and adjusted Rand index (ARI) as percentages. Higher values indicate better clustering. Repeated-run results are summarized by their mean and standard deviation over three seeds. A scalar entry indicates that no standard deviation is reported. Differences between percentage scores are expressed in percentage points. The synthetic experiments additionally report the connectivity diagnostic (CONN) [34], for which higher values are preferred. IV-A Clustering quality Tables I–I show that LoomSC ranks first or second in all 1515 dataset–metric comparisons, including a shared second place for COIL-100 NMI. No other evaluated method maintains this ranking across all five datasets and three metrics. Its unweighted means are 80.17%80.17\% ACC, 85.98%85.98\% NMI, and 72.77%72.77\% ARI. These exceed the highest corresponding baseline means by 6.666.66, 4.094.09, and 10.0210.02 points. Thus, the compact formulation provides consistently strong clustering quality while retaining linear-in-n computation under the conditions stated above. TABLE I: Clustering quality — ACC (%) on the five benchmarks. Higher is better. Mean ± standard deviation over three seeds where available. Bold indicates the largest mean in each row. Non-deep methods Deep methods Dataset n k EnSC SSC-OMP LMVSC SGL S5C LSR A-DSSC SSSC PRO-DSC LoomSC Fashion-MNIST 5,000 10 61.35 ± 0.00 56.19 ± 1.36 59.18 54.47 ± 0.36 54.62 ± 1.99 59.52 ± 0.00 39.01 ± 2.20 42.01 ± 1.70 52.12 ± 2.92 65.33 ± 1.79 ORL 400 40 78.00 ± 1.59 73.75 ± 0.82 55.25 ± 0.00 62.25 ± 0.00 42.50 ± 0.00 82.50 ± 0.00 50.75 ± 0.00 46.50 ± 1.54 52.00 ± 7.91 86.17 ± 2.01 UMIST 575 20 72.52 ± 0.00 67.70 ± 0.00 60.62 ± 0.00 65.00 ± 0.00 40.87 ± 0.00 58.26 ± 0.00 52.52 ± 0.00 21.16 ± 1.31 81.39 ± 1.58 80.00 ± 2.19 EYaleB 2,432 38 88.19 ± 0.00 80.84 ± 1.34 21.50 ± 0.00 21.50 ± 0.00 22.41 ± 0.00 96.26 ± 0.00 61.88 ± 0.00 25.11 ± 2.27 77.63 ± 0.73 94.73 ± 1.26 COIL-100 7,200 100 67.50 ± 0.69 59.46 ± 0.00 49.50 ± 0.00 49.26 ± 0.00 42.04 ± 0.00 61.10 ± 0.00 48.57 ± 0.00 13.81 ± 1.31 82.31 ± 0.30 74.64 ± 1.34 TABLE I: Clustering quality — NMI (%) on the five benchmarks. Higher is better. Mean ± standard deviation over three seeds. Bold indicates the largest mean in each row. Non-deep methods Deep methods Dataset EnSC SSC-OMP LMVSC SGL S5C LSR A-DSSC SSSC PRO-DSC LoomSC Fashion-MNIST 60.29 ± 0.88 57.93 ± 0.41 64.37 ± 0.63 62.51 ± 0.42 54.82 ± 2.33 63.33 ± 0.00 53.83 ± 1.00 46.88 ± 0.45 56.45 ± 3.08 63.91 ± 0.82 ORL 89.25 ± 0.07 84.54 ± 0.37 74.53 ± 0.00 78.24 ± 0.00 64.76 ± 0.00 90.07 ± 0.00 77.93 ± 0.00 66.84 ± 1.34 71.71 ± 5.39 91.76 ± 0.88 UMIST 87.06 ± 0.00 63.61 ± 0.53 62.49 ± 0.00 56.84 ± 0.00 56.72 ± 0.00 75.40 ± 0.00 66.01 ± 0.00 32.89 ± 1.40 90.02 ± 0.20 87.83 ± 1.44 EYaleB 74.31 ± 0.03 83.17 ± 0.55 13.37 ± 0.00 6.44 ± 0.00 37.92 ± 0.00 95.44 ± 0.00 75.61 ± 0.00 34.24 ± 2.42 82.38 ± 0.91 95.47 ± 0.43 COIL-100 90.95 ± 0.02 57.11 ± 0.16 74.69 ± 0.00 74.54 ± 0.00 68.03 ± 0.00 85.21 ± 0.00 74.31 ± 0.00 37.80 ± 1.08 95.75 ± 0.05 90.95 ± 0.41 TABLE I: Clustering quality — ARI (%) on the five benchmarks. Higher is better. Mean ± standard deviation over three seeds. Bold indicates the largest mean in each row. Non-deep methods Deep methods Dataset EnSC SSC-OMP LMVSC SGL S5C LSR A-DSSC SSSC PRO-DSC LoomSC Fashion-MNIST 43.21 ± 1.04 42.00 ± 0.16 45.35 ± 0.87 42.75 ± 0.22 40.27 ± 1.50 46.53 ± 0.00 22.24 ± 1.29 19.41 ± 1.49 37.76 ± 4.07 51.89 ± 1.43 ORL 69.50 ± 0.55 60.22 ± 0.74 39.22 ± 0.00 47.60 ± 0.00 22.66 ± 0.00 74.10 ± 0.00 21.49 ± 0.00 15.32 ± 4.98 35.55 ± 8.75 78.05 ± 2.67 UMIST 66.91 ± 0.00 34.81 ± 0.47 32.12 ± 0.00 25.19 ± 0.00 26.63 ± 0.00 46.89 ± 0.00 22.99 ± 0.00 1.07 ± 0.31 79.36 ± 0.65 75.55 ± 2.69 EYaleB 19.67 ± 0.06 53.91 ± 5.11 1.87 ± 0.00 -0.40 ± 0.00 13.65 ± 0.00 92.19 ± 0.00 21.67 ± 0.00 2.16 ± 0.31 64.15 ± 1.28 91.76 ± 0.90 COIL-100 61.38 ± 1.84 21.98 ± 0.13 43.00 ± 0.00 43.93 ± 0.00 35.05 ± 0.00 54.03 ± 0.00 7.84 ± 0.00 1.59 ± 0.04 80.63 ± 1.28 66.60 ± 0.66 LoomSC achieves higher ACC than every scalable non-deep baseline on all five datasets. The improvements over the strongest method in this group range from 3.983.98 to 8.178.17 points. The gains over sampling-based methods are larger. Relative to the best ACC among LMVSC, SGL, and S5C, LoomSC improves by 23.9223.92 points on ORL and 25.1425.14 points on COIL-100. The corresponding gains are 15.0015.00 points on UMIST and 72.3272.32 points on EYaleB. These comparisons show that a compact model need not be restricted to the clustering quality of fixed-representation sampling methods. LoomSC retains compactness while jointly refining the latent features and their self-expression factors. On Fashion-MNIST, LoomSC obtains the highest ACC of 65.33±1.79%65.33± 1.79\% and ARI of 51.89±1.43%51.89± 1.43\%. These exceed the strongest competing results by 3.983.98 and 5.365.36 points, respectively. Its NMI of 63.91±0.82%63.91± 0.82\% is within 0.460.46 points of the highest value. The agreement between ACC and ARI supports the quality of the partition beyond a single evaluation criterion. On ORL, LoomSC leads all three metrics with 86.17±2.01%86.17± 2.01\% ACC, 91.76±0.88%91.76± 0.88\% NMI, and 78.05±2.67%78.05± 2.67\% ARI. It improves upon LSR, the strongest alternative for each metric, by 3.673.67, 1.691.69, and 3.953.95 points. These results are obtained with only ten images per subject. On EYaleB, LoomSC reaches 94.73%94.73\% ACC, 95.47%95.47\% NMI, and 91.76%91.76\% ARI. Its NMI closely matches LSR at 95.44%95.44\%, and their ARI scores differ by only 0.430.43 points. Thus, the thin-factor representation can match or improve upon dense least-squares self-expression on several clustering measures without storing all pairwise coefficients. The comparison with PRO-DSC assesses whether this efficiency is compatible with deep clustering quality. LoomSC exceeds PRO-DSC on all three metrics for Fashion-MNIST, ORL, and EYaleB. The ACC gains are 13.2113.21, 34.1734.17, and 17.1017.10 points, respectively. Across the five datasets, its average improvements over PRO-DSC are 11.0811.08 ACC points, 6.726.72 NMI points, and 13.2813.28 ARI points. Therefore, the advantage over the non-deep baselines is accompanied by a strong comparison with a deep full-affinity method. PRO-DSC attains the highest scores on UMIST and COIL-100, where LoomSC ranks second on all three metrics. On UMIST, LoomSC reaches 80.00%80.00\% ACC compared with 81.39%81.39\% for PRO-DSC. On COIL-100, it obtains 74.64%74.64\% ACC, 90.95%90.95\% NMI, and 66.60%66.60\% ARI. Its COIL-100 ACC exceeds every non-deep baseline, with gains of 7.147.14 points over EnSC and 13.5413.54 points over LSR. These results establish competitive deep clustering performance on both datasets. As quantified next, LoomSC attains this performance at a small fraction of PRO-DSC’s measured runtime and without its non-linear-time full-affinity spectral stage. IV-B Computational efficiency Table IV reports the measured wall-clock times. LoomSC requires 0.2390.239 s on ORL and 0.2030.203 s on UMIST. Its mean time is 2.2402.240 s on EYaleB, 9.0089.008 s on Fashion-MNIST, and 25.31125.311 s on COIL-100. No competing method combines higher mean ACC with lower mean runtime on any of these datasets. Thus, LoomSC lies on the empirical ACC–runtime Pareto frontier for every benchmark among the evaluated implementations. TABLE IV: Measured wall-clock time (s). Lower is better. Bold indicates the smallest mean in each row. Non-deep methods Deep methods Dataset n EnSC SSC-OMP LMVSC SGL S5C LSR A-DSSC SSSC PRO-DSC LoomSC Fashion-MNIST 5,000 3219.752 ± 46.948 5.447 ± 0.681 16.965 ± 0.179 17.203 ± 0.219 8.280 ± 1.000 6.260 ± 0.157 35.612 ± 0.584 2231.314 ± 69.877 563.338 ± 1.392 9.008 ± 0.987 ORL 400 7.953 ± 0.514 1.856 ± 0.058 0.929 ± 0.031 0.958 ± 0.003 7.440 ± 0.419 0.049 ± 0.001 0.699 ± 0.005 7.084 ± 0.183 102.797 ± 2.671 0.239 ± 0.089 UMIST 575 24.366 ± 26.899 0.498 ± 0.012 1.086 ± 0.033 1.263 ± 0.043 10.249 ± 4.946 0.077 ± 0.002 1.026 ± 0.027 9.774 ± 1.423 78.689 ± 0.172 0.203 ± 0.003 EYaleB 2,432 43.843 ± 0.398 12.666 ± 0.975 9.124 ± 0.063 9.221 ± 0.083 471.768 ± 24.647 2.389 ± 0.542 13.918 ± 0.132 2286.421 ± 94.488 313.648 ± 2.670 2.240 ± 0.581 COIL-100 7,200 50.858 ± 1.523 1319.586 ± 164.863 21.421 ± 0.143 22.498 ± 0.390 4709.568 ± 169.023 49.888 ± 5.451 175.328 ± 1.551 2409.599 ± 245.492 7443.627 ± 4.992 25.311 ± 9.976 The distinction from PRO-DSC concerns the complete clustering pipeline rather than network training alone. PRO-DSC reduces training costs through mini-batches, but its published inference constructs a full pairwise affinity before spectral clustering [28]. Explicitly storing this affinity requires O(n2)O(n^2) memory for a clustering set of size n. LoomSC instead retains P∈ℝn×mP ^n× m and L∈ℝdz×mL ^d_z× m during representation learning. Its quadratic feature map transfers spectral computation to dimension q=m(m+1)/2q=m(m+1)/2, which is independent of n for fixed m. Consequently, it removes the full-affinity bottleneck rather than only accelerating its construction. Section I-G establishes linear time and memory in n under fixed dimensions and iteration budgets. The reported runtimes show substantial savings relative to the deep baseline. LoomSC is 62.5×62.5× faster than PRO-DSC on Fashion-MNIST while improving ACC by 13.2113.21 points. On ORL, it is 430.1×430.1× faster and improves ACC by 34.1734.17 points. On EYaleB, the corresponding speedup is 140.0×140.0×, together with a 17.1017.10-point ACC gain. Hence, on three datasets, the computational reduction accompanies higher clustering quality rather than a loss of accuracy. The datasets where PRO-DSC achieves its strongest results further illustrate the scalability advantage. On UMIST, LoomSC attains 80.00%80.00\% ACC in 0.2030.203 s, compared with 81.39%81.39\% in 78.68978.689 s for PRO-DSC. This is a 387.6×387.6× speedup with an ACC difference of 1.391.39 points. On COIL-100, LoomSC obtains the second-highest ACC while reducing the measured time from 7,443.6277,443.627 to 25.31125.311 s, a 294.1×294.1× speedup. Across all five datasets, the speedups over PRO-DSC range from 62.5×62.5× to 430.1×430.1×. Therefore, competitive deep clustering performance is retained without the runtime of the full-affinity reference. LoomSC also improves the quality–cost balance relative to scalable non-deep methods. On Fashion-MNIST, it exceeds EnSC by 3.983.98 ACC points with a 357.4×357.4× speedup. On COIL-100, it improves ACC over EnSC and LSR by 7.147.14 and 13.5413.54 points while approximately halving their runtimes. Compared with LMVSC, the fastest method on COIL-100, LoomSC requires only 3.893.89 additional seconds and improves ACC by 25.1425.14 points. These results show that LoomSC combines runtimes close to or below those of sampling methods with substantially stronger clustering. Its compact pipeline also remains far faster than the deep full-affinity reference. Fig. 1: Synthetic-data comparison with ten subspaces and nin_i samples per subspace. Panels show (a) ACC, (b) connectivity, (c) CPU runtime, and (d) GPU runtime. “Ours” denotes LoomSC. IV-C Scalability Analysis We first vary the sampling density in a controlled union-of-subspaces setting. Following the synthetic-data configuration in [20], we generate k=10k=10 subspaces with random dimensions between 66 and 1212 in an ambient space of dimension dx=784d_x=784. Each subspace contains ni∈30,50,100,200,500,1000n_i∈\30,50,100,200,500,1000\ samples. The total sample count is therefore n=10nin=10n_i, ranging from 300300 to 10,00010,000. Figure 1 compares ACC, connectivity, and CPU and GPU runtimes. LoomSC maintains 99%99\%–100%100\% ACC across the entire density range. Its connectivity remains high, at 0.680.68–0.750.75. Several competing methods attain high ACC at sufficiently large sample counts but have much smaller connectivity scores. Thus, the two diagnostics provide complementary evidence: LoomSC produces accurate partitions while retaining strong graph connectivity under the reported measure. This behavior is consistent with the goal of preserving informative sample relations through the factorized representation. The comparison with the SDSNet [35] reference is particularly informative at low sampling densities. With 3030 samples per subspace, LoomSC reaches 99%99\% ACC, compared with about 85%85\% for SDSNet. With 5050 samples per subspace, its ACC is 100%100\%, compared with about 87%87\%. At the largest size, the two approaches both retain near-perfect accuracy and high connectivity. LoomSC nevertheless reduces the measured runtime from 143143 to 3939 s in both the CPU and GPU panels. The improvement is visible within each hardware setting rather than arising only from a CPU-to-GPU comparison. The same experiment contrasts LoomSC with the full-affinity deep baseline. PRO-DSC also reaches near-perfect ACC at higher sampling densities. However, its CPU runtime is substantially larger at the highest density. LoomSC combines high clustering accuracy across the full density range with lower computational cost at the largest size in this experiment. Figure 2 extends the sample-size experiment from 1,0001,000 to 500,000500,000 points. LoomSC maintains 99.8%99.8\% ACC or higher throughout this 500500-fold increase. Its runtime grows from about 5757 s at 50,00050,000 samples to 362362 s at 500,000500,000 samples. At the largest size, it retains near-perfect accuracy in approximately six minutes. S5C reaches comparable accuracy but requires about 1,0541,054 s. LMVSC has a lower runtime of about 125125 s, with ACC near 94%94\%. SGL also reaches near-perfect ACC at large sizes, whereas LoomSC maintains this accuracy regime across the full tested range, including the smallest datasets. Fig. 2: Synthetic-data scalability from 1,0001,000 to 500,000500,000 samples: (a) runtime and (b) clustering accuracy. The horizontal-axis symbol N denotes the total sample count n. The reference curves indicate linear and quadratic growth. “Ours” denotes LoomSC. The large-sample runtime trend is consistent with the linear-in-n bounds in Section I-G. For fixed factor width, the quadratic feature dimension q=m(m+1)/2q=m(m+1)/2 does not grow with the number of samples. Neither factor learning nor the reduced spectral stage requires storage of all pairwise affinities. With fixed dimensions and iteration budgets, the sample-dependent work remains linear. The experiment shows that this design supports increasing the dataset size without a corresponding loss of clustering accuracy. IV-D Convergence of the Affinity Matrix Figure 3 examines the relative Frobenius change between successive affinity iterates. The diagnostic divides the Frobenius norm of the difference by the norm of the previous iterate. This difference measures the stability of the sample relationships used for clustering, rather than the scale of the matrix alone. Fig. 3: Relative Frobenius change ‖At−At−1‖F/‖At−1‖F\|A_t-A_t-1\|_F/\|A_t-1\|_F on (a) synthetic datasets with ten subspaces and (b) real datasets. Here, AtA_t denotes the affinity matrix in (24) at outer iteration t, and nin_i is the number of samples in synthetic subspace i. On synthetic data, all six sampling densities exhibit a short initial transient followed by gradual stabilization. The relative change falls to a few 10−310^-3 within roughly the first 100100 outer iterations and approaches 10−310^-3 by the end of the run. The curves reach similar terminal scales despite the increase from 3030 to 1,0001,000 samples per subspace. Thus, the larger synthetic datasets retain the same qualitative stabilization pattern. The real-data curves show similarly small late-stage changes. UMIST falls below 10−310^-3 after about 100100 iterations. Fashion-MNIST reaches a low baseline with temporary excursions, while EYaleB settles around a few 10−310^-3. The dominant pattern is a rapid reduction in the update magnitude followed by sustained small changes. This supports the practical stability of the learned affinities during joint training. Proposition 3 provides a complementary descent result for the exact factor updates with fixed latent features. The observed curves show that small affinity changes are also attained when these updates are coupled with representation learning. (a) Epoch 10 (b) Epoch 100 (c) Epoch 500 Fig. 4: Learned sample relations on a ten-class EYaleB subset at epochs 1010, 100100, and 500500. Class boundaries are overlaid. (a) Epoch 10 (b) Epoch 100 (c) Epoch 500 Fig. 5: Learned sample relations on a ten-class ORL subset at epochs 1010, 100100, and 500500. Class boundaries are overlaid. (a) Epoch 10 (b) Epoch 100 (c) Epoch 500 Fig. 6: Learned sample relations on synthetic data with ten subspaces at epochs 1010, 100100, and 500500. Subspace boundaries are overlaid. (a) n=200n=200 (b) n=300n=300 (c) n=500n=500 Fig. 7: Learned sample relations on synthetic data with ten subspaces and increasing total sample count at epoch 500500. Subspace boundaries are overlaid. IV-E Affinity-matrix structure Figures 4–6 visualize the learned sample relations at epochs 1010, 100100, and 500500 after pretraining. The face-data panels use ten-class subsets, and the synthetic panels contain ten subspaces. Samples are grouped by their classes, with class boundaries overlaid. These plots examine whether the representation retains class-aligned relations throughout optimization. On EYaleB, the diagonal blocks are already visible after 10 epochs and remain clearly identifiable at later epochs. The within-class regions contain more pronounced structure than most off-block regions. ORL also retains recognizable diagonal blocks despite having only ten images per subject. This visual organization is consistent with the strong face-clustering scores in Tables I–I. In particular, the displayed relations extend beyond individual self-connections and include associations between distinct samples from the same class. The synthetic panels exhibit an even clearer separation between dense diagonal blocks and the surrounding regions. This organization persists across the training stages (i.e., 1010, 100100, and 500500 epochs). Figure 7 compares n=200n=200, 300300, and 500500 at epoch 500500. The 10 diagonal blocks remain identifiable as more samples are added, and the larger panels show a clearer visual contrast between within-subspace structure and off-block regions. Together with the quantitative connectivity results, these visualizations support the ability of the affinity matrix to retain structured sample relations as the dataset grows. V Conclusion This work presents LoomSC, a scalable subspace clustering framework that combines projector-based self-expression with exact spectral reduction. Jointly learning latent representations and two thin factors avoids explicit sample-by-sample coefficient matrices. A nonnegative quadratic affinity preserves the projector’s support and enables normalized spectral clustering without forming the full affinity or sample Laplacian. For fixed dimensions and iteration budgets, the complete pipeline has linear time and memory complexity in the number of samples. Our analysis quantifies the projector approximation to regularized least squares. It also identifies conditions under which the full row-space projector preserves subspace memberships and ensures within-subspace connectivity. Across five image-clustering benchmarks, LoomSC achieves the highest average ACC, NMI, and ARI among the evaluated methods. Synthetic experiments extend to 500,000500,000 samples while maintaining at least 99.8%99.8\% accuracy throughout the tested range. These findings demonstrate that compact self-expression can combine strong clustering quality with scalability across both representation learning and spectral assignment. References [1] Y. Ma, H. Derksen, W. Hong, and J. Wright (2007) Segmentation of multivariate mixed data via lossy data coding and compression. IEEE Trans. Pattern Anal. Mach. Intell. 29 (9), p. 1546–1562. External Links: Document Cited by: §I. [2] S. Rao, R. Tron, R. Vidal, and Y. Ma (2010) Motion segmentation in the presence of outlying, incomplete, or corrupted trajectories. IEEE Trans. Pattern Anal. Mach. Intell. 32 (10), p. 1832–1845. External Links: Document Cited by: §I. [3] R. Vidal (2011) Subspace clustering. IEEE Signal Process. Mag. 28 (2), p. 52–68. External Links: Document Cited by: §I. [4] T. Ding, D. Lim, R. Vidal, and B. D. Haeffele (2022) Understanding doubly stochastic clustering. In Proc. Int. Conf. Mach. Learn. (ICML), Proceedings of Machine Learning Research, Vol. 162, p. 5153–5165. Cited by: §I. [5] R. Basri and D. W. Jacobs (2003) Lambertian reflectance and linear subspaces. IEEE Trans. Pattern Anal. Mach. Intell. 25 (2), p. 218–233. External Links: Document Cited by: §I. [6] C. Lu, J. Feng, Z. Lin, T. Mei, and S. Yan (2019) Subspace clustering by block diagonal representation. IEEE Trans. Pattern Anal. Mach. Intell. 41 (2), p. 487–501. External Links: Document Cited by: §I, §I. [7] E. Elhamifar and R. Vidal (2009) Sparse subspace clustering. In Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 2790–2797. External Links: Document Cited by: §I. [8] E. Elhamifar and R. Vidal (2013) Sparse subspace clustering: algorithm, theory, and applications. IEEE Trans. Pattern Anal. Mach. Intell. 35 (11), p. 2765–2781. External Links: Document Cited by: §I, §I. [9] C. Lu, H. Min, Z. Zhao, L. Zhu, D. Huang, and S. Yan (2012) Robust and efficient subspace segmentation via least squares regression. In Eur. Conf. Comput. Vis. (ECCV), Lecture Notes in Computer Science, Vol. 7578, p. 347–360. External Links: Document Cited by: TABLE V, §I, §I, §I-A, §I-A, §IV. [10] G. Liu, Z. Lin, and Y. Yu (2010) Robust subspace segmentation by low-rank representation. In Proc. Int. Conf. Mach. Learn. (ICML), p. 663–670. Cited by: §I. [11] D. Luo, F. Nie, C. Ding, and H. Huang (2011) Multi-subspace representation and discovery. In Machine Learning and Knowledge Discovery in Databases, Lecture Notes in Computer Science, Vol. 6912, p. 405–420. External Links: Document Cited by: §I. [12] B. Nasihatkon and R. Hartley (2011) Graph connectivity in sparse subspace clustering. In Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 2137–2144. External Links: Document Cited by: §I. [13] X. Peng, L. Zhang, and Z. Yi (2013) Scalable sparse subspace clustering. In Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 430–437. Cited by: §I. [14] X. Peng, H. Tang, L. Zhang, Z. Yi, and S. Xiao (2016) A unified framework for representation-based subspace clustering of out-of-sample and large-scale data. IEEE Trans. Neural Netw. Learn. Syst. 27 (12), p. 2499–2512. External Links: Document Cited by: §I. [15] Z. Kang, W. Zhou, Z. Zhao, J. Shao, M. Han, and Z. Xu (2020) Large-scale multi-view subspace clustering in linear time. In Proc. AAAI Conf. Artif. Intell. (AAAI), Vol. 34, p. 4412–4419. External Links: Document Cited by: TABLE V, Appendix F, §I, §I, §I-A, §IV. [16] Z. Kang, Z. Lin, X. Zhu, and W. Xu (2022) Structured graph learning for scalable subspace clustering: from single view to multiview. IEEE Trans. Cybern. 52 (9), p. 8976–8986. External Links: Document Cited by: TABLE V, Appendix F, §I, §I, §I-A, §IV. [17] S. Wang, B. Tu, C. Xu, and Z. Zhang (2014) Exact subspace clustering in linear time. In Proc. AAAI Conf. Artif. Intell. (AAAI), Vol. 28, p. 2113–2120. External Links: Document Cited by: §I. [18] C. You, D. P. Robinson, and R. Vidal (2016) Scalable sparse subspace clustering by orthogonal matching pursuit. In Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 3918–3927. Cited by: TABLE V, Appendix F, §I, §I-A, §IV. [19] C. You, C. Li, D. P. Robinson, and R. Vidal (2016) Oracle based active set algorithm for scalable elastic net subspace clustering. In Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 3928–3937. Cited by: TABLE V, Appendix F, Appendix F, §I, §I-A, §IV. [20] Y. Chen, C. Li, and C. You (2020) Stochastic sparse subspace clustering. In Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 4155–4164. Cited by: TABLE V, Appendix F, §I, §I-A, §IV, §IV-C. [21] A. Adler, M. Elad, and Y. Hel-Or (2015) Linear-time subspace clustering via bipartite graph modeling. IEEE Trans. Neural Netw. Learn. Syst. 26 (10), p. 2234–2246. External Links: Document Cited by: §I. [22] D. Lim, R. Vidal, and B. D. Haeffele (2020) Doubly stochastic subspace clustering. Note: arXiv:2011.14859 Cited by: TABLE V, Appendix F, §I, §I-A, §IV. [23] S. Zhang, C. You, R. Vidal, and C. Li (2021) Learning a self-expressive network for subspace clustering. In Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 12393–12403. Cited by: §I, §I-B. [24] Z. Kang, H. Xu, B. Wang, H. Zhu, and Z. Xu (2019) Clustering with similarity preserving. Neurocomputing 365, p. 211–218. External Links: Document Cited by: §I. [25] P. Ji, T. Zhang, H. Li, M. Salzmann, and I. Reid (2017) Deep subspace clustering networks. In Adv. Neural Inf. Process. Syst. (NeurIPs), Vol. 30, p. 24–33. Cited by: §I, §I, §I-B. [26] J. Zhang, C. Li, C. You, X. Qi, H. Zhang, J. Guo, and Z. Lin (2019) Self-supervised convolutional subspace clustering network. In Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 5473–5482. Cited by: §I, §I, §I-B. [27] J. Lv, Z. Kang, X. Lu, and Z. Xu (2021) Pseudo-supervised deep subspace clustering. IEEE Trans. Image Process. 30, p. 5252–5263. External Links: Document Cited by: §I, §I-B. [28] X. Meng, Z. Huang, W. He, X. Qi, R. Xiao, and C. Li (2025) Exploring a principled framework for deep subspace clustering. In Proc. Int. Conf. Learn. Represent. (ICLR), Cited by: TABLE V, Appendix F, Appendix F, §I, §I-B, §IV, §IV-B. [29] S. Matsushima and M. Brbic (2019) Selective sampling-based scalable sparse subspace clustering. In Adv. Neural Inf. Process. Syst. (NeurIPs), Vol. 32. Cited by: TABLE V, Appendix F, Appendix F, §I-A, §IV. [30] C. Fettal, L. Labiod, and M. Nadif (2023) Scalable attributed-graph subspace clustering. In Proc. AAAI Conf. Artif. Intell. (AAAI), Vol. 37, p. 7559–7567. Cited by: §I-A. [31] D. Arthur and S. Vassilvitskii (2007) k-means++: the advantages of careful seeding. In Proc. 18th Annu. ACM-SIAM Symp. Discrete Algorithms (SODA), p. 1027–1035. Cited by: Appendix E, §I-D. [32] D. P. Kingma and J. Ba (2015) Adam: a method for stochastic optimization. In Proc. Int. Conf. Learn. Represent. (ICLR), Cited by: Appendix E, §I-D. [33] A. Y. Ng, M. I. Jordan, and Y. Weiss (2001) On spectral clustering: analysis and an algorithm. In Adv. Neural Inf. Process. Syst. (NeurIPs), Vol. 14, p. 849–856. Cited by: §I-F. [34] Y. Chen, C. Li, and C. You (2020) Stochastic sparse subspace clustering. In Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR), p. 4154–4163. Cited by: §IV. [35] N. Mrabah, M. Bouguessa, and S. Sami (2025) Scalable deep subspace clustering network. In IEEE International Conference on Data Science and Advanced Analytics (DSAA), Cited by: §IV-C. Appendix A Spectral Analysis of the LSR Approximation The normal equation of (2) is (X⊤X+2λIn)C=X⊤X.(X X+2λ I_n)C=X X. (38) The coefficient matrix is positive definite for λ>0λ>0, so the solution is unique. Since X⊤X=Vxdiag(β1,…,βr)Vx⊤X X=V_xdiag( _1,…, _r)V_x , the vectors in VxV_x diagonalize the solution on the row space of X. Along its jjth eigenvector direction, the normal equation becomes (βj+2λ)τj=βj( _j+2λ) _j= _j. In the orthogonal complement, the solution has an eigenvalue of zero. This proves Proposition 1. The LSR solution and CrowC_row have the same eigenvectors on the row space. The invariance of the Frobenius norm under orthonormal changes of coordinates therefore gives ‖CLSR⋆−Crow‖F2=∑j=1r(1−τj)2,\|C _LSR-C_row\|_F^2= _j=1^r(1- _j)^2, (39) which proves (6). For the rank-m projector Vx,mVx,m⊤V_x,mV_x,m , the first m eigenvalues are one and the remaining eigenvalues are zero. Splitting the squared difference over these two sets proves (7). These identities quantify the spectral approximation used to motivate the factor model. Appendix B Proofs for the Factor Model For fixed Z∈ℝdz×nZ ^d_z× n and feasible P∈ℝn×mP ^n× m, decompose the residual as Z−LP⊤=Z(In−PP⊤)+(ZP−L)P⊤.Z-LP =Z(I_n-P )+(ZP-L)P . (40) The two terms are orthogonal in the Frobenius inner product because (In−PP⊤)P=0(I_n-P )P=0. Moreover, ‖(ZP−L)P⊤‖F2=‖ZP−L‖F2\|(ZP-L)P \|_F^2=\|ZP-L\|_F^2 since P⊤P=ImP P=I_m. Taking squared norms proves (14). The second term in that identity is strictly convex in L. Thus, L=ZPL=ZP is its unique minimizer. The equality LP⊤=ZPP⊤LP =ZPP holds after this update. For fixed L, expansion gives (20). Let Z⊤L=UPΣPVP⊤Z L=U_P _PV_P have the dimensions specified in (21). The trace inequality yields tr(P⊤Z⊤L)≤tr(ΣP).tr(P Z L) ( _P). (41) The feasible matrix P=UPVP⊤P=U_PV_P attains equality. This proves the global optimality of the Procrustes update. When Z⊤LZ L is rank deficient, zero singular directions must be included so that both singular bases have m columns. Different orthonormal completions can yield different minimizers. To prove Proposition 2, eliminate L and expand the remaining objective: ‖Z−ZPP⊤‖F2=‖Z‖F2−tr(P⊤Z⊤ZP).\|Z-ZPP \|_F^2=\|Z\|_F^2-tr(P Z ZP). (42) The variational characterization of eigenvalues maximizes the trace over an m-dimensional dominant eigenspace of Z⊤Z Z. Its eigenvalues are σ12,…,σrz2 _1^2,…, _r_z^2, followed by zeros. Subtracting their largest m values from ‖Z‖F2\|Z\|_F^2 gives (16). When m<rzm<r_z, a spectral tie at the truncation boundary can make the optimal subspace nonunique. When m=rzm=r_z, the minimum residual is zero. Hence, every global minimizer satisfies ZPP⊤=ZZPP =Z, which implies range(Z⊤)⊆range(P)range(Z ) (P). Both spaces have dimension rzr_z, so they coincide. Consequently, C=PP⊤=VzVz⊤=Z†ZC=P =V_zV_z =Z Z is the unique row-space projector, where † denotes the Moore–Penrose pseudoinverse. Its orthonormal basis P remains nonunique. If m>rzm>r_z, the same inclusion holds, but range(P)range(P) also contains m−rzm-r_z directions orthogonal to the row space. This proves the remaining assertions of Proposition 2. For Proposition 3, each exact block update minimizes the same fixed-Z objective with its other block held constant. Neither update can increase that objective. Its successive values are nonnegative and nonincreasing, so they converge. Appendix C Subspace Preservation and Connectivity Reorder the samples by subspace and write Z=[Z1,…,Zk]Z=[Z_1,…,Z_k], where Za∈ℝdz×naZ_a ^d_z× n_a, ∑a=1kna=n _a=1^kn_a=n, and span(Za)=Saspan(Z_a)=S_a. Under Proposition 4, SaS_a has dimension ra>0r_a>0 and the subspaces form a direct sum. Choose a basis matrix Ωa∈ℝdz×ra _a ^d_z× r_a and write Za=ΩaΞaZ_a= _a _a, where Ξa∈ℝra×na _a ^r_a× n_a has full row rank. Define Ω =[Ω1,…,Ωk]∈ℝdz×rz, =[ _1,…, _k] ^d_z× r_z, (43) Ξ =blkdiag(Ξ1,…,Ξk)∈ℝrz×n. =blkdiag( _1,…, _k) ^r_z× n. The direct-sum condition makes Ω full column rank. Therefore, Z=ΩΞZ= and Ξ=Ω†Z = Z, which imply row(Z)=row(Ξ)row(Z)=row( ). Their orthogonal row-space projectors coincide: Z†Z Z Z =ΆΠ= (44) =blkdiag(Z1†Z1,…,Zk†Zk). =blkdiag(Z_1 Z_1,…,Z_k Z_k). The second equality uses the block structure of Ξ and row(Ξa)=row(Za)row( _a)=row(Z_a). At a global fixed-Z optimum with m=rzm=r_z, Appendix B gives C=Z†ZC=Z Z. Equation (44) therefore proves subspace preservation. For connectivity, the stated sample condition can be written rank(Z)+rank(Zℬ)>rarank(Z_A)+rank(Z_B)>r_a (45) for every partition of a class’s indices into two nonempty sets A and ℬB. Here Z_A and ZℬZ_B contain the corresponding columns of ZaZ_a. Indeed, their spans sum to SaS_a. Equality of their rank sum with rar_a is equivalent to those spans forming a direct sum. Set Ca=Za†ZaC_a=Z_a Z_a. Suppose the graph defined by its nonzero off-diagonal entries is disconnected. There is a nontrivial partition ,ℬA,B for which CaC_a is block diagonal. Since row(Ca)=row(Za)row(C_a)=row(Z_a), this row space then splits into two spaces supported on the respective coordinate sets. Coordinate projection onto A has dimension rank(Z)rank(Z_A), which also equals the rank of the corresponding diagonal block of CaC_a. The analogous equality holds for ℬB. Adding the two block ranks gives rank(Z)+rank(Zℬ)=rarank(Z_A)+rank(Z_B)=r_a, contradicting (45). Conversely, equality in this rank sum gives a direct-sum partition. Applying (44) to its two groups makes CaC_a block diagonal. Thus, the sample condition is also necessary for connectivity of the full row-space projector graph. A singleton containing a nonzero sample has a connected one-vertex graph and satisfies the condition vacuously. Every diagonal entry of CaC_a is positive when all samples are nonzero. A zero diagonal in an orthogonal projector forces the corresponding row and column to vanish. The identity Za=ZaCaZ_a=Z_aC_a would then make the corresponding sample zero, a contradiction. Positive row rescaling and entrywise squaring preserve the off-diagonal support. The affinity specified in the method therefore has exactly k connected components under these assumptions. The eigenspace of eigenvalue one of its normalized affinity is spanned by the degree-weighted component indicators. Any orthonormal basis of this eigenspace differs from the normalized indicator basis by an orthogonal rotation. Row normalization maps every sample of one component to the same direction. Distinct components yield distinct directions. Consequently, the ideal embedding admits zero-distortion k-means assignments equal to the components, up to a label permutation. Appendix D Properties of the Quadratic Affinity and Spectral Reduction Affinity and degrees The definition of svecsvec gives svec(R)⊤svec(S)=tr(R⊤S)svec(R) svec(S)=tr(R S) for symmetric R,S∈ℝm×mR,S ^m× m. Applying this identity to the rank-one matrices pipi⊤p_ip_i and pjpj⊤p_jp_j proves (26). Hence A is positive semidefinite. Its entrywise nonnegativity follows from (24). Let ei∈ℝne_i ^n be the iith standard basis vector. Since C2=C^2=C, δi _i =ei⊤CC⊤ei=ei⊤Cei=‖pi‖22, =e_i C e_i=e_i Ce_i=\|p_i\|_2^2, (46) ∑iδi _i _i =tr(C)=m. =tr(C)=m. The diagonal entries of an orthogonal projector belong to [0,1][0,1], so 0≤δi≤10≤ _i≤ 1. Moreover, δi=0 _i=0 if and only if pi=0p_i=0, which makes the entire corresponding affinity row zero. Removing these rows, therefore, leaves the degrees of all active samples unchanged. Proof of Proposition 5 Let s=rank(F+)s=rank(F_+). Choose s rows of F+F_+ spanning its row space. The corresponding s rank-one matrices span all matrices pipi⊤p_ip_i . In particular, their span contains ∑i∈ℐ+pipi⊤=P⊤P=Im. _i _+p_ip_i =P P=I_m. (47) A linear combination of s rank-one matrices has rank at most s. Thus, m≤sm≤ s. The upper bound follows from the dimensions of F+F_+. Since Δ is positive definite, B and F+F_+ have the same rank. If Gv=νvGv=ν v for ν>0ν>0, then (BB⊤)(Bv)=B(Gv)=νBv(B )(Bv)=B(Gv)=ν Bv and ‖Bv‖22=ν‖v‖22\|Bv\|_2^2=ν\|v\|_2^2. Conversely, a positive eigenpair of BB⊤B maps to an eigenpair of G by multiplication by B⊤B . For orthonormal eigenvectors, va⊤B⊤Bvbνaνb=1a=b. v_a B Bv_b _a _b=1\a=b\. (48) This proves the correspondence, multiplicities, and orthonormality. Leading vector and spectral convention Let ω=svec(Im)ω=svec(I_m). We have F+⊤1n+=ωF_+ 1_n_+=ω, F+ω=(δi)i∈ℐ+F_+ω=( _i)_i _+, and ‖ω‖22=m\|ω\|_2^2=m. Therefore, Gω=F+⊤Δ−1F+ω=F+⊤1n+=ω.Gω=F_+ ^-1F_+ω=F_+ 1_n_+=ω. (49) This proves (30). Let A++∈ℝn+×n+A_++ ^n_+× n_+ be the principal submatrix of A indexed by ℐ+I_+. The normalized affinity BB⊤B is similar to the row-stochastic matrix Δ−1A++ ^-1A_++. Its eigenvalues have a magnitude of at most one. Positive semidefiniteness further restricts them to [0,1][0,1]. Thus, the retained eigenvalue is dominant. Locking this vector remains valid when the eigenvalue one is repeated. Every active row then contains the strictly positive coordinate δi/m _i/m, which justifies the row normalization in (31). Proposition 6 (Conditional recovery on the active graph). Suppose A++A_++ has exactly k connected components. The row-normalized embedding in (31), computed from exact eigenvectors, has one distinct row per component. An exact solution of k-means therefore recovers those components up to label permutation. Proof. Let Υ∈0,1n+×k ∈\0,1\^n_+× k be the component membership matrix. The eigenvalue-one eigenspace of BB⊤B has dimension k and an orthonormal basis Δ1/2Υ(Υ⊤ΔΥ)−1/2O,O⊤O=Ik. ^1/2 ( )^-1/2O, O O=I_k. (50) The diagonal matrix Υ⊤ΔΥ contains the positive component volumes. Within component a, each row of this basis is a positive scalar multiple of the aath row of O. Row normalization removes that scalar. The k resulting rows are distinct because the rows of O are orthonormal. Clustering these rows into k nonempty groups has minimum within-cluster sum of squares zero, attained exactly by the component partition up to permutation. ∎ Implicit reduced operator Using (25) and the Frobenius inner-product identity, for each original index i∈ℐ+i _+, the corresponding row of F+vF_+v equals pi⊤smat(v)pip_i smat(v)p_i. Multiplying by F+⊤Δ−1F_+ ^-1 yields (33). This proves the implicit operator identity. Appendix E Implementation Conventions The implementation specifies the network architecture, preprocessing, initialization, random seeds, Adam settings, and the iteration budgets in Algorithm 1. Each training update follows the stated order. When gradients are accumulated over sample chunks, their contributions in (19) are summed before the parameter update. Joint training starts with zero Adam moments. Sample-based initialization For k-means++ seeding of L, select the first index uniformly from 1,…,n\1,…,n\. Let S be the indices already selected. For each remaining index, compute its squared distance to the closest selected column, diseed=minj∈‖zi−zj‖22,Pr(i∣)=diseed∑j∉djseed.d_i^seed= _j \|z_i-z_j\|_2^2, (i )= d_i^seed _j d_j^seed. (51) Use this distribution until m distinct indices have been selected. When its denominator is zero, sample uniformly from the remaining indices. Put the selected columns into L in selection order. This is the squared-distance seeding rule [31]. Adam update For completeness, let wtw_t be the vector containing the current network parameters and let gtg_t be its gradient of the specified objective. These are vectorizations used only to define the optimizer. Let ata_t and btb_t be its first- and second-moment vectors. With a0=b0=0a_0=b_0=0, coefficient inputs 0≤ρ1,ρ2<10≤ _1, _2<1, step size ηt>0 _t>0, and ϵAdam>0 _Adam>0, one update is at a_t =ρ1at−1+(1−ρ1)gt, = _1a_t-1+(1- _1)g_t, (52) bt b_t =ρ2bt−1+(1−ρ2)(gt⊙gt), = _2b_t-1+(1- _2)(g_t g_t), a^t a_t =at/(1−ρ1t),b^t=bt/(1−ρ2t), =a_t/(1- _1^t), b_t=b_t/(1- _2^t), wt+1 w_t+1 =wt−ηta^t⊘(b^t+ϵAdam). =w_t- _t\, a_t ( b_t+ _Adam). Products ⊙ , divisions ⊘ , and the square root are elementwise. These equations specify Adam [32], without an added weight-decay penalty. For pretraining, gtg_t differentiates (13); for joint training, it differentiates (17) with fixed factors. The moment vectors have the same dimension as wtw_t. Factor accuracy The thin decomposition in (21) uses m orthonormal columns, completing the singular bases when zero singular values occur. After each factor update, we check the residual ‖P⊤P−Im‖F\|P P-I_m\|_F at the chosen working precision. Spectral accuracy and reproducibility The eigensolver retains the analytical mode in (30) and computes the remaining dominant modes on its orthogonal complement. For prescribed residual, orthogonality, and positive-eigenvalue tolerances εeig>0 _eig>0, ϵorth>0 _orth>0, and εrank>0 _rank>0, respectively, we accept computed eigenpairs satisfying max1≤j≤k‖Gvj−νjvj‖2 _1≤ j≤ k\|Gv_j- _jv_j\|_2 ≤εeig, ≤ _eig, (53) max1≤a,b≤k|va⊤vb−a=b| _1≤ a,b≤ k|v_a v_b-1\a=b\| ≤ϵorth, ≤ _orth, νk _k >εrank. > _rank. The iteration budget or working precision is increased as needed to meet these criteria. Eigenvector ordering and random seeds are recorded for reproducibility. When an eigenvalue is repeated at the truncation boundary, the solver fixes a reproducible orthonormal basis of the selected modes. The active set contains every strictly positive degree in (27). Rows with zero degree are assigned by (32). Final clustering Let ι(i) (i) denote the row of F+F_+ corresponding to the original sample index i∈ℐ+i _+. The same ordering is used in H and H H. With spectral centers μa∈ℝk _a ^k, the clustering objective and Lloyd updates are minγi,μa _\ _i\,\ _a\ ∑i∈ℐ+∥H^ι(i):⊤−μγi∥22, _i _+\| H_ (i): - _ _i\|_2^2, (54) γi _i ←argmin1≤a≤k∥H^ι(i):⊤−μa∥22, ← _1≤ a≤ k\| H_ (i): - _a\|_2^2, μa _a ←1|a|∑i∈aH^ι(i):⊤. ← 1|J_a| _i _a H_ (i): . Use Rkm≥1R_km≥ 1 recorded k-means++ restarts on these rows. Apply the seeding rule (51) to the spectral rows, selecting k seeds among the n+n_+ active samples. Each restart takes at most TkmT_km Lloyd steps. Ties between centers use the smallest center index. Repair an empty cluster before its mean update by moving the sample with the largest current squared residual from a cluster containing at least two samples, with ties resolved by the original sample index. Then recompute affected centers. Stop early when assignments do not change. Keep the restart with minimum within-cluster sum of squares. This procedure returns k nonempty groups for n+≥kn_+≥ k, as ensured by (29). Appendix F Complexity Comparison with Baselines Table V separates representation learning from spectral computation. The bounds count sequential arithmetic operations. The numbers of training passes and optimization iterations are held fixed to isolate sample-size scaling. Eigensolver budgets are retained explicitly. To compare dense, sparse, and implicit spectral operators on the same basis, define sp(dop,ω)=Teig(ω+dopseig)+Reig(dopseig2+seig3).C_sp(d_op,ω)=T_eig(ω+d_ops_eig)+R_eig(d_ops_eig^2+s_eig^3). (55) Here, dopd_op is the operator dimension and ω bounds the cost of one matrix–vector product. The solver-specific budgets TeigT_eig, seigs_eig, and ReigR_eig have the meanings given in Section I-G. For an explicitly stored sparse graph, ω is its number of nonzero entries. A dense graph instead has ω=O(n2)ω=O(n^2). Thus, partial eigensolvers retain the benefit of sparse affinities. The common final k-means cost, O(Rkm(Tkm+1)nk2)O(R_km(T_km+1)nk^2), is added once to the sum of the two time columns for the complete pipeline. TABLE V: Time and representation-storage complexity of the compared methods. All entries are inside O(⋅)O(·). The two time columns are additive, with the common final k-means cost added once. Fixed optimization budgets are suppressed. Spectral budgets are included through (55). The last column counts the principal coefficient, factor, and affinity matrices, not total working memory. Bounds are obtained from the original algorithms under the implementation conventions stated below. Method Deep Representation / affinity learning Spectral computation Matrix storage LMVSC [15] No ndxm+nm3nd_xm+nm^3 nm2+m3nm^2+m^3 nmnm SGL [16] No ndxm+nm3nd_xm+nm^3 nm2+m3nm^2+m^3a nmnm S5C [29] No n(dxm+cℓ1)n(d_xm+c_ _1) sp(n,nm)C_sp(n,nm) nmnm EnSC [19] No n2dx+ncENn^2d_x+nc_EN sp(n,ns)C_sp(n,ns) nsns SSC-OMP [18] No n2dxsn^2d_xs sp(n,ns)C_sp(n,ns) nsns A-DSSC [22] No C+chk+eactC_C+C_chk+e_act sp(n,eact)C_sp(n,e_act) ℳC+eactM_C+e_act SSSC [20] No Rdropρn2dxsR_dropρ n^2d_xs sp(n,edrop)C_sp(n,e_drop) edrope_drop LSR [9] No n2dx+minn3,ndx2+dx3n^2d_x+ \n^3,nd_x^2+d_x^3\b sp(n,n2)C_sp(n,n^2) n2n^2 PRO-DSC [28] Yes (n/nb)cPRO(nb)+ncnet+n2dz aligned &(n/n_b)c_PRO(n_b)\\[-2.0pt] & 9.24994pt+nc_net+n^2d_z aligned sp(n,n2)C_sp(n,n^2) n2n^2 LoomSC Yes ncnet+ndzm+nm2+m3nc_net+nd_zm+nm^2+m^3 nm+sp(q,nq)+nqk aligned &nm+C_sp(q,nq)\\[-2.0pt] & 9.24994pt+nqk aligned nm+dzmnm+d_zm For LoomSC, q=m(m+1)/2q=m(m+1)/2. The table reports spectral computation using the operator in (33), without assembling G. s bounds the nonzeros per sample for EnSC and SSC-OMP, and per dropout subproblem for SSSC. For SSSC, ρ is the expected retained dictionary fraction and RdropR_drop is the number of dropout dictionaries, giving edrop=minn2,nsRdrope_drop= \n^2,nsR_drop\. The SSSC row gives the expected sequential computational cost at a fixed consensus budget. cℓ1c_ _1 and cENc_EN are the costs of one restricted Lasso and elastic-net solve, respectively. The former uses at most m representatives; the latter depends on the active-set size. For A-DSSC, CC_C and ℳCM_C denote the chosen coefficient solver’s computational cost and matrix storage. chkC_chk includes support initialization and global checks, and eacte_act bounds the entries processed by its support-restricted transport solver. nbn_b is PRO-DSC’s training batch size, and cPROc_PRO is defined in (56). Anchor-method spectral costs use a direct thin SVD of the n×mn× m factor. aSGL updates this representation within its graph-learning loop. bLSR uses the cheaper of the sample-space and feature-space linear solves. Both routes explicitly form the coefficient matrix. PRO-DSC [28] trains on mini-batches and constructs its final affinity over the complete clustering set. Combining its regularizer costs with neural computation and matrix products gives the following per-batch computational-cost bound: cPRO(nb)= c_PRO(n_b)= nbcnet+nb2dz+minnb3,dz3 n_bc_net+n_b^2d_z+ \n_b^3,d_z^3\ (56) +sp(nb,nb2). +C_sp(n_b,n_b^2). For 1≤nb≤n1≤ n_b≤ n, the n/nbn/n_b multiplier counts batches in a full training pass. A fixed number of Sinkhorn iterations is absorbed into the matrix-product term. The separate n2dzn^2d_z term accounts for constructing the full-set affinity at inference. The per-sample network cost cnetc_net is evaluated for each method’s own architecture. The regression-based methods differ in how they reduce optimization cost. S5C [29] uses a selection budget of m and solves the final restricted regression for every sample. Its bound includes both selection and regression costs. For EnSC [19], the active-set solver is left explicit because the original algorithm permits different restricted solvers. A-DSSC [22] likewise permits different coefficient estimators. Each support-restricted transport evaluation is linear in eacte_act. A complete support check additionally costs O(n2)O(n^2) with stored coefficients, or O(dxn2)O(d_xn^2) when feature-space LSR entries are evaluated on demand. For LSR [22] with dx<nd_x<n, a feature-space solve reduces coefficient construction to O(ndx2+dx3+n2dx)O(nd_x^2+d_x^3+n^2d_x). Explicit assembly and storage of the full coefficient matrix remain quadratic in n. The matrix-storage column isolates the representation responsible for affinity scaling. Total working memory also includes the input data, spectral embedding, eigensolver workspace, and regression-solver workspace. Deep methods additionally retain latent features and network state. For LoomSC, these terms are specified in (37), with the corresponding eigensolver workspace. Quadratic feature rows can be generated as needed, so an n×qn× q feature matrix need not be stored. The comparison identifies the source of LoomSC’s scalability advantage. LMVSC [15], SGL [16], and S5C [29] already have linear sample-size dependence for fixed representative and solver budgets. LoomSC retains this dependence while jointly learning nonlinear features and the two self-expression factors. EnSC [19] and SSC-OMP [18] can reduce coefficient storage through sparsity, but their dictionary scans still compare each target with the full dataset. SSSC reduces the dictionary size by dropout, giving the explicit Rdropρn2R_dropρ n^2 dependence in the table [20]. PRO-DSC [28] instead improves the training stage through mini-batches, while its full-affinity spectral inference retains quadratic storage. LoomSC avoids explicitly constructing either n×n× n matrix. Its q-dimensional reduction is exact for the defined squared affinity. For fixed m and solver budgets, both its learned factors and its reduced spectral computation scale linearly with n. Thus, LoomSC combines deep representation learning with scalability through the final clustering stage.