Research Radar — 2026-11-03

Methodology Review: Variable projection: solve the easy parameters before searching the hard ones

English · About 10 minutes

Main reading inspected in author full text; related papers primarily through publisher abstracts. No new biological experiment or benchmark reproduction. First-online, issue and digitization dates distinguished.

The intuition

Imagine fitting a curve whose shape depends on a decay rate, while its vertical offset and height are also unknown. A general optimizer can search all three quantities together. But once you choose the decay rate, finding the best offset and height is an ordinary linear least-squares problem. Variable projection exploits that asymmetry: let a reliable inner calculation solve the easy coefficients exactly, and reserve nonlinear search for the parameters that genuinely change the shape. This is a strategy for organizing a fitting problem, not a new biological model.

The main reading by O’Leary and Rust explains how to implement this separation, including numerical treatment of the linear subproblem. The reduced objective asks a precise question: for this proposed nonlinear shape, how small can the residual become after choosing the best linear coefficients? Search that objective rather than asking the outer optimizer to rediscover the coefficients at every step. This differs from alternating one tentative coefficient update with one shape update: the inner problem is solved to its required accuracy for each outer evaluation. It does not remove the need to choose a sensible model. [R1]

A small worked example

Consider three invented measurements at times zero, one and two: their values are 3, 2 and 1.5. Fit \(y_i\approx b+a q^{t_i}\). Here \(i\) labels an observation, \(t_i\) is its time, \(y_i\) its value, \(b\) a baseline, \(a\) an amplitude and \(q\) the fraction of the decaying component retained per time step. Restrict \(q\) to be positive and no greater than one. This is a teaching example, not experimental evidence. The important distinction is that the model is linear in \(a\) and \(b\) when \(q\) is held fixed.

Try \(q=1/2\) first. Choosing \(b=1\) and \(a=2\) predicts 3, 2 and 1.5 exactly, so the sum of squared residuals is zero. Now try \(q=1/4\). Solving the two-coefficient least-squares problem gives \(b=3/2\) and \(a=32/21\). The predictions are \(127/42\), \(79/42\) and \(67/42\). Subtracting these from the observations gives residuals \(-1/42\), \(5/42\) and \(-2/21\). Their squared sum is \(1/42\). Thus the outer search can compare two numbers, zero and one forty-second, without separately searching amplitude and baseline for either trial.

One more trial exposes a failure mode. At \(q=1\) the decaying component is constant, so the model depends only on the sum of amplitude and baseline. The best constant is the observation mean, \(13/6\), and the squared residual sum is \(7/6\). Infinitely many coefficient pairs have that same sum. A fitting routine may return one pair, but that does not make the two coefficients individually identifiable. As the decay factor approaches one, distinguishing a slowly decaying signal from a constant offset also becomes fragile. A good curve fit and trustworthy parameter estimates are different accomplishments.

To reproduce the arithmetic, form one column of ones and one column containing the chosen decay factor raised to each time. Solve for the two coefficients, calculate three predictions, then square and add their residuals. Repeat after slightly perturbing the observations. Ask whether the curve remains similar while the coefficient estimates change. This exercise separates three questions often compressed into a single optimization score: whether the calculation is correct, whether the chosen family can describe the observations, and whether the available sampling times can distinguish its parameters.

The necessary mathematics

In the general problem, \(y\in\mathbb{R}^{m}\) is the observation vector, \(m\) is the number of measurements, \(\theta\in\mathbb{R}^{d}\) contains \(d\) nonlinear parameters, and \(A(\theta)\in\mathbb{R}^{m\times p}\) is a design matrix with \(p\) columns. The coefficient vector \(c\in\mathbb{R}^{p}\) is linear. The Euclidean norm, denoted \(\lVert\cdot\rVert_2\), measures residual size. The joint objective is

For each fixed shape, choose its least-squares coefficients before comparing it with other shapes. This is the separation formalized in the foundational paper. [R2]

Write \(A(\theta)^+\) for the Moore–Penrose pseudoinverse. A minimum-norm inner solution and the reduced objective are

