Subspace Regularization: Methods & Insights

• Subspace Regularization (SSR) is a family of methods that leverage inherent low-dimensional structures to regularize computations across optimization, inverse problems, and representation learning.
• SSR techniques include projection methods, subspace recycling, and explicit priors on subspace coefficients, each tailored to specific problem domains for improved efficiency and accuracy.
• Practical applications of SSR span unconstrained optimization, hyperspectral image super-resolution, and streaming 3D reconstruction, demonstrating significant gains in convergence and computational performance.

Searching arXiv for papers and exact SSR variants to ground the article. arxiv_search(query="subspace regularization SSR arXiv", max_results=10) arxiv_search(query="all:\"Subspace projection regularization\" OR all:\"Subspace Recycling-based Regularization\" OR all:\"New subspace minimization conjugate gradient methods based on regularization model\"", max_results=10) arxiv_search(query="all:\"A preconditioned Krylov subspace method for linear inverse problems with general-form Tikhonov regularization\" OR all:\"Hybrid LSMR algorithms for large-scale general-form regularization\"", max_results=10) Subspace Regularization (SSR) denotes a family of methods that impose, exploit, or preserve low-dimensional subspace structure as part of regularization, search-direction construction, or latent-state correction. The term is not standardized across the literature: it has been used for $p$-regularized two-dimensional subspace minimization in nonlinear conjugate-gradient optimization, Krylov and projection schemes for Bayesian and general-form inverse problems, nonconvex low-rank/sparse representation learning for subspace clustering, a training-free self-expressive operator for streaming 3D reconstruction, and subspace-coefficient tensor regularization for hyperspectral image super-resolution (Zhao et al., 2020, Li, 2023, Brbić et al., 2018, Deng et al., 16 Mar 2026, Zhang et al., 5 Aug 2025). Across these uses, SSR typically acts by restricting computation to a low-dimensional subspace, regularizing coefficients defined in that subspace, or projecting evolving estimates back toward a subspace-consistent trajectory.

1. Terminological scope and unifying structure

The acronym “SSR” is best understood as an umbrella label rather than a single algorithm. In different domains, the “subspace” may be a search plane $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$, a Krylov solution space, a union-of-subspaces representation space, a Grassmannian latent-state manifold, or a spectral coefficient subspace for tensor reconstruction (Zhao et al., 2020, Li, 2023, Brbić et al., 2018, Deng et al., 16 Mar 2026, Zhang et al., 5 Aug 2025).

Context	Regularized object	Mechanism
Unconstrained optimization	search direction	$p$-regularized subproblem in a two-dimensional subspace
Linear inverse problems	approximate solution	projection onto Krylov or augmented subspaces
Subspace clustering	representation matrix $C$	low-rank and sparse penalties
Streaming 3D reconstruction	latent recurrent state	self-expressive affinity correction
HSI super-resolution	subspace coefficient tensor $C$	low-rankness and smoothness on coefficients

This suggests a common operational pattern: regularization is achieved not only by adding penalties to a full-space objective, but also by choosing a computational subspace whose geometry encodes prior structure. A common misconception is that SSR always refers to an explicit penalty term. In fact, some SSR methods are projection procedures with early stopping, some are augmented iterative schemes, and some are training-free inference-time operators rather than optimization-based estimators.

2. $p$-regularized subspace minimization in unconstrained optimization

In nonlinear optimization, SSR appears as a subspace minimization conjugate-gradient framework based on the local model

$$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$

where $x_k\in\mathbb R^n$, $g_k=\nabla f(x_k)$, $B_k$ is a symmetric positive-definite approximation to $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$0, $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$1, and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$2 is chosen dynamically (Zhao et al., 2020). The model may also use a scaled norm $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$3 with $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$4, $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$5, and two choices are emphasized: $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$6 and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$7.

The defining construction is the restriction of $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$8 to the two-dimensional subspace

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$9

In this basis, the model becomes a $p$0 $p$1-regularized subproblem. The unique global minimizer $p$2 satisfies

$p$3

with $p$4 defined by the corresponding norm equation. For $p$5 and $p$6, explicit formulas are obtained for $p$7, $p$8, and $p$9. The resulting direction

$C$0

satisfies the sufficient descent condition

$C$1

using the condition that $C$2 is bounded away from zero.

Globalization is obtained by a modified nonmonotone Wolfe-type line search of Zhang–Hager: $C$3 with $C$4, and reference levels updated through

$C$5

Under assumptions that $C$6 is $C$7 on a level set containing $C$8, $C$9 is $C$0-Lipschitz, the nonmonotone Wolfe conditions hold, and $C$1 and $C$2 stay positive and bounded, one obtains $C$3, $C$4, and hence $C$5. If $C$6 is convex, satisfies the global error bound $C$7, and $C$8 are uniformly bounded above, then

