Variable Projection: Separable Optimization

• Variable Projection (VarPro) is a computational approach that leverages separable structures to reduce dimensionality and enhance numerical efficiency.
• It eliminates linear variables by solving their conditional subproblem exactly, simplifying the optimization of the remaining nonlinear parameters.
• Modern VarPro techniques extend to nonsmooth, large-scale, and manifold settings, improving robustness in inverse problems and machine learning applications.

Variable Projection (VarPro) is a computational and inferential principle for solving parameter estimation and function learning problems where the unknowns exhibit a partially linear (separable) structure. In its classical form, VarPro eliminates or “projects out” a subset of variables (typically those entering linearly) by solving the corresponding conditional problem exactly, reducing the remaining task to an optimization or inference over nonlinear parameters only. This dimensionality reduction enhances numerical efficiency, stability, and, in recent advances, theoretical convergence and adaptation properties in high-dimensional settings. Modern VarPro methodologies encompass both convex and nonconvex models, structured low-rank settings, kernel and tensor algebra frameworks, as well as the design of scalable second-order and Riemannian algorithms.

1. Classical VarPro Formulation and Theoretical Foundations

The core setting for VarPro is the separable nonlinear least squares problem: $\min_{x, y} \ f(x, y)$ where $x$ (linear parameters) and $y$ (nonlinear parameters) enter the function $f$ in a manner permitting, for fixed $y$, an explicit or efficient solution for $x$ (often via linear least squares). The classical Golub–Pereyra framework performs the following sequence:

• For a given $y$, solve $\bar{x}(y) = \arg\min_x f(x, y)$ (analytically or by iterative methods).
• Substitute $\bar{x}(y)$ back, yielding the reduced function $F(y) = f(\bar{x}(y), y)$.
• Optimize $F(y)$ (using, e.g., Gauss–Newton) to obtain the solution.

The analytic derivatives for $F(y)$ inherit a key property: the chain rule takes a simplified form, as the stationarity of the inner problem implies

$\nabla_y F(y) = \nabla_y f(\bar{x}(y), y)$

without requiring explicit differentiation of $\bar{x}(y)$ (Leeuwen et al., 2016).

Recent theoretical studies analyze the local convergence of classical and approximated VarPro algorithms; notably, Kaufman’s simpler Jacobian update yields similar local convergence to Golub–Pereyra’s full analytic Jacobian when the residual is small, as shown by Taylor expansions and Lipschitz bounds (Chen et al., 21 Feb 2024). For large residuals, curvature can be poorly approximated, necessitating robustified Hessian corrections (see Section 3).

2. Extensions to Non-Smooth and Large-Scale Problems

Traditional VarPro assumes smoothness and closed-form projections. Modern applications demand:

• Non-smooth regularization: The introduction of ℓ₁, TV, or group norms (for sparsity, edge-preservation, or robust estimation) breaks smoothness. By strong convexity and suitable smoothing (using, e.g., majorization-minimization or dualization), the projected function retains differentiability and the minimality over the projected variable is unique. Proximal and projected gradient algorithms then exploit these properties, and error inexactness (from iterative solvers) can be adaptively controlled to preserve overall convergence (Leeuwen et al., 2016).
• Large-scale inverse and imaging problems: In very high-dimensional settings, computing Jacobians and Hessians exactly is prohibitive. Iterative techniques (e.g., LSQR) are applied for inexact inner projections, with stopping criteria guided by bounds on residual and step size (Español et al., 13 Feb 2024). Error bounds demonstrate that with tolerances decreasing across outer iterations, the method converges at a comparable rate to exact Gauss–Newton VarPro.

VarPro has also been extended to nonsmooth over-parameterized problems, where the original non-smooth objective is recast as a smooth but nonconvex function via redundant variables (e.g., Hadamard decomposition for group sparsity). The VarPro reduction in these lifted spaces both improves conditioning and renders algorithms robust to dimensionality (Poon et al., 2022).

3. Structured and Nonlinear VarPro: Low-Rank, Affine and Fully Coupled Problems

VarPro applies naturally to structured low-rank approximation (SLRA) and weighted low-rank approximation (WLRA):

• In affine SLRA, constraints are imposed so that linearly parameterized structures (e.g., mosaic Hankel matrices $S(p)=S_0+\sum p_k S_k$) are projected onto sets with rank at most $r$. Here, VarPro reformulates a high-dimensional constrained problem to a low-dimensional manifold optimization (often Grassmannian), with cost, gradient, and Hessian approximations efficiently leveraging block and Toeplitz structure for scalability (Usevich et al., 2012).
• In weighed low-rank matrix approximation (WLRA), the cost is inherently overparameterized and the Jacobian is always rank-deficient due to the non-uniqueness of low-rank factors. Explicit formulas are developed for the gradient and Hessian; updating only along “active” directions (those changing the column space) turns VarPro into a Riemannian optimization on the Grassmann manifold, closely linking it to geometric approaches (Terray, 6 May 2025).
• For fully nonlinear, coupled problems (e.g., in conformal surface flattening or bilinear PDEs), a nonlinear extension of VarPro is employed. Here, neither block can be eliminated in closed form, but a nested inner-outer iterative scheme (iterative projection) solves the “more linear” block by repeated nonlinear minimization, then updates the coupled block. While computationally intensive, this nonlinear VarPro extension yields higher-quality solutions in problems with strong coupling restrictions (Miki, 26 Feb 2025).

