Hybrid LS/Gradient Descent Methods

• Hybrid least squares/gradient descent methods are optimization schemes that integrate closed-form least-squares insights with iterative gradient updates.
• They employ techniques such as SVD preconditioning, recursive least-squares, and orthogonal projections to enhance convergence and stability.
• These methods are scalable and versatile, with applications ranging from deep learning operator training to solving large-scale inverse problems.

Searching arXiv for the cited papers to ground the article in current literature. Searching (Chang et al., 7 Feb 2026) A hybrid least squares/gradient descent method is an optimization scheme that combines least-squares structure—residual maps, normal equations, Jacobians, Gauss–Newton matrices, recursive least-squares recursions, or exact least-squares subsolves—with gradient-based iteration, momentum, adaptive steps, or stochastic-gradient updates. In recent arXiv work, the term covers several distinct but mathematically related designs: SVD-preconditioned descent for nonlinear least squares, one-pass algorithms bridging recursive least-squares and orthogonal gradient descent, randomized conjugate-gradient least-squares, Krylov–regularization hybrids, and alternating least-squares/gradient training for operator networks (Chang et al., 7 Feb 2026, Min et al., 2022, Zeng et al., 24 May 2026, Yang, 2024, Choi et al., 21 Aug 2025, Choi et al., 8 Jul 2026).

1. Structural principle

The common template starts from a least-squares objective and then modifies a gradient-based update by injecting problem-specific quadratic structure. In nonlinear least squares one writes

$$\min_{\theta\in\mathbb{R}^m} f(\theta):=\frac12\|F(\theta)\|_2^2,$$

with gradient and Hessian

$$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$

Classical Gauss–Newton and Levenberg–Marquardt exploit $J_F^\top J_F$, whereas plain gradient descent uses only $J_F^\top F$. Hybrid methods typically keep an update of the form

$$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$

but choose $P_t$ from least-squares structure rather than from a generic first-order rule (Chang et al., 7 Feb 2026).

The same pattern appears outside standard nonlinear least squares. In stochastic least-squares regression, accelerated schemes retain stochastic-gradient updates but tune them to the covariance structure of the quadratic loss, improving the bias decay while preserving the minimax-optimal variance term $O(d\sigma^2/n)$ (Jain et al., 2017). In distributed linear least squares, iterative pre-conditioning learns an approximate inverse of $A^\top A$ and then uses it inside a gradient-descent outer loop (Chakrabarti et al., 2020). In operator learning, the last layer of DeepONets and MIONets is linear in its parameters, so those parameters can be updated by exact or alternating least-squares solves while the hidden parameters are updated by Adam or gradient descent (Choi et al., 21 Aug 2025, Choi et al., 8 Jul 2026).

A plausible unifying description is that hybrid methods partition the optimization variables or directions into components that admit closed-form least-squares treatment and components that remain better handled by iterative descent.

2. Principal algorithmic families

The recent literature organizes naturally into several families. They differ in where the least-squares structure enters—preconditioning, projection, variance reduction, or alternating exact solves—but all preserve an iterative descent backbone.

Family	Least-squares ingredient	Gradient ingredient
Spectral/preconditioned descent	$J_F^\top J_F$, Fisher-type matrices, or $A^\top A$ define preconditioners	Preconditioned GD, Adam, or distributed GD
Orthogonality/RLS hybrids	Recursive LS or minimum-norm interpolation constraints	Orthogonal gradient projection or stochastic-gradient updates
Krylov/projection hybrids	CGLS/LSMR projections, sketched normal equations, inner LSQR solves	Conjugate-gradient-style or projected descent
Alternating LS/GD in operator nets	Exact or alternating LS solves for last-layer parameters	Adam/GD for hidden branch and trunk parameters
Reduced-objective LS/GD methods	Linear parameters eliminated by LS, leaving nonlinear variables	Gradient descent on reduced objective

Representative instances are "SVD-Preconditioned Gradient Descent Method for Solving Nonlinear Least Squares Problems" (Chang et al., 7 Feb 2026), "One-Pass Learning via Bridging Orthogonal Gradient Descent and Recursive Least-Squares" (Min et al., 2022), "Randomized conjugate gradient least squares" (Zeng et al., 24 May 2026), "Hybrid LSMR algorithms for large-scale general-form regularization" (Yang, 2024), "A Novel Gradient Descent Least Squares (GDLS) Algorithm for Efficient SMV Gridless Line Spectrum Estimation with Applications in Tomographic SAR Imaging" (Shi et al., 2022), "Hybrid Least Squares/Gradient Descent Methods for DeepONets" (Choi et al., 21 Aug 2025), and "Hybrid Least Squares/Gradient Descent Methods for MIONets" (Choi et al., 8 Jul 2026).

3. Canonical mechanisms