$C$9

so the decay is $p$0-linear.

The numerical study uses 145 unconstrained problems from the CUTEr collection and reports performance in $p$1, $p$2, $p$3, and $p$4. Against CG_DESCENT 5.3, CGOPT, SMCG_BB, and SMCG_Conic, SSR with $p$5 solves $p$6 of problems, wins on $p$7 in $p$8, on $p$9 in $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$0, and is fastest on $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$1 in CPU time; for strongly ill-conditioned problems it often reduces $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$2 by $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$3–$$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$4 versus CG_DESCENT or CGOPT. In this literature, “subspace regularization” therefore refers to regularized search-direction generation inside a dynamically chosen low-dimensional plane rather than to parameter shrinkage in a statistical sense.

3. Projection-based SSR for Bayesian and general-form inverse problems

In inverse problems, SSR commonly denotes projection of a large ill-posed problem onto a sequence of solution subspaces, usually Krylov spaces generated by Golub–Kahan-type processes, with the iteration index acting as the regularization parameter or with additional regularization applied to the projected problem (Li, 2023, Li, 2023, Yang, 2024).

For Bayesian linear inverse problems with Gaussian noise $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$5 and Gaussian prior $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$6, the objective is

$$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$7

The method of subspace projection regularization (SPR) introduces the weighted Hilbert spaces $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$8 and $$m_k(s)\coloneqq g_k^T s+\frac12 s^T B_k s+\frac{\sigma_k}{p}\|s\|^p,$$9, with adjoint

$x_k\in\mathbb R^n$0

and builds the Krylov subspace

$x_k\in\mathbb R^n$1

via a generalized Golub–Kahan bidiagonalization. Writing $x_k\in\mathbb R^n$2, the projected problem becomes

$x_k\in\mathbb R^n$3

with $x_k\in\mathbb R^n$4. LSQR-type recurrences update $x_k\in\mathbb R^n$5, while discrepancy principle, L-curve, and generalized cross-validation select the stopping index. The filtered-GSVD expansion

$x_k\in\mathbb R^n$6

formalizes semi-convergence and the preference for dominant GSVD modes.

A related general-form Tikhonov formulation considers

$x_k\in\mathbb R^n$7

with $x_k\in\mathbb R^n$8. The preconditioned Golub–Kahan bidiagonalization (pGKB) method introduces

$x_k\in\mathbb R^n$9

works in the $g_k=\nabla f(x_k)$0-inner product, and generates a basis $g_k=\nabla f(x_k)$1 for a $g_k=\nabla f(x_k)$2-orthonormal solution subspace. Restricting $g_k=\nabla f(x_k)$3 yields the reduced problem

$g_k=\nabla f(x_k)$4

Two hybrid variants regularize the small problem using WGCV or a secant update driven by the discrepancy principle. The associated analysis shows a filtered-GSVD form with filter factors

$g_k=\nabla f(x_k)$5

and hence the expected semi-convergence behavior.

Hybrid LSMR places the same projection idea in the discrepancy-form problem

$g_k=\nabla f(x_k)$6

with $g_k=\nabla f(x_k)$7 from Golub–Kahan bidiagonalization of $g_k=\nabla f(x_k)$8. The projected Tikhonov problem is

$g_k=\nabla f(x_k)$9

and the discrepancy-form correction is

$B_k$0

The inner least-squares problem is solved by LSQR, and the condition number satisfies

$B_k$1

so conditioning improves with $B_k$2. The reported large-scale experiments indicate that the best regularized solution is as accurate as that of JBDQR, while wall-clock time is $B_k$3–$B_k$4 lower.

The empirical evidence emphasizes efficiency and early stopping. In the Bayesian SPR study, the gravity problem with $B_k$5 and noise $B_k$6 yields semi-convergence at $B_k$7 with relative error $B_k$8, versus a generalized-hybrid competitor reaching $B_k$9 at $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$00; in PRspherical with $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$01, $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$02, and noise $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$03, SPR gives $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$04 and error $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$05, while the comparison method reaches $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$06 at $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$07. At the same time, the literature also records limitations: GCV may fail by over-stopping in some cases, and hybrid projected methods can over-smooth when the projected parameter choice is not well matched to the noise model.

4. Recycling, decomposition, and restricted coercivity

A second inverse-problem line uses SSR to mean augmentation or decomposition of the solution space itself. In subspace recycling-based regularization, one starts with a finite-dimensional augmentation subspace $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$08, lets $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$09, and defines projectors $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$10 and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$11 with

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$12

The exact solution is split as

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$13

where the projected component is stably recovered from noisy data by

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$14