4. Manifold, Kernel, and Tensor Generalizations

Emergent advances generalize VarPro to optimization on non-Euclidean spaces:

• Manifold-based VarPro: Many modern problems require optimizing over spaces constrained by orthogonality, subspace, or invariance properties (e.g., Grassmannians or O(n)). For example, polynomial ridge approximations and tensor decompositions often decouple subspace (U) or transform (M) parameters from linear coefficients; VarPro solves for the coefficients exactly, then applies Gauss–Newton or Riemannian gradient steps on the manifold parameter. Efficient update rules, geodesic-based step size selection, and orthogonality-adapted gradients enable scalable algorithms for problems including data-driven surrogate modeling and optimal matrix-mimetic tensor algebras (Hokanson et al., 2017, Newman et al., 11 Jun 2024).
• Kernel methods and projection operators: In scalable variable selection for multi-view learning, VarPro-inspired strategies employ projection operators to measure high-dimensional (possibly nonlinear) correlations and iteratively select features by projecting out redundancy. Updates leveraging SVDs and kernel expansions ensure low per-iteration complexity even at very large scale (Szedmak et al., 2023).

5. Recent Algorithmic Innovations: Large Residuals, Robustness, Initialization-Free and Deep Learning

Recent research addresses robustness, large residuals, and integration with modern learning architectures:

• Large residual problems: In separable optimization, significant residuals can distort the Gauss–Newton Hessian. The VPLR algorithm (Variable Projection for Large Residuals) incorporates an adaptive correction term, updated via a secant-like condition akin to quasi-Newton updates, to compensate for the neglected second-order residual coupling, demonstrated to yield improved convergence in challenging settings (Chen et al., 21 Feb 2024).
• Initialization-free and large-scale optimization: In bundle adjustment for structure-from-motion, power series expansions applied within a VarPro framework (PoVar) avoid full matrix inversion, perform updates efficiently, and exhibit a wide basin of convergence even from arbitrary initializations. By further incorporating Riemannian optimization in homogeneous coordinates, such methods scale to thousands of cameras and points without prior initialization (Weber et al., 8 May 2024).
• Neural networks and PDEs: VarPro principles pervade modern machine learning applications. Feature-extracting VP layers in neural networks enable compact, interpretable models for signal processing, outperforming standard architectures on tasks such as ECG classification (Kovács et al., 2020). In deep DNN training, Gauss-Newton VarPro methods extend classical projection to nonquadratic (e.g., cross-entropy) losses, projecting out the final affine layer and optimizing feature-extracting layers more efficiently than SGD or ADAM (Newman et al., 2020). In regression and inverse PDEs, VarPro paired with neural network representations or randomized networks yields efficient, high-accuracy, and flexible solvers (Dong et al., 2022, Dong et al., 2022).

6. Statistical and Model-Independent Variable Selection via VarPro

VarPro also underpins entirely model-independent variable selection methodologies:

• The “Variable Priority” (VarPro) approach evaluates the importance of features not by retraining predictors or simulating knockoffs, but by comparing local sample averages within rule-defined partitions of the data (e.g., leaves of a decision tree) to those with feature-specific constraints released. This approach is shown, under minimal and general regularity assumptions, to possess a consistent filtering property for noise variables: the importance score asymptotically vanishes if the target function is conditionally independent of the candidate variable set. Empirical studies document its robustness to high-dimensionality and strong feature correlations, with competitive or superior geometric mean (gmean) and area-under-curve (AUC) measures compared to lasso, permutation importance, and other state-of-the-art methods (Lu et al., 13 Sep 2024).

7. Impact, Limitations, and Future Directions

Variable Projection has driven advances in diverse fields: numerical linear algebra, inverse problems, high-dimensional statistics, machine learning, signal processing, tensor analysis, computer graphics, and PDE-constrained optimization. Its ability to reduce problem dimension, improve conditioning, and provide explicit analytic derivatives has enabled the development of scalable, robust, and high-quality algorithms for classical and emerging problems.

Principal limitations include computational cost in fully nonlinear or strongly coupled settings (where no block can be eliminated in closed form), intrinsic non-uniqueness (due to model overparameterization or invariance), and, for certain structured problems (e.g., with missing data or large-scale weighting), the need for regularization or careful algorithmic design to ensure numerical stability and differentiability.

Future directions are likely to focus on hybridizing VarPro with stochastic, manifold, and meta-learning approaches; extending analytic and convergence guarantees to more general nonconvex and over-parameterized regimes; deep integration into scientific machine learning workflows; and further leveraging its ability to provide theoretically justified and computationally efficient methods for non-smooth, large-scale, and high-dimensional optimization.