A first mechanism is spectral preconditioning of the gradient. In SPGD, the Jacobian SVD $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$0 is used to define

$$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$1

yielding

$$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$2

The Adam variant treats $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$3 as the effective gradient and then applies first- and second-moment adaptation; for cross-entropy the paper replaces $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$4 by $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$5, where $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$6 is built from softmax Fisher blocks (Chang et al., 7 Feb 2026).

A second mechanism is recursive least-squares geometry enforced through gradient projections. ORFit projects the current loss gradient onto the orthogonal complement of the span of past sensitive directions, then chooses the step size so that the new point is interpolated in one step. For overparameterized linear models it solves

$$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$7

and is equivalent to EW-RLS with $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$8, $$\nabla f(\theta)=J_F(\theta)^\top F(\theta),\qquad \nabla^2 f(\theta)=J_F(\theta)^\top J_F(\theta)+\sum_{i=1}^n F_i(\theta)\nabla^2F_i(\theta).$$9. This connects minimum-norm interpolation, recursive least squares, and orthogonal gradient descent in a single update rule (Min et al., 2022).

A third mechanism is constraint-corrected gradient construction. RCGLS starts from the least-squares problem

$J_F^\top J_F$0

and interprets classical CGLS as the solution of a constrained correction problem that keeps the next search direction close to the negative gradient while enforcing $J_F^\top J_F$1-conjugacy. The randomized version replaces the full gradient by a sketched gradient $J_F^\top J_F$2, then adds the conjugate correction term $J_F^\top J_F$3. The result is gradient-related, but its geometry is explicitly inherited from the least-squares normal equations (Zeng et al., 24 May 2026).

A fourth mechanism is exact least-squares elimination of linear parameters. GDLS for gridless line spectrum estimation writes

$J_F^\top J_F$4

eliminates amplitudes by

$J_F^\top J_F$5

and then minimizes the reduced least-squares objective

$J_F^\top J_F$6

by gradient descent over the continuous frequencies (Shi et al., 2022).

A fifth mechanism is alternating least-squares inside neural-operator training. DeepONets are linear in the last-layer parameters of the branch network, so those parameters solve a ridge-type LS problem while the hidden branch and trunk parameters are updated by Adam. MIONets generalize this by multilinearity: each branch last layer is optimized in turn by ALS, with the other branch last layers frozen, and hidden/trunk parameters are optimized by gradient descent between ALS sweeps (Choi et al., 21 Aug 2025, Choi et al., 8 Jul 2026).

4. Convergence and acceleration properties

The strongest analytical results arise when the least-squares structure regularizes the difficult directions of the gradient flow. For SPGD, local linear convergence is proved under standard regularity assumptions. The comparison with plain GD is explicit: near a zero-residual solution, GD has a factor governed by

$J_F^\top J_F$7

whereas SPGD has a factor of the form

$J_F^\top J_F$8

so the degradation depends linearly rather than quadratically on the Jacobian condition number. A modified AMSGrad-style SPGD–Adam with

$J_F^\top J_F$9

and monotone second moments is shown to have vanishing average gradient norm under boundedness and smoothness assumptions (Chang et al., 7 Feb 2026).

RCGLS converges linearly in expectation, with a contraction bound that is always at least as good as the corresponding randomized coordinate-descent bound because the factor $J_F^\top F$0 satisfies $J_F^\top F$1. The method also inherits a natural variance-reduction property from the sketched gradient $J_F^\top F$2: as $J_F^\top F$3, the estimator variance automatically vanishes (Zeng et al., 24 May 2026).

In stochastic least-squares regression, the accelerated stochastic-gradient method of "Accelerating Stochastic Gradient Descent For Least Squares Regression" achieves the minimax-optimal variance term $J_F^\top F$4 while improving the bias contraction rate from $J_F^\top F$5 to $J_F^\top F$6. This refutes the blanket claim that acceleration is unusable in noisy stochastic approximation for least squares (Jain et al., 2017).

In distributed linear least squares, iterative pre-conditioning converges linearly in general and achieves superlinear convergence when the least-squares problem has a unique solution. The analysis compares favorably with traditional GD and with the accelerated methods discussed in that paper (Chakrabarti et al., 2020).

ORFit supplies a different convergence statement: in overparameterized linear models, its iterate is exactly the point to which sufficiently small-step-size multi-pass SGD would converge, and the nonlinear NTK analysis extends the same minimum-norm interpretation to highly overparameterized networks (Min et al., 2022).

5. Computational realizations and scalability

A recurrent theme is that hybridization is only useful if the least-squares substructure can be exploited without destroying scalability. SPGD addresses the cost of $J_F^\top F$7 by approximating $J_F^\top F$8 with Lanczos iteration in Krylov subspaces

$J_F^\top F$9

using only Jacobian-vector and vector-Jacobian products. This avoids explicit formation of $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$0 and makes the method viable for models with hundreds of thousands of parameters (Chang et al., 7 Feb 2026).