and the complementary component is reconstructed by applying any regularization $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$15 to the projected operator $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$16 and right-hand side $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$17 (Ramlau et al., 2020). Under the compatibility and projection assumptions, the augmented scheme

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$18

is itself a regularization whenever $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$19 regularizes the projected system. Augmented Landweber and augmented steepest descent then inherit discrepancy-principle stopping while often accelerating semi-convergence.

The same paper gives an explicit augmented gradient update: $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$20 with steepest-descent step

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$21

The numerical examples include a $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$22 Gaussian blur operator with additive uniform noise scaled so $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$23, adaptive optics wavefront reconstruction, and atmospheric image deconvolution. In the image deconvolution study, using the first $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$24 eigenvectors of the telescope’s diffraction-limited PSF operator accelerates residual and error reduction, with diminishing returns beyond $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$25.

A more abstract convex analysis of regularized least squares studies

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$26

through the decomposition

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$27

Writing $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$28, one can express the solution set via the conjugate $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$29. The range component is unique,

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$30

and the full solution set is

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$31

Existence, compactness, and uniqueness are then separated rather than conflated. Compactness is characterized by

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$32

while uniqueness holds iff

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$33

The notion of restricted coercivity is central: $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$34 is coercive on a subspace $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$35 if

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$36

equivalently if there exist $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$37 and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$38 such that

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$39

Applied to $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$40, restricted coercivity prevents escape along null directions and yields a precise interpretation of how regularization constrains the kernel slice of the solution set (Xue et al., 28 Jul 2025).

Together, these works broaden SSR beyond “solve a small projected problem.” They show that subspace regularization may mean preserving a useful recycled component, or decomposing the ambient space so that the range part and kernel part of the solution can be regularized and analyzed separately.

5. Representation-level SSR: self-expression, clustering, and Grassmannian state correction

In representation learning, SSR often regularizes a coefficient matrix or a latent trajectory by exploiting self-expressive subspace structure. Two prominent instances are low-rank sparse subspace clustering and training-free state correction for streaming 3D reconstruction (Brbić et al., 2018, Deng et al., 16 Mar 2026).

In low-rank sparse subspace clustering, the data matrix $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$41 is assumed to lie approximately in a union of unknown low-dimensional subspaces, and one seeks $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$42 satisfying

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$43

Conventional LRSSC solves

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$44

The SSR variants replace the nuclear and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$45 norms by tighter nonconvex regularizers. GMC-LRSSC uses the generalized minimax-concave penalty

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$46

and leads to firm-thresholding proximal maps. S$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$47-LRSSC uses the exact quasi-norms

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$48

with hard-thresholding proximal operators

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$49

The joint $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$50 proximal map is approximated by a proximal average, and both algorithms are solved by ADMM with convergence statements based on boundedness, KKT conditions, and semi-algebraicity. On synthetic data, S$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$51-LRSSC attains the lowest clustering error at low noise and large $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$52, while GMC-LRSSC is superior under higher noise. On Extended Yale B, GMC-LRSSC reports $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$53 and S$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$54-LRSSC reports $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$55 for $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$56, compared with $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$57. Each ADMM iteration costs $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$58, and both methods converge in $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$59 iterations.

Streaming 3D reconstruction uses SSR in a different sense. The latent persistent state $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$60 is modeled as an orthonormal basis of an $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$61-dimensional subspace, hence a point on the Grassmann manifold $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$62. Distances between subspaces are measured by

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$63

Over a local window, the sequence is assumed to satisfy the self-expressive relation

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$64

where $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$65. Rather than solving the nuclear-norm formulation at inference time, SSR forms a local window

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$66

computes normalized dot-product affinities

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$67

and corrects only the newest state by

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$68

The operator is training-free, plug-and-play, introduces no new learned weights, and has per-frame complexity $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$69 for similarity and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$70 for the final weighted sum, with $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$71 small, e.g. $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$72. On long synthetic Sintel sequences it reduces Absolute Rel from $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$73 and raises $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$74 from $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$75; on Bonn (dynamic), Abs Rel falls $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$76 and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$77 rises $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$78; on KITTI, Abs Rel drops $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$79 and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$80 rises $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$81. Camera-pose results include TUM-Dynamics ATE $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$82 m and ScanNet ATE $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$83 m. A notable caveat is that in sparse-view settings it underperforms slightly because affinities can become degenerate when $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$84 window length.

These two uses of SSR share the self-expressive prior, but differ sharply in implementation. In clustering, self-expression defines a nonconvex optimization problem over $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$85; in streaming reconstruction, it defines an analytic inference-time correction operator on a temporal state window.

6. Subspace-coefficient tensor regularization for hyperspectral image super-resolution

