Variable projection
Variable projection is a reduction method for optimization problems in which some unknown parameters enter a model linearly while the remaining parameters enter nonlinearly. It is principally associated with separable nonlinear least squares, where eliminating the linear parameters produces a lower-dimensional objective function over the nonlinear parameters. The reduced problem has the same attainable residual norms as the original problem, provided that the eliminated linear subproblem is solved consistently.
The method was introduced by Gene H. Golub and Victor Pereyra in their 1973 analysis of differentiating pseudoinverses and solving nonlinear least-squares problems whose variables separate. Its name refers to the projection of the observed data onto the orthogonal complement of a parameter-dependent linear model space. Although the method removes variables from the explicit optimization, their influence remains encoded in the changing projection operator.
Mathematical formulation
Let (y\in\mathbb{R}^m) be an observation vector, let (\alpha\in\mathbb{R}^p) denote nonlinear parameters, and let (c\in\mathbb{R}^n) denote linear coefficients. A separable model has the form
[ f(\alpha,c)=\Phi(\alpha)c, ]
where the columns of the matrix (\Phi(\alpha)) are basis functions evaluated at the observation points. The associated least-squares problem is
[ \min_{\alpha,c}\left|y-\Phi(\alpha)c\right|_2^2. ]
For fixed (\alpha), the minimizing linear coefficient vector is
[ c(\alpha)=\Phi(\alpha)^\dagger y, ]
where (\Phi(\alpha)^\dagger) is the Moore–Penrose inverse. Substitution into the original residual yields
[ r(\alpha)
y-\Phi(\alpha)\Phi(\alpha)^\dagger y
P_{\Phi(\alpha)}^\perp y. ]
Here,
[ P_{\Phi(\alpha)}^\perp
I-\Phi(\alpha)\Phi(\alpha)^\dagger ]
is the orthogonal projector onto the complement of the column space of (\Phi(\alpha)). Variable projection therefore replaces the joint problem with the reduced problem
[ \min_\alpha \left|P_{\Phi(\alpha)}^\perp y\right|_2^2. ]
Every value of the reduced objective equals the minimum residual norm attainable by the original model at the corresponding nonlinear parameter value. Consequently, minimizers of the reduced problem generate minimizers of the full problem after reconstruction of (c(\alpha)), subject to the usual qualifications concerning nonuniqueness when (\Phi(\alpha)) is rank deficient.
The reduction does not amount to treating the linear coefficients as fixed. Instead, it incorporates their complete least-squares response to every change in (\alpha). This distinction separates variable projection from alternating optimization, in which the two parameter blocks are updated through successive conditional steps.
Differential structure
The reduced residual depends on (\alpha) through both the basis matrix and its pseudoinverse. When (\Phi(\alpha)) has constant column rank in a neighborhood and is differentiable, the differential of the residual can be written as
[ dr
-P_\Phi^\perp(d\Phi)c
(\Phi^\dagger)^{T}(d\Phi)^{T}r, ]
where (c=\Phi^\dagger y), (r=P_\Phi^\perp y), and the dependence on (\alpha) has been suppressed. Applying this identity to each component of (\alpha) produces the Jacobian matrix required by reduced Gauss–Newton methods.
The first term measures how the fitted model changes outside its current column space. The second term accounts for movement of that column space relative to the existing residual. Their geometric roles are distinct even though both arise from differentiating the same orthogonal projection.
Linda Kaufman examined a reduced Jacobian approximation in 1975 that omits the second term. The omitted term is proportional to the residual and therefore disappears at an exact fit. Near a small-residual solution, the approximation produces a Gauss–Newton model closely related to the exact reduced model, while away from such a solution the two curvature approximations can differ materially.
The gradient of the reduced objective
[ F(\alpha)=\frac{1}{2}|r(\alpha)|_2^2 ]
is
[ \nabla F(\alpha)=J(\alpha)^T r(\alpha), ]
with (J) denoting the reduced residual Jacobian. Because (r) is orthogonal to the columns of (\Phi), several terms generated by direct differentiation cancel in the objective gradient even though they remain relevant to the full residual Jacobian.
Weighted observations
A weighted least-squares problem has the form
[ \min_{\alpha,c} \left(y-\Phi(\alpha)c\right)^T W \left(y-\Phi(\alpha)c\right), ]
where (W) is a symmetric positive-definite weight matrix. If (W=L^TL), the problem is equivalent to ordinary least squares for the transformed quantities
[ \widetilde y=Ly, \qquad \widetilde\Phi(\alpha)=L\Phi(\alpha). ]
The reduced residual in whitened coordinates is consequently
[ \widetilde r(\alpha)
P_{\widetilde\Phi(\alpha)}^\perp \widetilde y. ]
You Watanabe established the invariant form of this weighted reduction in 1978 by expressing elimination as projection under the (W)-inner product rather than as a coordinate-dependent premultiplication. In that formulation, the fitted-value operator is
[ P_{\Phi,W}
\Phi \left(\Phi^T W\Phi\right)^\dagger \Phi^T W, ]
and the residual operator is (I-P_{\Phi,W}). The formulation also identifies whitening and variable elimination as equivalent transformations whenever the weight matrix is fixed and positive definite.
When the weights depend on the nonlinear parameters, whitening contributes additional derivative terms. The reduced objective still represents exact elimination of the linear coefficients, but its differential structure no longer follows solely from changes in the model basis.
Rank and regularity
Smooth variable projection theory generally assumes that the rank of (\Phi(\alpha)) remains constant throughout the region under consideration. Under this condition, the pseudoinverse and the corresponding projectors are differentiable whenever the basis matrix is differentiable. The reduced objective then has the regularity required by standard local analyses of nonlinear least squares.
At a rank-changing parameter value, the pseudoinverse can fail to be differentiable and may become unbounded in directions associated with a vanishing singular value. The projector itself can also change discontinuously if the dimension of the represented subspace changes. Such points reflect a loss of model identifiability rather than a defect created by the elimination procedure.
A singular value decomposition exposes this structure through
[ \Phi=U\Sigma V^T. ]
Small singular values indicate directions in the linear coefficient space that have little effect on the fitted observations. Truncating those directions changes the model subspace and therefore defines a regularized or rank-restricted variant of the original problem rather than an algebraically identical implementation.
Åke Björck analyzed these rank-sensitive least-squares structures through orthogonal factorizations and perturbation theory, clarifying the connection between numerical rank, pseudoinverse continuity, and the stability of projected residuals. This analysis applies directly to the linear subproblems embedded within variable projection.
Statistical interpretation
In a regression model with independent Gaussian errors of common variance, minimizing the residual sum of squares is equivalent to maximizing the likelihood. Eliminating (c) produces the profile likelihood for (\alpha), with the linear coefficients functioning as nuisance parameters. The reduced objective therefore preserves the likelihood maximizer of the original separable model.
The dimensional reduction does not remove uncertainty associated with the linear coefficients. After estimating (\alpha), covariance calculations for the full parameter vector include the dependence of (c(\alpha)) on the nonlinear parameters. Treating the reconstructed coefficients as independent of (\alpha) omits cross-covariance terms generated by the shared model geometry.
In generalized least squares, the same interpretation holds with an inverse error covariance matrix replacing the ordinary Euclidean metric. If the covariance model itself contains unknown nonlinear parameters, its determinant and parameter dependence also enter the likelihood. A weighted residual norm alone then represents only part of the statistical objective.
Computational characteristics
Variable projection reduces the explicit nonlinear search dimension from (p+n) to (p), while requiring a linear least-squares solution for each evaluation of the reduced residual. Its computational structure is therefore governed by repeated factorizations of (\Phi(\alpha)) and by derivatives of the parameterized basis.
The normal-equation representation
[ c= \left(\Phi^T\Phi\right)^{-1}\Phi^T y ]
is algebraically valid when (\Phi) has full column rank, but it squares the spectral condition number of the basis matrix. Orthogonal decompositions represent the same projection without forming (\Phi^T\Phi), and the singular value decomposition additionally represents rank deficiency. These are numerical properties of the embedded least-squares problem rather than changes to the variable-projection objective.
The reduction is most consequential when a model contains many linear coefficients but comparatively few nonlinear parameters. Typical mathematical forms include sums of exponentials with unknown decay rates, spectral models with nonlinear line locations, and rational approximations whose coefficients separate after fixing pole parameters. In each case, the projection removes amplitudes or other linear coefficients while retaining the nonlinear geometry that determines the model subspace.
Relation to other elimination methods
Variable projection is a special case of partial minimization, in which an objective (G(\alpha,c)) is replaced by
[ F(\alpha)=\inf_c G(\alpha,c). ]
Its distinguishing feature is that the inner minimization is a linear least-squares problem and therefore has a projector-based representation. This structure provides explicit residual derivatives and connects the reduced problem to orthogonal geometry.
The method is also related to the Schur complement. In a local quadratic model of the full least-squares problem, elimination of the linear parameter increment produces a reduced system for the nonlinear increment. Under constant-rank assumptions, that reduced system corresponds to projecting the nonlinear sensitivity matrix away from directions already represented by the linear basis.
Unlike block coordinate descent, variable projection does not define a sequence of alternating parameter estimates. The linear block remains at its conditional optimum throughout the reduced formulation, so every reduced residual is already orthogonal to the current linear model space.