ORFit uses incremental PCA to compress the basis of past sensitive directions. The memory footprint drops from $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$1 to $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$2, with per-step overhead $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$3, while still approximating the orthogonality constraints required for one-pass learning (Min et al., 2022).

RCGLS develops an equivalent sparse implementation that avoids full-dimensional vector operations by maintaining implicit representations of the iterates and applying $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$4 or $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$5 only to sketched pieces such as $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$6 and $$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$7. This is central for sparse matrices, where vanilla RCGLS would otherwise be impractical (Zeng et al., 24 May 2026).

Hybrid LSMR introduces a two-level architecture. The outer Golub–Kahan/LSMR step constructs a rank-$$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$8 projected normal-equation problem; the inner correction solves

$$\theta_{t+1}=\theta_t-\alpha_t P_t \nabla f(\theta_t),$$9

with LSQR. The paper proves that the condition number of $P_t$0 decreases with $P_t$1, and states that it is generally enough to set $P_t$2 in LSQR (Yang, 2024).

DeepONet and MIONet hybrids turn an infeasible gigantic least-squares system into small structured problems. DeepONet uses Kronecker products and a generalized Sylvester equation for the last-layer solve; MIONet generalizes this with Khatri–Rao products and tensor permutation matrices, then solves one branch last layer at a time by ALS. The practical consequence is that LS updates can be inserted into training at modest cost relative to Adam epochs (Choi et al., 21 Aug 2025, Choi et al., 8 Jul 2026).

6. Applications, limitations, and interpretation

Hybrid least squares/gradient descent methods have been deployed across markedly different problem classes. SPGD is tested on function approximation, PDE solving, and image classification on the CIFAR-10 dataset (Chang et al., 7 Feb 2026). ORFit is evaluated on a Rotated MNIST streaming-regression task and on overparameterized linear and NTK-linearized settings (Min et al., 2022). RCGLS extends naturally to ridge regression and sparse large-scale linear problems (Zeng et al., 24 May 2026). Hybrid LSMR targets ill-posed inverse problems with general-form regularization (Yang, 2024). GDLS is motivated by gridless line spectrum estimation and TomoSAR imaging (Shi et al., 2022). DeepONet and MIONet LSGD methods target supervised and physics-informed operator learning for advection, diffusion–reaction, and Poisson-type equations (Choi et al., 21 Aug 2025, Choi et al., 8 Jul 2026). The Riemannian LS inverse eigenvalue method interprets lift-and-projection as gradient descent in an induced metric and includes an application to a $P_t$3 sparse Mn$P_t$4 Hamiltonian (Riley et al., 10 Apr 2025).

Several misconceptions are corrected by this literature. First, a hybrid least-squares/gradient method is not synonymous with “run Gauss–Newton instead of gradient descent.” Some hybrids use a square-root inverse of $P_t$5 rather than a full inverse, some insert least-squares solves only on a subset of parameters, and some use least-squares only to build control variates or memory summaries (Chang et al., 7 Feb 2026, Nobile et al., 28 Jul 2025). Second, the least-squares component need not appear as explicit normal equations; it may appear as a projection constraint, a recursive filter, a Sylvester equation, or a reduced objective after eliminating linear variables (Min et al., 2022, Choi et al., 21 Aug 2025, Shi et al., 2022). Third, the acceleration is not unconditional. In ill-conditioned Bratu regimes, the extrapolation-accelerated preconditioned gradient methods studied in "A New Combination of Preconditioned Gradient Descent Methods and Vector Extrapolation Techniques for Nonlinear Least-Squares Problems" did not consistently improve convergence, whereas restarted GNKS remained robust (Mouhssine, 31 Jan 2026).

The main limitations are structural. SPGD requires Jacobian-based spectral approximations and uses regularization and AMSGrad-style modifications to avoid Adam-type non-convergence issues (Chang et al., 7 Feb 2026). ORFit’s IPCA compression inevitably loses some information and can induce forgetting (Min et al., 2022). Hybrid LSMR reduces to plain LSMR when $P_t$6, so its distinctive regularization role appears only for $P_t$7 (Yang, 2024). DeepONet LSGD depends on linearity of the last branch layer and on operator factorizations that preserve the LS structure; MIONet LSGD further assumes product-form data layouts and tensor factorizability (Choi et al., 21 Aug 2025, Choi et al., 8 Jul 2026).

Taken together, these methods show that least-squares structure can be exploited at many granularities: as curvature information, as exact subproblem solvability, as recursive memory, as variance-reduction regression, or as tensorized linear algebra inside neural operators. The recurring mathematical idea is to extract the parts of the model for which least-squares geometry is informative or tractable, and to embed that information inside a gradient-based outer iteration rather than replacing descent altogether.