A Generalized Variable Projection Algorithm for Least Squares Problems in Atmospheric Remote Sensing

TL;DR

The paper generalizes the variable projection (VP) algorithm to separable least squares problems with multiple right-hand sides that may vary in size, addressing limitations of the classical Golub–LeVeque formulation. It develops MRHS extensions, including naive, Golub–LeVeque, and Kaufman variants, and extends them to differing dataset lengths and model setups, all implemented in Python. Through synthetic and real OCO-2 radiance data for CO$_2$ retrieval, the study shows MRHS VP achieves comparable accuracy to standard nonlinear LS while dramatically speeding up computation as the number of RHS grows. The results suggest broad applicability to remote sensing and other scientific fitting problems with similar separable structures, especially when forward-model evaluations are expensive.

Abstract

This paper presents a solution for efficiently and accurately solving separable least squares problems with multiple datasets. These problems involve determining linear parameters that are specific to each dataset while ensuring that the nonlinear parameters remain consistent across all datasets. A well-established approach for solving such problems is the variable projection algorithm introduced by Golub and LeVeque, which effectively reduces a separable problem to its nonlinear component. However, this algorithm assumes that the datasets have equal sizes and identical auxiliary model parameters. This article is motivated by a real-world remote sensing application where these assumptions do not apply. Consequently, we propose a generalized algorithm that extends the original theory to overcome these limitations. The new algorithm has been implemented and tested using both synthetic and real satellite data for atmospheric carbon dioxide retrievals. It has also been compared to conventional state-of-the-art solvers, and its advantages are thoroughly discussed. The experimental results demonstrate that the proposed algorithm significantly outperforms all other methods in terms of computation time, while maintaining comparable accuracy and stability. Hence, this novel method can have a positive impact on future applications in remote sensing and could be valuable for other scientific fitting problems with similar properties.

Figures (7)

• Figure 1: Four exemplary soundings of frame 1728.0 from the OCO-2 level 1b (L1B) measurement product Crisp, each displaying radiance spectra in units [$erg /(\s\; \cm\squared\; sr\;\cm\tothe{-1})$] from both the strong and weak bands with 809.0 and 651.0 spectral pixels each.
• Figure 2: Total column average dry air mole fraction (i.e., “concentration”) of carbon dioxide, denoted as $x\ce{CO2}$ in the unit (parts per million), from the level 2 (L2) retrieval product of OCO-2 L2dataCrisp2 for one exemplary frame (cf. Figure \ref{['fig:spectra']}), including eight soundings.
• Figure 3: Distribution plots of the relative errors (difference between exact and fitted $\alpha_{\ce{CO2}}$ results) on the right and the corresponding standard deviations from the exact solution on the left for both fitting setups single and MRHS with increasing signal-to-noise ratios (SNRs).
• Figure 4: Mean sigma of regression for both single and MRHS fits for different noisy spectra. The dashed lines correspond to fitted hyperbolas.
• Figure 5: Vertical box plots (with horizontal offset for better distinction) of the confidence bounds of the $\alpha_{\ce{CO2}}$ parameter achieved by each method for different numbers of datasets.