The star labels the chosen inner solution; \(F\) is the best attainable squared residual at that shape. Minimizing \(F\) over the nonlinear parameters is equivalent to minimizing the joint objective. Geometrically, the fitted vector is the orthogonal projection of the observations onto the column space of the current design matrix. At deficient rank the coefficient vector need not be unique, even though the fitted vector is. Smooth derivative arguments need rank qualifications; the original paper explicitly treats constant rank. [R2]

Assumptions and failure modes

These formulas describe unconstrained, unweighted squared-error fitting. If amplitudes must be nonnegative, or the loss is a count likelihood rather than squared error, the inner problem changes: do not retain the pseudoinverse formula simply because the model contains linear-looking coefficients. The main reading supplies computational guidance for the least-squares setting, including SVD-based treatment of the design matrix. A numerically stable solver cannot manufacture information absent from nearly indistinguishable columns. Record parameter constraints, scaling and rank diagnostics rather than treating a successful solver exit as scientific validation. [R1]

Reducing the number of nonlinear search variables does not guarantee that the remaining objective is convex, globally solved or cheap. Osborne analyzes separable least squares and Gauss–Newton behavior, including small-residual behavior and large-data asymptotic arguments under additive-error assumptions. Those conditions matter; they are not a blanket convergence promise for arbitrary noisy data. An outer method still needs initialization and termination choices. Repeated inner solves can also dominate runtime when the design matrix is large. The appropriate comparison is accuracy and total computational cost for the actual problem, not parameter count alone. [R3]

Verified applications

Mullen and van Stokkum report variable-projection developments for matrix-valued measurements and constrained coefficients, with case studies spanning time-resolved spectroscopy, microscopy and chromatography–mass spectrometry. These are established measurement applications rather than a biological connection invented for this column. The inspected publisher abstract supports that scope, but it does not independently establish performance on a particular current omics assay. The paper appeared online in September 2008 and in a 2009 issue; those are versions of one reading, not two discoveries. Read it for the modeling structure and extensions, not as a universal endorsement of exponential fitting. [R4]

Possible transfers — speculative

A possible omics transfer is a time-course model with shared nonlinear temporal shapes and feature-specific linear amplitudes. This is explicitly speculative here. Start with a small, independently justified model and ask whether sampling times distinguish its shapes before fitting thousands of features. Compare held-out predictions with a simpler baseline, and inspect sensitivity to initialization and the observation transformation. A squared-error working model on transformed measurements needs justification; this review does not claim that it reproduces the count-generating process, identifies causal regulation, or validates any autonomous experimental-design system.

What to try and what to read next

Use variable projection when the scientific model already contains a genuinely linear block inside a smaller nonlinear problem. Begin with the three-point exercise, including its rank-deficient endpoint, before considering a large dataset. Read the main paper for implementation, the 1973 paper for the foundation, Osborne for convergence qualifications, and the measurement paper for established applications. A practical stopping rule for a first prototype is not merely a small residual: you should be able to explain which parameters were eliminated, how their inner problem was solved, and which parameter combinations the data cannot distinguish.

Reading list

• [R1] Primary reading: Variable projection for nonlinear least squares problems. Dianne P. O’Leary and Bert W. Rust (2012). Main reading: start with the model and implementation discussion in sections 1–2, not the whole software specification. Author manuscript is openly accessible; the journal issue year is 2013, distinct from first publication.

• [R2] Related reading: The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate. Gene H. Golub and Victor Pereyra (1973). Foundational derivation: distinguish elimination of coefficients from alternating updates, and notice the constant-rank qualification.

• [R3] Related reading: Separable least squares, variable projection, and the Gauss-Newton algorithm. M. R. Osborne (2007). Convergence perspective: reduction in dimension is not a promise of globally optimal fitting or universally faster convergence.

• [R4] Related reading: The variable projection algorithm in time-resolved spectroscopy, microscopy and mass spectrometry applications. Katharine M. Mullen and Ivo H. M. van Stokkum (2008). Established measurement applications and extensions; inspected evidence is the abstract, not an independent reproduction of its case studies.