In hyperspectral image super-resolution, SSR is realized by regularizing subspace coefficients rather than the full hyperspectral tensor. The target HR-HSI $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$86 is represented as

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$87

where $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$88 is an orthonormal spectral basis with $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$89, $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$90 is the coefficient tensor, and $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$91 (Zhang et al., 5 Aug 2025). In matricized form,

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$92

Given LR-HSI $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$93 and co-registered HR-MSI $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$94, the degradations are

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$95

with spatial blur $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$96, down-sampler $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$97, and spectral response $\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$98. The data-fidelity term is

$\Omega_k=\operatorname{span}\{-g_k,d_{k-1}\}$99

Regularization is imposed on clustered coefficient tensors $p$00, obtained by grouping $p$01 MSI patches and transferring the same grouping to the $p$02-band coefficient patches. The unified tensor regularizer JLRST is

$p$03

where $p$04, $p$05, and $p$06 is the circulant difference of $p$07. The mode-3 logarithmic tensor nuclear norm is

$p$08

with $p$09. Its stated role is to penalize small singular values heavily while preserving large ones, thereby reducing the bias toward rank-overestimation inherent in convex TNN.

Optimization is performed by ADMM using auxiliary tensors $p$10 and $p$11. The $p$12-update yields the Sylvester equation

$p$13

the $p$14-update admits closed form via singular-value shrinkage in the Fourier domain, and the $p$15-update solves

$p$16

efficiently via FFT. Under boundedness and diminishing multiplier increments, every accumulation point satisfies the KKT conditions of the original constrained problem.

The experimental validation uses Pavia University, Indian Pines, CAVE “Balloons,” and University of Houston, with metrics PSNR, SSIM, ERGAS, SAM, and UIQI. With subspace dimension $p$17 and cluster count $p$18, JLRST attains the highest PSNR and SSIM, including $p$19 dB over the second-best method on Indian Pines and $p$20 dB on Balloons, while also obtaining the lowest ERGAS and SAM and the highest UIQI across three of four datasets. Ablation with any $p$21 degrades performance, the convergence curves stabilize within $p$22 iterations, and runtime is reported as moderate, including $p$23 s on Pavia.

This formulation makes explicit a recurring SSR design choice: rather than regularizing the full ambient object, one first constructs a low-dimensional subspace representation and then places structurally richer priors on the subspace coefficients. In this case, global low-rankness and local smoothness are enforced after spectral dimension reduction, which is both a modeling assumption and a computational device.

7. Recurring themes, distinctions, and limitations

Across the literatures surveyed here, SSR has three recurrent meanings. First, it can denote regularized minimization inside a low-dimensional search subspace, as in $p$24-regularized conjugate-gradient methods. Second, it can denote projection of an ill-posed problem onto a sequence of solution subspaces, where regularization is delivered by filtering, early stopping, or hybrid projected Tikhonov solves. Third, it can denote explicit priors on subspace coefficients or latent subspace trajectories, including self-expressive, low-rank, sparse, and smoothness-inducing constructions (Zhao et al., 2020, Li, 2023, Brbić et al., 2018, Deng et al., 16 Mar 2026, Zhang et al., 5 Aug 2025).

The main distinctions concern what is being regularized. In optimization, SSR regularizes the step computation. In inverse problems, it regularizes the solution reconstruction by restricting or augmenting the recoverable space. In clustering and representation learning, it regularizes a coefficient matrix or sequence affinity structure. In hyperspectral reconstruction, it regularizes coefficient tensors after spectral subspace factorization. Treating these as a single method would therefore be misleading.

The limitations are likewise domain-specific. Projection methods exhibit semi-convergence and depend critically on stopping criteria such as discrepancy principle, L-curve, or GCV; GCV can fail by over-stopping, and WGCV can over-smooth in some settings (Li, 2023, Li, 2023). Self-expressive streaming correction can underperform in sparse-view settings because affinities become degenerate when the temporal window is poorly matched to the sequence (Deng et al., 16 Mar 2026). Nonconvex clustering penalties offer tighter approximations to rank and sparsity, but their tractability relies on ADMM splitting, proximal approximations, and convergence conditions rather than on global convexity (Brbić et al., 2018). Subspace-factorized tensor models depend on the quality of the chosen spectral basis and patch grouping (Zhang et al., 5 Aug 2025).

The broad significance of SSR is therefore methodological rather than terminological. It identifies subspace structure as a regularization vehicle: one can regularize by selecting a subspace, by solving within it, by decomposing the ambient space into identifiable and null components, or by enforcing priors on subspace coefficients. The literature shows that these choices lead to distinct algorithmic families with different guarantees, failure modes, and computational advantages, even though they share the same name.