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

Adelina Bärligea
Remote Sensing Technology Institute
German Aerospace Center
82234 Oberpfaffenhofen
& Philipp Hochstaffl
Remote Sensing Technology Institute
German Aerospace Center
82234 Oberpfaffenhofen
& Franz Schreier
Remote Sensing Technology Institute
German Aerospace Center
82234 Oberpfaffenhofen

###### 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.

_Keywords_ separable least squares $\cdot$ nonlinear optimization $\cdot$
python $\cdot$ inverse problems $\cdot$ trace gas retrieval $\cdot$
atmospheric composition $\cdot$ carbon dioxide $\cdot$ infrared spectroscopy

##  1 Introduction

One of the fundamental tasks in scientific computing is to find the unknown
vector of parameters $\boldsymbol{x}\in\mathbb{R}^{n+p}$, which, for given
data $(y_{i},t_{i})$, $i=1,\ldots m$, $m\geq n+p$, and model function
$\eta(\mathbf{x},t)$, solves the least squares problem

$\min_{\boldsymbol{x}}\sum_{i=1}^{m}(y_{i}-\eta(\boldsymbol{x},t_{i}))^{2}=\min_{\boldsymbol{x}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\eta}(\boldsymbol{x})\bigr{\|}_{2}^{2},$
(1)

where the vectors $\boldsymbol{\eta}\in\mathbb{R}^{m}$ and
$\boldsymbol{y}\in\mathbb{R}^{m}$ have the entries
$[\boldsymbol{\eta}(\boldsymbol{x})]_{i}=\eta(\boldsymbol{x},t_{i})$ and
$[\boldsymbol{y}]_{i}=y_{i}$, respectively. If the model function is
nonlinear, the minimization problem is solved iteratively by using step-length
or trust-region methods, such as the Levenberg–Marquardt algorithm (Moré,
1978).

In separable least squares problems, the model function $\eta$ is a linear
combination of nonlinear functions $\varphi_{j}$, $j=1,\ldots,n$, i.e.,

$\eta(\boldsymbol{\alpha},\boldsymbol{\beta},t)=\sum_{j=1}^{n}\beta_{j}\varphi_{j}(\boldsymbol{\alpha},t),$
(2)

and the vector of parameters $\boldsymbol{x}$ is split into a vector of linear
parameters $\boldsymbol{\beta}\in\mathbb{R}^{n}$, and a vector of nonlinear
parameters, $\boldsymbol{\alpha}\in\mathbb{R}^{p}$.

In 1973, Golub and Pereyra (1973) proposed a powerful method for solving such
problems, called variable projection (VP). Specifically, the problem is
reduced to a nonlinear least squares problem involving $\boldsymbol{\alpha}$
only, and a linear least squares problem involving $\boldsymbol{\beta}$ only.
Thus, the dimension of the nonlinear problem to be solved is reduced from
$n+p$ to $p$. Later on, Golub and LeVeque (1979) extended the method to the
case of multiple datasets $\boldsymbol{y}_{k}$ with $k=1,\ldots,s$. In this
type of problem, a vector of linear parameters $\boldsymbol{\beta}_{k}$ is
associated with a dataset $\mathbf{y}_{k}$, while the vector of nonlinear
parameters $\boldsymbol{\alpha}$ corresponds to all datasets.

One example of this type of problem can be found in the retrieval of the
atmospheric composition from passive remote sensing measurements (cf.
(Hochstaffl et al., 2022)). Here, an atmospheric radiative transfer model with
molecular composition parameters (nonlinear) and reflectivity parameters
(linear) is fitted to spectral radiance measurements to determine the
underlying atmospheric state parameters. For long-lived molecules, such as CO2
or CH4, which are homogeneously spread throughout the atmosphere, it is
possible to fit several (nearby) observations for one concentration value
while the surface reflectivity differs for each measurement, making this a
separable problem with multiple right-hand sides.

Fitting, for example, $$4$\times$4$=$16$$ radiance spectra with three distinct
linear reflectivity variables simultaneously for two nonlinear molecular
concentration variables, the total number of unknowns for conventional solvers
adds up to $50$ ($=$16$\cdot$3$+$2$$). However, with a separable approach,
such as variable projection, one can reduce this size by a factor of $25$.

However, the variable projection algorithm established by Golub and LeVeque
(1979) to solve such problems with multiple right-hand sides is based on the
assumption that all datasets have the same lengths and corresponding nonlinear
models.

In this article, the theory by Golub and LeVeque (1979) is, therefore,
elaborated and further enhanced for scenarios such as differently sized
datasets or varying nonlinear model setups (see Section 2). The modifications
were not only necessary in order to apply their method to the example
described above, but this new algorithm could also be useful for other
scientific fitting problems. The method was implemented in Python (see Section
3) and applied to the least squares problem arising in atmospheric trace gas
retrieval (see Section 4). In Section 5, experimental results are discussed
and concluded.

##  2 Theoretical Background

Assuming a model that is a linear combination of nonlinear functions, such as
(2), one can define the arising separable least squares problem as follows:

$\min_{\boldsymbol{\alpha},\boldsymbol{\beta}}\sum_{i=1}^{m}\Bigl{[}y_{i}-\sum_{j=1}^{n}\beta_{j}\varphi_{j}(\boldsymbol{\alpha},t_{i})\Bigr{]}^{2}=\min_{\boldsymbol{\alpha},\boldsymbol{\beta}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\beta}\bigr{\|}_{2}^{2},$
(3)

where the matrix $\boldsymbol{\Phi}(\boldsymbol{\alpha})\in\mathbb{R}^{m\times
n}$ has the entries
$[\boldsymbol{\Phi}(\boldsymbol{\alpha})]_{ij}=\varphi_{j}(\boldsymbol{\alpha},t_{i})$.

###  2.1 The Variable Projection Method

The minimization problem in (2), written as

$\min_{\boldsymbol{\alpha}}\Bigl{(}\min_{\boldsymbol{\beta}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\beta}\bigr{\|}_{2}^{2}\Bigr{)},$
(4)

is reduced to two least squares problems. Golub and Pereyra (1973) pointed out
that for any fixed $\boldsymbol{\alpha}$, one considers the merely linear
least squares problem

$\min_{\boldsymbol{\beta}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\beta}\bigr{\|}_{2}^{2}.$
(5)

which is solved by

$\boldsymbol{\beta}(\boldsymbol{\alpha})=\boldsymbol{\Phi}^{\dagger}(\boldsymbol{\alpha})\boldsymbol{y},$
(6)

where
$\boldsymbol{\Phi}^{\dagger}=(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi})^{-1}\boldsymbol{\Phi}^{\top}$
is the generalized inverse of the matrix $\boldsymbol{\Phi}$. Second, for the
functional

$\mathcal{F}(\boldsymbol{\alpha})=\min_{\boldsymbol{\beta}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\beta}\bigr{\|}_{2}^{2}=\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\Phi}^{\dagger}(\boldsymbol{\alpha})\boldsymbol{y}\bigr{\|}_{2}^{2},$
(7)

which is obtained by substituting the expression of
$\boldsymbol{\beta}(\boldsymbol{\alpha})$ given by Equation (6) into the
residual function, one considers the nonlinear least squares problem

$\min_{\boldsymbol{\alpha}}\mathcal{F}(\boldsymbol{\alpha})=\min_{\boldsymbol{\alpha}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\Phi}^{\dagger}(\boldsymbol{\alpha})\boldsymbol{y}\bigr{\|}_{2}^{2}=\min_{\boldsymbol{\alpha}}\bigl{\|}\boldsymbol{P}_{\boldsymbol{\Phi}(\boldsymbol{\alpha})}^{\perp}\boldsymbol{y}\bigr{\|}_{2}^{2},$
(8)

where
$\boldsymbol{P}_{\boldsymbol{\Phi}(\boldsymbol{\alpha})}^{\perp}=\boldsymbol{I}_{m\times
m}-\boldsymbol{P}_{\boldsymbol{\Phi}(\boldsymbol{\alpha})}\in\mathbb{R}^{m\times
m},$ and
$\boldsymbol{P}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})=\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\Phi}^{\dagger}(\boldsymbol{\alpha})\in\mathbb{R}^{m\times
m}$ is the orthogonal projection operator onto the column space of the matrix
$\boldsymbol{\Phi}(\boldsymbol{\alpha})$. Thanks to this formulation, the
method is called variable projection (VP).

Thus, the method of solution involves, in the first step, the computation of
the minimizer
$\widehat{\boldsymbol{\alpha}}=\arg\min_{\boldsymbol{\alpha}}\mathcal{F}(\boldsymbol{\alpha})$,
and in the second step, the computation of the vector of linear parameters as
$\widehat{\boldsymbol{\beta}}=\boldsymbol{\Phi}^{\dagger}(\widehat{\boldsymbol{\alpha}})\boldsymbol{y}$.
Golub and Pereyra (1973) showed that this solution method yields the correct
result, under the assumption that $\boldsymbol{\Phi}(\boldsymbol{\alpha})$ has
a locally constant rank in the neighborhood of
$\widehat{\boldsymbol{\alpha}}$. In order to minimize a nonlinear problem,
such as (7), the derivative of the objective function, with respect to the
unknown variables, is needed. If the rank of matrix
$\boldsymbol{\Phi}(\boldsymbol{\alpha})$ was not constant across the points
throughout the iteration where its derivative is calculated, the pseudo-
inverse $\boldsymbol{\Phi^{\dagger}}(\boldsymbol{\alpha})$ would not be a
continuous function and, therefore, not differentiable.

This separation technique has powerful advantages over conventional
algorithms, which all follow from the fact that the nonlinear functional
$\mathcal{F}(\boldsymbol{\alpha})$ only depends on $\boldsymbol{\alpha}$ (see
also the review by Golub and Pereyra (2003)). Since $\boldsymbol{\alpha}$ is a
vector of length $p$ and $\boldsymbol{\beta}$ of length $n$, the VP method
effectively turns the original nonlinear problem of $n+p$ variables into one
of only $p$. This reduction of the parameter space of the problem results in a
smaller size of the problem’s Jacobian and, therefore, requires less time for
its computation. Moreover, O’Leary and Rust (2013) pointed out that the
reduced parameter space may also lead to a reduced number of local minimizers,
making it more likely to find the global minimum instead of a local one. To
summarize, the VP solver is a lot more efficient and also converges better
than conventional methods with no separation (Ruhe and Wedin, 1980; Golub and
Pereyra, 2003). Moreover, the need for a smaller initial guess vector can
generally lead to a better-conditioned and more stable problem.

For solving the nonlinear least squares problem (8), one needs to calculate
the partial derivatives
$\partial\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}/\partial\alpha_{l}$. In
this regard, Golub and Pereyra (1973) provide the computational formula

$\frac{\partial\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}}{\partial\alpha_{l}}=-\Bigl{[}\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}\frac{\partial\boldsymbol{\Phi}}{\partial\alpha_{l}}\boldsymbol{\Phi}^{\dagger}+\Bigl{(}\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}\frac{\partial\boldsymbol{\Phi}}{\partial\alpha_{l}}\boldsymbol{\Phi}^{\dagger}\Bigr{)}^{\top}\Bigr{]},$
(9)

which is proved in the Appendix.

Kaufman (1975) proposed a modification of the original method by making use of
a QR decomposition of the matrix $\boldsymbol{\Phi}$:

$\boldsymbol{\Phi}=\boldsymbol{Q}\boldsymbol{R}=(\begin{array}[]{cc}\boldsymbol{Q}_{1}&\boldsymbol{Q}_{2})\left(\begin{array}[]{c}\boldsymbol{R}_{1}\\\
\boldsymbol{0}_{(m-n)\times
n}\end{array}\right)=\boldsymbol{Q}_{1}\boldsymbol{R}_{1},\end{array}$ (10)

where $\boldsymbol{Q}\in\mathbb{R}^{m\times m}$,
$\boldsymbol{Q}_{1}\in\mathbb{R}^{m\times n}$, and
$\boldsymbol{Q}_{2}\in\mathbb{R}^{m\times(m-n)}$ are matrices with orthonormal
columns (i.e., $\boldsymbol{Q}^{\top}\boldsymbol{Q}=\boldsymbol{I}_{m\times
m}$, $\boldsymbol{Q}_{1}^{\top}\boldsymbol{Q}_{1}=\boldsymbol{I}_{n\times n}$,
$\boldsymbol{Q}_{2}^{\top}\boldsymbol{Q}_{2}=\boldsymbol{I}_{(m-n)\times(m-n)}$)
and $\boldsymbol{R}_{1}\in\mathbb{R}^{n\times n}$ is an upper triangular, non-
singular matrix. In this case, the generalized inverse
$\boldsymbol{\Phi}^{\dagger}$, satisfying the relation
$\boldsymbol{\Phi}^{\dagger}\boldsymbol{\Phi}=\boldsymbol{I}_{n\times n}$,
computes as

$\boldsymbol{\Phi}^{\dagger}=(\begin{array}[]{cc}\boldsymbol{R}_{1}^{-1}&\boldsymbol{0}_{n\times(m-n)})\left(\begin{array}[]{c}\boldsymbol{Q}_{1}^{\top}\\\
\boldsymbol{Q}_{2}^{T}\end{array}\right)=\boldsymbol{R}_{1}^{-1}\boldsymbol{Q}_{1}^{\top}.\end{array}$
(11)

From the invariance of the $2$-norm under orthogonal transformations and
$\boldsymbol{Q}^{\top}\boldsymbol{\Phi}=\boldsymbol{Q}^{\top}\boldsymbol{Q}\boldsymbol{R}=\boldsymbol{R}$),
it follows that

$\displaystyle\mathcal{F}(\boldsymbol{\alpha})$
$\displaystyle=\min_{\boldsymbol{\beta}}\bigl{\|}\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\beta}\bigr{\|}_{2}^{2}$
$\displaystyle=\min_{\boldsymbol{\beta}}\bigl{\|}\boldsymbol{Q}^{\top}(\boldsymbol{\alpha})(\boldsymbol{y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\beta})\bigr{\|}_{2}^{2}$
$\displaystyle=\min_{\boldsymbol{\beta}}\Bigl{\|}\left(\begin{array}[]{c}\boldsymbol{Q}_{1}^{\top}(\boldsymbol{\alpha})\boldsymbol{y}\\\
\boldsymbol{Q}_{2}^{\top}(\boldsymbol{\alpha})\boldsymbol{y}\end{array}\right)-\left(\begin{array}[]{c}\boldsymbol{R}_{1}(\boldsymbol{\alpha})\\\
\boldsymbol{0}_{(m-n)\times
n}\end{array}\right)\boldsymbol{\beta}\Bigr{\|}_{2}^{2}$ (16)
$\displaystyle=\min_{\boldsymbol{\beta}}\bigl{\|}\boldsymbol{Q}_{1}^{\top}(\boldsymbol{\alpha})\boldsymbol{y}-\boldsymbol{R}_{1}(\boldsymbol{\alpha})\boldsymbol{\beta}\bigr{\|}_{2}^{2}+\bigl{\|}\boldsymbol{Q}_{2}^{\top}(\boldsymbol{\alpha})\boldsymbol{y}\bigr{\|}_{2}^{2}.$
(17)

Since the optimal $\boldsymbol{\beta}$ for any given $\boldsymbol{\alpha}$ is

$\boldsymbol{\beta}(\boldsymbol{\alpha})=\boldsymbol{\Phi}^{\dagger}(\boldsymbol{\alpha})\boldsymbol{y}=\boldsymbol{R}_{1}^{-1}(\boldsymbol{\alpha})\boldsymbol{Q}_{1}^{\top}(\boldsymbol{\alpha})\boldsymbol{y},$
(18)

we find
$\mathcal{F}(\boldsymbol{\alpha})=\bigl{\|}\boldsymbol{Q}_{2}^{T}(\boldsymbol{\alpha})\boldsymbol{y}\bigr{\|}_{2}^{2}$,
showing that the nonlinear least squares problem reduces to

$\min_{\boldsymbol{\alpha}}\mathcal{F}(\boldsymbol{\alpha})=\min_{\boldsymbol{\alpha}}\bigl{\|}\boldsymbol{Q}_{2}^{\top}(\boldsymbol{\alpha})\boldsymbol{y}\bigr{\|}_{2}^{2}.$
(19)

Moreover, for derivative calculations, Kaufman (1975) proposed the simplified
formula

$\frac{\partial\boldsymbol{Q}_{2}^{T}}{\partial\alpha_{l}}=-\boldsymbol{Q}_{2}^{\top}\frac{\partial\boldsymbol{\Phi}}{\partial\alpha_{l}}\boldsymbol{\Phi}^{\dagger},$
(20)

which is justified in the Appendix. This was shown to save function and
gradient evaluation costs and, therefore, reduce the computing time per
iteration, which is why this simplified version of (9) became well
established.

In conclusion, the minimization problems (8) and (19) are equivalent, but the
size of the matrix $\boldsymbol{Q}_{2}^{\top}$ is smaller than that of the
matrix $\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}$. Consequently, an
algorithm for solving Equation (19) should be more efficient.

###  2.2 Multiple Right-Hand Sides (MRHS)

As Golub and LeVeque (1979) pointed out, there can be optimization problems
where multiple sets of data
$\boldsymbol{y_{1}},\boldsymbol{y_{2}},\ldots,\boldsymbol{y_{s}}$ are to be
fit to a model function, such as (2). If the model parameters are to vary for
each set $\boldsymbol{y_{k}}$, this will just result in $s$ distinct separable
problems, such as (3). There are, however, cases where only the linear
variables are specific to each dataset, while the nonlinear variables have to
hold for all available data simultaneously. By exploiting separability, the
sizes of such minimization problems can be reduced from $sn+p$ to just $p$
unknown variables, which is even greater than with only one dataset.

Separable problems with $s$ right-hand side(s) (RHS) can be posed as the
minimization of

$\min_{\boldsymbol{\alpha},\boldsymbol{B}}\
\left\|\boldsymbol{Y}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{B}\right\|_{\text{F}}^{2},$
(21)

where the nonlinear parameter vector is $\boldsymbol{\alpha}\in\mathbb{R}^{p}$
and the matrix is $\boldsymbol{B}\in\mathbb{R}^{n\times s}$, containing the
linear parameter vectors $\boldsymbol{\beta_{k}}$ for each RHS
$\boldsymbol{y_{k}}$ with $k=1,\ldots,s$ as its columns, using the Frobenius
norm $\left\|\cdot\right\|_{\text{F}}$. Here, the data matrix
$\boldsymbol{Y}=(\boldsymbol{y_{1}}\ldots\boldsymbol{y_{s}})\in\mathbb{R}^{m\times
s}$ is fit to a single model of the form
$\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{B}$, where
$\boldsymbol{\Phi}(\boldsymbol{\alpha})\in\mathbb{R}^{m\times n}$ is as before
(cf. Equation (3)).

####  2.2.1 Naive Approach

The first intuitive approach, which was also mentioned by Golub and LeVeque
(1979), is to reformulate (21) using the matrix

$\boldsymbol{G}(\boldsymbol{\alpha})=\begin{pmatrix}\boldsymbol{\Phi}(\boldsymbol{\alpha})&&\ldots&\boldsymbol{0}\\\
&\boldsymbol{\Phi}(\boldsymbol{\alpha})&&\vdots\\\ \vdots&&\ddots&\\\
\boldsymbol{0}&\ldots&&\boldsymbol{\Phi}(\boldsymbol{\alpha})\end{pmatrix}\in\mathbb{R}^{ms\times
ns},$ (22)

and the vectors $\boldsymbol{\tilde{y}}=\begin{pmatrix}\boldsymbol{y_{1}}\\\
\vdots\\\ \boldsymbol{y_{s}}\end{pmatrix}\in\mathbb{R}^{ms}$ and
$\boldsymbol{\tilde{\beta}}=\begin{pmatrix}\boldsymbol{\beta_{1}}\\\ \vdots\\\
\boldsymbol{\beta_{s}}\end{pmatrix}\in\mathbb{R}^{ns}$, such that a problem of
the original vectorial form

$\min_{\boldsymbol{\alpha},\boldsymbol{\tilde{\beta}}}\
\left\|\boldsymbol{\tilde{y}}-\boldsymbol{G}(\boldsymbol{\alpha})\boldsymbol{\tilde{\beta}}\right\|^{2},$
(23)

arises. This formulation allows for solving the separable problem with
multiple RHS, by means of the variable projection algorithm already discussed.
However, even for a moderate number of datasets, the matrix
$\boldsymbol{G}(\boldsymbol{\alpha})$ becomes overly large. Moreover, the
sparse structure and the fact that all diagonal blocks in (22) are the same
can be better utilized in the earlier formulation (21).

####  2.2.2 Golub–LeVeque Approach

Therefore, Golub and LeVeque (1979) suggested a different approach: Starting
from (21), one can, in the same manner as Equation (7), exploit the problem’s
separability and reduce it to a purely nonlinear minimization problem of the
form

$\min_{\boldsymbol{\alpha}}\
\left\|\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})\boldsymbol{Y}\right\|_{\text{F}}^{2},$
(24)

with the orthogonal projector
$\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})=\boldsymbol{I}-\boldsymbol{\Phi}(\boldsymbol{\alpha})\boldsymbol{\Phi^{\dagger}}(\boldsymbol{\alpha})$
as before. Acknowledging that the Frobenius norm in (24) is equivalent to the
$2$-norm of a vector function $\boldsymbol{z}(\boldsymbol{\alpha})$, set up as

$\boldsymbol{z}(\boldsymbol{\alpha})=\begin{pmatrix}\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})\>\boldsymbol{y_{1}}\\\
\vdots\\\
\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})\>\boldsymbol{y_{s}}\end{pmatrix}\in\mathbb{R}^{ms},$
(25)

this can be minimized with any of the established methods for nonlinear least
squares problems, such as the Levenberg–Marquardt algorithm, which is used in
(Golub and LeVeque, 1979). The Jacobian matrix of
$\boldsymbol{z}(\boldsymbol{\alpha})$ can be calculated analogously by
defining its $l{\text{th}}$ column as

$\frac{\partial\boldsymbol{z}(\boldsymbol{\alpha})}{\partial\alpha_{l}}=\begin{pmatrix}\frac{\partial\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})}{\partial\alpha_{l}}\>\boldsymbol{y_{1}}\\\
\vdots\\\
\frac{\partial\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})}{\partial\alpha_{l}}\>\boldsymbol{y_{s}}\end{pmatrix}\in\mathbb{R}^{ms},$
(26)

where
$\frac{\partial\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})}{\partial\alpha_{l}}$
is exactly as in (9).

####  2.2.3 Kaufman Approach

Since the Frobenius norm is invariant under the orthogonal transformation, the
above method can also be written in terms of the QR decomposition of
$\boldsymbol{\Phi}(\boldsymbol{\alpha})$ outlined in Section 2.1. This
approach based on Kaufman (1975) can, therefore, be established by replacing
the orthogonal projector
$\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}(\boldsymbol{\alpha})\in\mathbb{R}^{m\times
m}$ in all of the above equations of the Golub–LeVeque method, with the
smaller orthogonal matrix $\boldsymbol{Q_{2}}^{\top}\in\mathbb{R}^{(m-n)\times
m}$ derived in (10). Both versions of the approach are included in the Python
implementation outlined in Section 3 and are, therefore, the subjects of the
numerical experiments performed in Section 4.

###  2.3 Extensions to the Golub–LeVeque Approach

Working with real measurements, which can be subject to errors or missing data
points, it is not always possible to use datasets that all have the exact same
number of data points. In this case, each dataset $\boldsymbol{y_{k}}$ has a
specific size $m_{k}$, such that the matrix depiction used in (21) does not
hold any longer.

Moreover, the nonlinear model functions stored in the matrix
$\boldsymbol{\Phi}(\boldsymbol{\alpha})$ often depend on further auxiliary
parameters, which may vary for each individual measurement (e.g., observation
angle). Thus, there can be different
$\boldsymbol{\Phi_{k}}(\boldsymbol{\alpha})$ of size $m_{k}\times n$.

None of these aspects were explicitly mentioned by Golub and LeVeque (1979) or
Kaufman and Sylvester (1992), who later simplified the Golub–LeVeque method
further, with the assumption of equal lengths of the $\boldsymbol{y_{k}}$. It
was only mentioned by Kaufman (2010) that the TIMP package by Mullen and van
Stokkum (2007) allows for differently sized data vectors, but not for
differently constructed $\boldsymbol{\Phi}$s.

For the implementation introduced in this paper, both of these changes are
taken into account. It can be seen that the special structure of vector (25)
and matrix (26) can be exploited to rewrite

$\boldsymbol{z}(\boldsymbol{\alpha})=\begin{pmatrix}\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi_{1}}}(\boldsymbol{\alpha})\>\boldsymbol{y_{1}}\\\
\vdots\\\
\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi_{s}}}(\boldsymbol{\alpha})\>\boldsymbol{y_{s}}\end{pmatrix}\in\mathbb{R}^{m_{1}+\ldots+m_{s}}$
(27)

as a stack of (differently sized) vectors and, likewise, the $l{\text{th}}$
column of its Jacobian

$\frac{\partial\boldsymbol{z}(\boldsymbol{\alpha})}{\partial\alpha_{l}}=\begin{pmatrix}\frac{\partial\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi_{1}}}(\boldsymbol{\alpha})}{\partial\alpha_{l}}\>\boldsymbol{y_{1}}\\\
\vdots\\\
\frac{\partial\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi_{s}}}(\boldsymbol{\alpha})}{\partial\alpha_{l}}\>\boldsymbol{y_{s}}\end{pmatrix}\in\mathbb{R}^{(m_{1}+\ldots+m_{s})s}\
\mathrm{.}$ (28)

This necessary modification naturally increases the computational expense
compared to the original Golub–LeVeque method, as matrix $\boldsymbol{\Phi}$
and its derivative have to be calculated $s$ times instead of once. This also
reduces the gain in efficiency that one would have had from exchanging
$\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}$ with $\boldsymbol{Q_{2}}^{\top}$,
as the QR decomposition now needs to be calculated for every single matrix
$\boldsymbol{\Phi}_{k}$. However, this modified VP method for multiple right-
hand sides is still significantly more efficient than using a standard
nonlinear optimizer for the same problem, as will be shown in Section 4.

##  3 Implementation

The first FORTRAN implementation of a variable projection algorithm that
allows for multiple right-hand sides was VARP2 (LeVeque, 1985), which was
developed by Randy LeVeque. It is a modification of the subroutine VARPRO
(Bolstad, 1977) for classical separable problems with only a single RHS. Both
work with user-provided Jacobians. Another efficient implementation of a
solver for multiple datasets is the TIMP package by Mullen and van Stokkum
(2007), written in the statistical computing language R, where the Jacobian is
calculated by finite difference approximations.

This work was partly motivated by the fact that many scientific computing
problems are nowadays implemented in Python. However, the only modern
implementation of a VP algorithm is a MATLAB code by O’Leary and Rust (2013),
which does not allow for multiple right-hand sides, as discussed in Section
2.3. Here, an implementation is introduced where the standard VP algorithm by
O’Leary and Rust (2013) is enhanced for multiple datasets, as described in
Section 2.

The MATLAB code by O’Leary and Rust (2013) was specifically designed to be
brief and easily understandable, so it would be well-suited for translation.
O’Leary and Rust (2013) argued that many of the established implementations,
such as VARPRO (Bolstad, 1977), and PORT library (Fox et al., 1978)
subroutines, such as NSF/NSG (Kaufman, 1975), are somewhat outdated today and
lack readability; thus, they proposed an efficient present-day implementation
written in an interpreted language, with the advantage that it can easily be
enhanced or modified. One reason for their code’s brevity is that they used
built-in MATLAB functions, such as lsqnonlin.m, to solve the nonlinear least
squares problem (7), and svd.m to solve the remaining linear problems via
singular value decomposition, instead of writing their own, making the
algorithm modular. Another advantage of their code is the variety of
statistical diagnostics it offers, some of which were also modified for
multiple datasets and used for the assessment in Section 4.

###  3.1 Algorithm for a Single Right-Hand Side

A short outline of a standard VP algorithm (such as the one by O’Leary and
Rust (2013)) can be formulated as this:

  1. 1.

The user supplies the measurement vector $\boldsymbol{y}$, the number of
linear parameters $n$, additional independent variables for the model, and an
initial guess for the nonlinear variables $\boldsymbol{\alpha_{0}}$. In
addition, a subroutine/function which will be called ADA (note that a detailed
description of how to set up and implement this subroutine for a variable
projection algorithm, such as the one presented here, can be found in (O’Leary
and Rust, 2013), Sections 2.3 and 2.4) needs to be provided, which, for the
given input $n$ and $\boldsymbol{\alpha}$, calculates the nonlinear model
matrix $\boldsymbol{\Phi}$ and its partial derivatives
$\partial\boldsymbol{\Phi}/\partial\alpha_{l}$.

  2. 2.

Calculate the pseudo-inverse $\boldsymbol{\Phi}^{\dagger}$ to generate the
variable projection functional from (7), which is dependent on a given
$\boldsymbol{\alpha}$.

  3. 3.

Use the partial derivatives $\partial\boldsymbol{\Phi}/\partial\alpha_{l}$ to
generate the Jacobian matrix
$\frac{\mathrm{d}\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}}{\mathrm{d}\boldsymbol{\alpha}}\>\boldsymbol{y}$
for a given $\boldsymbol{\alpha}$, as in (9).

  4. 4.

Minimize the variable projection functional (8) by using the results from
steps 2 and 3 and insert them into an already existing nonlinear least squares
solver to obtain the final result $\hat{\boldsymbol{\alpha}}$.

  5. 5.

Take $\hat{\boldsymbol{\alpha}}$ and calculate a final
$\boldsymbol{\Phi}^{\dagger}(\hat{\boldsymbol{\alpha}})$ to solve the
remaining linear parameters $\hat{\boldsymbol{\beta}}$, as in (6).

In the original implementation by O’Leary and Rust (2013), the linear
solution, including the computation of
$\boldsymbol{\Phi}^{\dagger}(\boldsymbol{\alpha})$, was calculated via a
singular value decomposition, which is a very robust direct method for solving
linear least squares problems. Moreover, instead of using the simplification
by Kaufman (1975) of dropping the second term in the Jacobian (9), O’Leary and
Rust (2013) argued that, in modern computers, the balance between the
computing time for extra iterations vs. the second term has tipped. In 1975,
when Kaufman (1975) proposed simplification, it used to be more efficient to
put up with more functional evaluations by saving matrix computations; this
could have changed today (O’Leary and Rust, 2013). To verify this argument, a
simplification by Kaufman (1975) is included as a subject for testing in
Section 4.3.

To solve the nonlinear problem in step 4, the  least_squares function (SciPy
v1.8.0 Manual, 06.02.2022) from the scipy.optimize package (Virtanen et al.,
2020) is used in the implementation. Its default solving method is the so-
called trust region reflective algorithm (Branch et al., 1999), which has been
shown to work efficiently within the VP algorithm for the example discussed in
Section 4 (Bärligea, 2022).

###  3.2 Modifications for Multiple Right-Hand Sides

In the following, the modifications required for multiple right-hand sides
$\boldsymbol{y_{1}},\ldots,\boldsymbol{y_{s}}$ with possibly varying lengths
$m_{1},\ldots,m_{s}$, are going to be discussed. For the naive approach, the
only necessary change of the above procedure arises in the user-supplied
function ADA from step 1; instead of one matrix
$\boldsymbol{\Phi}(\boldsymbol{\alpha})$ and its derivative, it has to return
the sparse matrix $\boldsymbol{G}(\boldsymbol{\alpha})$ as in (22), with the
possibly different $\boldsymbol{\Phi}_{k}(\boldsymbol{\alpha})$ on its
diagonal, and the corresponding derivatives are also sparse matrices with the
Jacobians $\partial\boldsymbol{\Phi}_{k}/\partial\boldsymbol{\alpha}$ on the
diagonal. After setting up $\boldsymbol{\tilde{y}}$ and
$\boldsymbol{\tilde{\beta}}$, as in Equation (23), the above steps can be
performed unchanged to obtain the results $\hat{\boldsymbol{\alpha}}$ and
$\tilde{\boldsymbol{\beta}}=(\hat{\boldsymbol{\beta_{1}}},\ldots,\hat{\boldsymbol{\beta_{s}}})^{\top}$.

This is different from the Golub–LeVeque approach described in Section 2.2.
Here, the user-defined method ADA remains the same as in step 1 of Section
3.1; however, instead of calculating (8) in step 2, the new vector function
$\boldsymbol{z}(\boldsymbol{\alpha})$ has to be set up, such as (27). The same
holds for the new derivative matrix
$\mathrm{d}\boldsymbol{z}(\boldsymbol{\alpha})/\mathrm{d}\boldsymbol{\alpha}$
from (28) in step 3. After this, in step 4,
$\left\|\boldsymbol{z}(\boldsymbol{\alpha})\right\|^{2}$ is minimized for
$\hat{\boldsymbol{\alpha}}$; consequently, in step 5, the final linear
parameter matrix $\hat{\boldsymbol{B}}$ can be derived by calculating its
$k{\text{th}}$ column as follows:

$\boldsymbol{\beta_{k}}=\boldsymbol{\Phi_{k}}^{\dagger}(\hat{\boldsymbol{\alpha}})\>\boldsymbol{y_{k}}\
\mathrm{.}$ (29)

Finally, for the Kaufman approach (Kaufman, 1975), it is important to note
that the pseudo-inverse computed in steps 2 and 5 is now calculated via the QR
decomposition (10), as in (11), while the remaining matrix
$\boldsymbol{Q_{2}}$ of the decomposition can be taken for setting up
$\boldsymbol{z}(\boldsymbol{\alpha})$ and
$\partial\boldsymbol{z}(\boldsymbol{\alpha})/\partial\boldsymbol{\alpha}$,
just as before, by replacing $\boldsymbol{P}^{\perp}_{\boldsymbol{\Phi}}$ in
(27) and (28) with $\boldsymbol{Q_{2}}^{\top}$. The simplified Jacobian (20)
is also adopted in the implementation of this approach.

##  4 Numerical Experiments

The goal of this section is to outline an example of a separable least squares
problem with multiple right-hand sides and use it to evaluate the performance
of the suggested VP algorithms. The tests were also set up to establish a
comparison to conventional nonlinear least squares (NLS) methods, which ignore
separability. Those were implemented using the least_squares function from the
scipy.optimize package (SciPy v1.8.0 Manual, 06.02.2022).

Without separation, the minimization problem of multiple right-hand sides can
be composed as follows:

$\min_{\boldsymbol{x}}\ \left\|\begin{pmatrix}\boldsymbol{y_{1}}\\\ \vdots\\\
\boldsymbol{y_{s}}\end{pmatrix}-\begin{pmatrix}\boldsymbol{\eta}(\boldsymbol{\beta_{1}},\boldsymbol{\alpha})\\\
\vdots\\\
\boldsymbol{\eta}(\boldsymbol{\beta_{s}},\boldsymbol{\alpha})\end{pmatrix}\right\|^{2}\qquad\text{with}\
\boldsymbol{x}=(\boldsymbol{\alpha},\boldsymbol{\beta_{1}},\ldots,\boldsymbol{\beta_{s}})^{\top}\
\mathrm{.}$ (30)

which is easily extendable to the case of varying sizes of
$\boldsymbol{y_{1}},\ldots,\boldsymbol{y_{s}}$, and correspondingly different
$\boldsymbol{\eta_{1}},\ldots,\boldsymbol{\eta_{s}}$.

###  4.1 Application: Trace Gas Retrieval

A real-world example for the described problem set can arise in the area of
remote sensing, more specifically, in the retrieval of atmospheric trace gas
concentrations from spectral radiance measurements. The concentration of
carbon dioxide (CO2) or methane (CH4)—both important greenhouse gases—can be
inferred from spectra observed in the short-wave infrared (SWIR). Such
measurements are often spaceborne in order to achieve global coverage of the
atmospheric composition.

For this paper, observations from the OCO-2 (Orbiting Carbon Observatory-2)
satellite by NASA (Crisp et al., 2004; Eldering et al., 2017; Jet Propulsion
Laboratory, California Institute of Technology, 2023) were used, which was
designed to monitor CO2 by measuring its absorption bands in the SWIR. In this
spectral region, a radiance measurement can be modeled by the radiative
transfer model (based on the Beer–Lambert law for molecular absorption
neglecting scattering) (Gimeno García et al., 2011)

$\hat{I}(\nu,\boldsymbol{r},\boldsymbol{\alpha})=\
\left(\sum_{j=0}^{n-1}r_{j}\>\nu^{j}\right)\cdot\mu_{\odot}\cdot
I_{0}(\nu)\cdot\exp\left(-\sum_{l=1}^{p}\alpha_{l}\>\tau_{l}^{\text{prior}}(\nu)\right)\
\otimes\ S(\nu),$ (31)

dependent on the wavenumber $\nu$. The term wavenumber, which represents the
inverse of the often used wavelength $\lambda$, has the common unit for the
SWIR region of ${\mathrm{cm}}^{-1}$, corresponding to
$${10}^{4}$/$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$. The two sets of fitting
parameters are $\boldsymbol{r}\in\mathbb{R}^{n}$ for linear parameters $r_{j}$
and $\boldsymbol{\alpha}\in\mathbb{R}^{p}$ for nonlinear parameters
$\alpha_{l}$.

The first term is a polynomial that approximates the wavenumber-dependent
surface reflectivity of the Earth at the measurement location. Factor
$\mu_{\odot}$ corresponds to $\cos(\theta_{\odot})$ (with $\theta_{\odot}$:
solar zenith angle) and accounts for the geometry of the measurement,
$I_{0}(\nu)$ is the incoming solar radiation (at the top of the atmosphere),
and

$\tau_{l}^{\text{prior}}(\nu)=\int_{s}n_{l}(s)\ \sigma_{l}(\nu,p(s),T(s))\
\mathrm{d}s$ (32)

is the total optical depth of the $l{\text{th}}$ molecule, which is the path
integral over the number density $n_{l}$, and its pressure- and temperature-
dependent cross-section $\sigma_{l}$. In trace gas retrieval, $\tau(\nu)$ is
the most important measure, as it is directly related to a molecule’s
concentration in the atmosphere on a given path $s$. Unfortunately, SWIR
observations do not provide enough information to retrieve the concentration
profile of a molecule. Therefore, the calculation of (32) is only possible
under prior assumptions of the atmospheric state (i.e., the temperature and
pressure profiles $p(s)$, $T(s)$, and the molecular number density, i.e.,
$n(s)$). Hence, a simple scaling factor $\alpha_{l}$ is fitted as

$\tau_{l}(\nu)=\alpha_{l}\cdot\tau_{l}^{\text{prior}}(\nu)$ (33)

in the forward model (31) to retrieve the “real” optical depth at the time and
place of measurement. Lastly, $S(\nu)$ is the spectral response function of
the sensor, which has to be convolved with the monochromatic radiance in order
to mimic a real measurement.

For trace gas retrieval, one has to consider all $p$ molecules that have non-
negligible absorbance in the measured spectral region. In the case of OCO-2
observations, the only relevant molecule apart from CO2 is H2O (water vapor),
meaning that there are two nonlinear fitting parameters. For the linear
parameters, it is common to use approximately three reflectivity coefficients
(depending on the size of the spectral interval). This means that for each
spectrum, the necessary fitting parameters are

$\boldsymbol{r}=(\>r_{0},\>r_{1},\>r_{2}\>)\quad\text{and}\quad\boldsymbol{\alpha}=(\>\alpha_{\ce{CO2}},\>\alpha_{\ce{H2O}}\>)\
\mathrm{,}$ (34)

Note that even though it is physically necessary to use all of these variables
in a fit, the only one of interest in this context is the molecular scaling
factor $\alpha_{l}$ of the molecule $l$ under scrutiny, which in this case is
$\alpha_{\ce{CO2}}$), as this alone contains the relevant information about
its atmospheric concentration.

This together with (31) clearly fits the criteria for a separable problem, and
a conventional VP algorithm from the PORT Mathematical Subroutine Library
(Dennis et al., 1981; Fox et al., 1978) has already been tested by Gimeno
García et al. (2011) and validated by Hochstaffl et al. (2018).

How is this an example of problems with multiple right-handed sides? Many
satellites have sensors that measure radiance simultaneously in several
spectral windows, e.g., OCO-2 observes the strong (around
$6250\text{\,}{\mathrm{cm}}^{-1}$) and weak absorption bands (around
$5000\text{\,}{\mathrm{cm}}^{-1}$) of CO2 (cf. Figure 1). Assuming consistent
model input, both spectra should deliver the same values for
$\boldsymbol{\alpha}$, but as surface reflectivity varies strongly for
different wavelength regions, they each have a specific reflectivity
polynomial and, therefore, $\boldsymbol{r}$. Thus, for every observation, two
spectral measurement windows (of different lengths) should be fitted
simultaneously.

This concept of multiple right-hand sides can also be transferred into a
spatial dimension: Some molecules, such as carbon dioxide or methane, are very
long-lived, so they are distributed relatively homogeneously in the
atmosphere. This means that observations from nearby locations should all
yield quite similar concentrations. Thus, they might as well be fit for one
$\alpha_{\ce{CO2}}$ at once. Note that the assumption made about atmospheric
carbon dioxide might not hold for all other absorbing molecules in the
observed spectral region, such as H2O, which has rather variable
concentrations across the globe. However, as their variations are less than
that of surface reflectivity, and no physical insight is sought from the fit
of the $\alpha_{\ce{H2O}}$ parameter, it can be seen as a mere auxiliary
parameter for completing the model and, therefore, be treated as a “constant”
nonlinear fitting variable for a group of neighboring spectra. Still, in this
case, the reflectivity coefficients, $r_{j}$, representative of the surface at
the place of measurement, are distinct for every geolocation and, therefore,
specific to each measured spectrum. Another possible linear model parameter,
which is distinct for each spectrum, would be a constant baseline correction
added to the model (31), as suggested by Gimeno García et al. (2011).

Figure 1: Four exemplary soundings of frame $1728$ from the OCO-2 level 1b
(L1B) measurement product (Crisp et al., 2021), each displaying radiance
spectra in units
[$$\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{/}\mathrm{(}\mathrm{s}\;\mathrm{c}\mathrm{m}\mathrm{{}^{2}}\;\mathrm{sr}\;\mathrm{c}\mathrm{m}^{-1}$)$]
from both the strong and weak bands with $809$ and $651$ spectral pixels each.

Finally, it is possible to combine both the spectral and spatial dimensions of
multiple RHS fittings in trace gas retrieval. The OCO-2 satellite, for
instance, always stores eight observations (“soundings”) in one so-called
“frame”, with the spatial coverage not exceeding
$24\text{\,}\mathrm{k}\mathrm{m}^{2}$. The concentration of carbon dioxide can
be assumed to fluctuate only minimally within such an area on the globe.
Figure 2 shows exemplary retrieval results of eight soundings from one OCO-2
frame. Most of the fluctuations in these results (all except for sounding
number 5) could be merely due to noise in the measurements, as the mean value
stays within the uncertainty for almost all of them. A fit over multiple
observations, as proposed, could, therefore, help to constrain the fluctuation
to a more reliable retrieval product.

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 (Gunson and Eldering, 2019a;
Crisp et al., 2020) for one exemplary frame (cf. Figure 1), including eight
soundings.

To summarize, the OCO-2 product (Crisp et al., 2021) allows for multiple RHS
fits of a combined $8$ (spatial) $\cdot$ $2$ (spectral) $=$ $16$ datasets (cf.
Section 4.3).

The following tests were conducted with a Python version of BIRRA (Beer
infrared retrieval algorithm) (Gimeno García et al., 2011), which has been
validated for the SWIR trace gas retrieval of CO by Hochstaffl et al. (2018).
In this code, the Jacobian matrix of the model function (31) is set up
analytically for the least squares fit, reducing numerical instabilities.
BIRRA is an extension of the radiative transfer model Py4CAtS (Python for
Computational Atmospheric Spectroscopy) (Schreier et al., 2019), which is used
to calculate the a priori total optical depths (32) needed in (31).

It must be noted that all retrievals conducted with the model described above
are only supposed to evaluate the methodology and algorithms and in no way
claim to represent full-fledged physical CO2 retrieval products, such as those
by Crisp et al. (2020). Moreover, this technique of fitting multiple spectra
measured within a certain spatial distance from each other is only reasonable
as long as the assumption holds that, within an order of magnitude of the
spatial resolution of the sensor, there are only small to no physically caused
fluctuations/gradients in the sought trace gas concentration(s). This is, of
course, not the case for localized emitters, such as power plants or biomass-
burning events.

###  4.2 Tests with Synthetic Data

The goal of this subsection is to show the conceptual and effective
differences of algorithms solving multiple RHS compared to the classical case
of solving one. This analysis was conducted on the basis of synthetic spectra.
Those are simulated radiance measurements generated with the radiative
transfer model Py4CAtS (Schreier et al., 2019). The benefit of using this in
tests is that, in the retrieval (i.e., the fitting process), there is no model
error and the exact solution is known. The only deviation from a “perfect” fit
is, therefore, controlled by adding noise to the modeled measurements.

In order to be representative of the later tests with real measurements, the
same numbers, sizes, and types (distinct in spatial or spectral dimensions) of
datasets were generated as the ones used by OCO-2. Moreover, for consistency,
all test retrievals were conducted using the same number of fitting
parameters, with $n=$3$$ linear ones per dataset and a total of $p=$2$$
nonlinear ones (cf. vectors in (34)). This allowed for the test cases
summarized in Table 1.

Table 1: Test cases for the analysis of fitting $1$ (denoted as single) or $16$ (denoted as MRHS) synthetically generated datasets, simultaneously, with both VP and NLS algorithms. Name |  #RHS |  #fits |  #Parameters per Fit |  Tested Algorithm Types
---|---|---|---|---
single |  1 |  16 |  $(1\cdot 3)+2=5$ |  VP & NLS
MRHS |  16 |  1 |  $(16\cdot 3)+2=50$ |  VP

For the MRHS case, the Golub–LeVeque VP algorithm (introduced in Section 2.2)
was used as a representative solver for multiple RHS problems. Tests with
synthetic spectra indicated that all mentioned MRHS methods (see Table 2)
yielded equal accuracy, confirming the theoretical proof offered by Golub and
Pereyra (1973) that the solutions found by a variable projection solver should
be equivalent to those of conventional nonlinear solving methods.

Table 2: Test cases for the analysis of real radiance spectra with both VP and NLS algorithms for multiple RHS. Name |  Algorithm |  Reference
---|---|---
VP GL |  Golub–LeVeque approach |  Section 2.2.2, Golub and LeVeque (1979)
VP KM |  Kaufman approach |  Section 2.2.3, Kaufman (1975)
VP naive |  Naive approach |  Section 2.2.1, Golub and LeVeque (1979)
NLS TRF |  Trust Region Reflective |  Branch et al. (1999)
NLS LM |  Levenberg–Marquardt |  Moré (1978)

For the single cases, a classical VP (based on O’Leary and Rust (2013)) and a
conventional NLS single-RHS solver (SciPy v1.8.0 Manual, 06.02.2022) (based on
Branch et al. (1999)) were tested.

First, the fitting precision of the VP MRHS solver was compared to that of
single solvers using spectra with signal-to-noise (SNR) ratios in the range of
$20$ to $500$. One measure of the goodness of fitting results is the relative
error,

$\epsilon=\frac{\alpha_{\text{true}}-\alpha_{\text{fit}}}{\alpha_{\text{true}}},$
(35)

compared to the true parameter values.

Figure 3 shows the distribution of these errors and the corresponding standard
deviations for different SNR values. The signal-to-noise ratios achieved by
satellites, including OCO-2, ranged between approximately $200$ and $800$ for
the frames used (Crisp et al., 2021). This broad variation of OCO-2’s SNR
comes from changes in the solar position and varying surface reflectivities
across the orbit. As expected, both solvers achieved improved precision for
increasing SNR values, since a fit becomes more accurate for less noisy data.
While both single methods (NLS and VP) showed equal performances, the VP MRHS
yielded standard deviations that were slightly worse. This trend is also
reflected in the distributions of the relative errors, which are always more
sharply distributed around zero for the single solvers than for the MRHS
solver. This behavior is not surprising since a less-dimensional residual
vector (coming from the shorter data vector) and fewer unknown parameters most
generally leave less freedom in the fit and, therefore, lead to more precise
fitting results. To phrase it differently: one can expect that, as the size of
the least squares problem increases, the number of possible local minima that
the fit can reach will behave accordingly.

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).

Considering the similarly shaped distributions in Figure 3 and the fact that
VP MRHS seems to improve at the same rate as NLS and VP single, MRHS fits can
be viewed as equally effective. In particular, at higher SNRs, both
methodologies achieve deviations from the exact results in such low orders of
magnitude that the precisions of their fitted $\alpha_{\ce{CO2}}$ values are
very comparable.

Another important measure for the accuracy of a fit is the standard deviation
of the residuals $r_{i}=y_{i}-\hat{y}_{i}$ (prediction errors with
$\boldsymbol{y}$: vector of observations, and $\boldsymbol{\hat{y}}$: fitted
model), also known as the sigma of regression, defined as

$\sigma=\frac{\left\|\boldsymbol{r}\right\|}{\sqrt{m-n-p}},$ (36)

where the numerator is the norm of the residual vector of the fit and the
denominator represents the number of degrees of freedom (the number of data
points minus the number of unknown variables). Of course, for the $s$ right-
hand side case, the number of linear parameters $n$ becomes $sn$, and the
number of data points becomes $ms$ (or $m_{1}+\ldots+m_{s}$), respectively.

Figure 4 shows the mean sigma of regression of all fits over the SNR. Here,
the development of single and MRHS is similar to that of the errors discussed
above. The mean sigmas produced by the MRHS fits are slightly higher than the
ones achieved by the single solvers. Again, this is intuitive: A single solver
is able to produce distinct nonlinear parameters for the noisy spectra and,
therefore, has more degrees of freedom to mimic the noisy spectra. An MRHS
solver, on the other hand, only has one set of nonlinear parameters for all
the spectra, leading to an overall larger deviation between the “observed” and
modeled data. This does not necessarily mean the latter is less accurate. On
the contrary, since it is less prone to including the specific noise of the
spectra into the fit, MRHS solvers could have a smoothing effect on otherwise
fluctuating retrieval results (see Figure 2).

Figure 4: Mean sigma of regression for both single and MRHS fits for different
noisy spectra. The dashed lines correspond to fitted hyperbolas.

The effect of possible “overfitting” might, however, only be relevant for high
noise levels (low SNR). As the sigma of regression is proportional to the
noisy radiance, which is proportional to
$(\boldsymbol{$1$}+\gamma/\text{SNR})$ (with $\boldsymbol{\gamma}$: normally
distributed random value), it decreases by the inverse of the SNR value for
both the single and MRHS fits (see fitted hyperbolas in Figure 4). In this
way, as the noise increases, the sigmas of regression become more similar,
such that for an SNR of $200$ and higher (representative of OCO-2), the
performance differences between the methodologies (single and MRHS) disappear.
Thus, for reasonably good data, there is no sacrifice in the precision or
accuracy of the produced results when fitting multiple RHS simultaneously
instead of fitting one by one.

Now that the differences between MRHS solvers and classical single solvers are
established, we need to analyze which MRHS algorithm (see Section 2.2) is the
best.

###  4.3 Tests with Real Measurements

The goal of this subsection is to assess the performance of the new enhanced
VP algorithms for multiple RHS described in Sections 2.3 and 3.2 (VP naive, VP
Golub–LeVeque, VP Kaufman) by comparing them to conventional NLS solvers.

The SciPy function used for the NLS reference approach allows the user to
choose between three different nonlinear least squares algorithms (explored
thoroughly for a single RHS VP algorithm by Bärligea (2022)). In order to
better judge the solvers’ performances, two of them were used in the tests:
the trust region reflective method (’TRF’) (Branch et al., 1999) and the
Levenberg–Marquardt method (’LM’) (Moré, 1978). While the former is the most
efficient, the latter can be considered the most robust (Bärligea, 2022),
which could be helpful for an increasing number of variables.

In this subsection, an analysis of the test cases listed in Table 2 is
conducted. For the assessment, a set of $18$ OCO-2 frames (all measured on the
25 of May 2020 on orbit 31366a in the nadir (downward view) acquisition mode
just above Australia, with a spacecraft altitude of approximately
$711\text{\,}\mathrm{km}$ (Gunson and Eldering, 2019b)) was used; each
included 8 observations in both spectral bands (cf. Figure 1), within an area
of $24\text{\,}\mathrm{k}\mathrm{m}^{2}$ measured along a ground track no
wider than $80\text{\,}\mathrm{km}$, labeled as cloud-free (no scattering),
above land (better reflectivity), and good quality, according to criteria
defined by Crisp et al. (2021). With those, a few hundred test fits were
performed, with the VP and NLS methods (see Table 2, with varying numbers of
RHS ranging from $2$ to the maximum available number of $16$ (only even
numbers due to the combination of the two spectral bands in one observation).
Again, the fits used $n=$3$$ linear parameters per spectrum and $p=$2$$
nonlinear parameters.

For the evaluation of accuracy, the sigma of regression, the R-Score measure,
the confidence bounds of the results, and the fitted residuals, were analyzed.
The sigma of regression $\sigma$ defined in Equation (36) turned out to be
equal for all of the tested methods (see the residual analysis below).

A second statistical quantity is the so-called R-score, defined as

$R=\
\frac{\sum_{i=1}^{M}(\hat{y}_{i}-\overline{\boldsymbol{y}})^{2}}{\sum_{i=1}^{M}(y_{i}-\overline{\boldsymbol{y}})^{2}},$
(37)

indicating the amount of variance (the mean of the measurements
$\overline{\boldsymbol{y}}$) accounted for by the fitted model
$\hat{\boldsymbol{y}}$. $M$ means the cumulative size $m_{1}+\ldots+m_{s}$ of
all the datasets. $R$ must be within $[0,1]$, and the best possible score a
fit can achieve would be $1$. In the experiments, all of the discussed methods
obtained R-scores of approximately $0.99$. The only difference could be
observed for VP GL and VP KM, which had average higher scores by
$0.02\text{\,}\mathrm{\char 37\relax}$ compared to all other methods, which is
negligible.

In order to calculate the confidence bounds of the retrieval results, the
covariance matrix

$\boldsymbol{C}=\sigma^{2}\cdot(\boldsymbol{H}^{\top}\boldsymbol{H})^{-1}\
\in\ \mathbb{R}^{(p+ns)\times(p+ns)}$ (38)

needs to be calculated, with $\boldsymbol{H}$ containing the partial
derivatives of the model function, with respect to the $p$ nonlinear and $sn$
linear parameters. For a VP method with multiple right-hand sides, it can be
composed as follows:

$\boldsymbol{H}=\begin{pmatrix}\frac{\mathrm{d}\boldsymbol{z}(\boldsymbol{\alpha})}{\mathrm{d}\boldsymbol{\alpha}}\big{|}_{\boldsymbol{\alpha}=\boldsymbol{\hat{\alpha}}}&\boldsymbol{G}(\boldsymbol{\hat{\alpha}})\end{pmatrix}\
\in\ \mathbb{R}^{M\times(p+ns)}\mathrm{.}$ (39)

Here, the first matrix is the $M\times p$ Jacobian of the purely nonlinear
function $\boldsymbol{z}(\boldsymbol{\alpha})$ defined in (26) and (27) with
respect to the nonlinear parameters $\boldsymbol{\alpha}$, and the second
matrix $\boldsymbol{G}(\boldsymbol{\hat{\alpha}})\in\mathbb{R}^{M\times ns}$,
defined in (22), is the Jacobian of all the linear parameters
$\boldsymbol{\beta}$ (cf. Equation (23)). For a confidence level of
$95\text{\,}\mathrm{\char 37\relax}$, one can then calculate the confidence
bound(s) (CB) of the retrieved parameters $\hat{\boldsymbol{x}}$ by

$\hat{\boldsymbol{x}}_{\text{CB}}=q\cdot\sqrt{\text{diag}(\boldsymbol{C})},$
(40)

for which $q$ represents the standard normal distribution quantile of
$\frac{1}{2}(1-0.95)$, and the diagonal elements of $\boldsymbol{C}$ are the
variances of the estimated parameters $\hat{\boldsymbol{x}}$.

Figure 5 shows the distribution and mean values of the calculated confidence
bounds for the $\alpha_{\ce{CO2}}$ parameter for an increasing number of RHS.
While the confidence bounds are already relatively small, they are decreasing
for an increasing number of datasets, similar to before with the increasing
SNR values (see Figure 4). This indicates that more data cause more accurate
fitting results. However, a small difference can be observed in Figure 5
between the naive VP method and the “others” (including VP GL, VP KM, NLS TRF,
and NLS LM, which all produced the same results). Apparently, the confidence
bounds of the results from the naive method, though decreasing, are slightly
worse than the rest. This is probably due to the different and more lavishly
calculated Jacobian matrix of the naive problem (23). One can, therefore,
argue that this is mainly a numerical issue and does not correspond to a lack
of accuracy of the VP naive solver. In light of the measures considered above,
the tested MRHS solvers all achieved equally accurate fits.

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.

This was also confirmed when the residuals of the fits were analyzed, which
turned out equally for all methods (cf. Table 2). The statistical diagnostics
for the residuals of one exemplary VP GL fit (representative of all methods,
including NLS) are shown in Figure 6.  Ideally, the errors between the fitted
model and measurements should be normally distributed. Due to noise and
outliers in the spectral data, this distribution may, however, deviate
slightly from a normal one. Yet, the fact that the residuals have their
highest density around zero indicates that all the algorithms conducted
reasonably good fits.

As for the robustness, all algorithms yielded convergence rates of
$100\text{\,}\mathrm{\char 37\relax}$ for decent initial guesses. For a
discussion on the impact of bad initial guesses, see O’Leary and Rust (2013),
who showed that the VP method ultimately converges more reliably than
conventional NLS algorithms. This is mostly due to the fact that the former
deal with a reduced nonlinear least squares problem needing only $p$ instead
of $p+n$ initial guesses, making the solver a lot more stable.

Figure 6: Statistical diagnostics for one exemplary VP fit. The plot on the
left shows the distribution of the residuals via a coarse histogram and the
continuous kernel density estimate (KDE). The plot on the right is a normal
probability plot, displaying the deviation of the residuals from a normal
distribution.

In the next step, the fitting times of all mentioned methods were analyzed to
compare their computational efficiency. Figure 7 shows the mean running times
for a fit for an increasing number of datasets. Here, the VP KM method is not
shown since its performance is similar to VP GL. For fairly small numbers of
datasets, the NLS algorithms were faster than the VP methods. This stems from
the fact that these algorithms are part of the SciPy package (Virtanen et al.,
2020), which is operationally optimized, whereas the proposed VP code was
originally made as a proof of concept and is not yet optimized in the same
manner. Still, this scheme changed drastically when more RHS were used. Table
3 shows the exact values for $2$, $4$, and $6$ datasets.

Figure 7: Comparison of the evolution of mean running times of NLS and VP methods for a fit over the number of used datasets. Table 3: Fitting times for VP methods compared to NLS methods for small numbers of datasets. The bold numbers mark the smallest mean values within a row. #RHS |  VP GL [s] |  VP Naive [s] |  NLS TRF [s] |  NLS LM [s]
---|---|---|---|---
2 |  3.97 |  4.3 |  1.95 |  2.46
4 |  7.83 |  9.18 |  6.40 |  8.02
6 |  11.59 |  14.78 |  13.51 |  17.0

For six RHS and more, the suggested VP GL algorithm not only becomes
significantly faster than the rest, but it is also the only method with
fitting times that increase linearly with the number of right-hand sides (see
Figure 7), while all the other tested methods exhibit an almost quadratic
evolution. This confirms that VP GL and VP KM are the most efficient methods
when it comes to dealing with the rising complexity of multiple RHS problems.
It also reveals the inferiority of the naive VP method compared to the ’good’
VP methods in every test. Even though the naive approach separates the problem
and should, therefore, be just as stable as the good approaches, the time it
needs for solving also rises quadratically with the number of fitting windows,
similar to the slower (and inferior) NLS solvers. This must be due to the
increasing size of the block diagonal matrix
$\boldsymbol{G}(\boldsymbol{\alpha})$ and the resulting extra costs for
calculating $\boldsymbol{G}^{\dagger}(\boldsymbol{\alpha})$ and the Jacobian
$\partial\boldsymbol{P}^{\perp}_{\boldsymbol{G}}/\partial\boldsymbol{\alpha}$.

Comparisons of the two “good” VP algorithms (VP GL and VP KM) showed that, in
all of the above categories, such as robustness or accuracy, the Kaufman
approach did equally as well as the Golub–LeVeque one. The only difference
could be found in the fitting times, for which the method by Kaufman (1975),
as predicted, was consistently faster. However, the relative improvements in
the running times remained below $1\text{\,}\mathrm{\char 37\relax}$ and are,
therefore, almost negligible. This confirms the point made by O’Leary and Rust
(2013) that Kaufman’s simplification does not necessarily pose a computational
benefit to modern computers anymore.

##  5 Discussion and Conclusions

Motivated by a real-world application in atmospheric remote sensing, a
variable projection algorithm was extended to multiple right-hand sides of
different sizes and nonlinear model setups. A modern MATLAB implementation by
O’Leary and Rust (2013) was translated into Python and modified according to
the theory presented in this paper. It incorporates the ideas of Golub and
LeVeque (1979) and Kaufman (1975) for solving separable nonlinear least
squares problems with multiple RHS.

Numerical tests using synthetic data demonstrate that simultaneous fittings
over multiple measurements maintain accuracy and precision compared to single
dataset solvers, with potential benefits in reducing “overfitting” (with noise
and outliers affecting the retrieval results) and fluctuations in the results.

A comprehensive comparison with conventional nonlinear least squares solvers
using real measurements from NASA’s OCO-2 satellite (Jet Propulsion
Laboratory, California Institute of Technology, 2023) indicated similar
accuracy among all algorithms. The most significant finding was that the
variable projection methods based on Golub and LeVeque (1979) and Kaufman
(1975) significantly outperformed all other methods in computing time,
particularly as the number of datasets increased. Thus, these algorithms are
deemed more efficient than conventional solvers. Furthermore, our experiments
indicate that a popular simplification proposed by Kaufman (1975) did not
yield significant performance improvements.

The algorithm presented in this article proved to be highly effective and
efficient. This indicates that the recommended modifications to the original
algorithm by Golub and LeVeque (1979) preserve its computational advantages.
Note, that the benefits arising from a fast solver must always be considered
in relation to the computational costs associated with the remaining part of
the overall task. In trace gas retrieval applications, the computation time
required for the forward model, i.e., radiative transfer Equation (31), can
significantly exceed that of solving the inherent least squares problem. Our
algorithm, thus, offers the most significant advantage in tasks where the
overall performance heavily relies on the fitting process. Consequently, we
endorse using this implementation not only for remote sensing, but also for
other scientific problems implemented in Python with similar characteristics.
The simultaneous fitting of more data can reduce fluctuations in the results,
which is highly desirable in some applications.

## Appendix

In this appendix, derivations for Equations (9) and (20) are presented (cf.
proof by Golub and Pereyra (1973)). Beginning with the full Jacobian (9): The
generalized inverse
$\boldsymbol{\Phi}^{\dagger}=(\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi})^{-1}\boldsymbol{\Phi}^{\top}$
satisfies the identities

$\displaystyle\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}\boldsymbol{\Phi}$
$\displaystyle=\boldsymbol{\Phi},$ (41)
$\displaystyle(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger})^{\top}$
$\displaystyle=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}.$ (42)

Note that the second identity follows from the fact that the inverse of the
symmetric matrix $\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}$ is also symmetric.
From Equation (42), we infer that the projector
$\boldsymbol{P}_{\boldsymbol{\Phi}}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}$
is symmetric; from Equation (41), we infer that

$\displaystyle\boldsymbol{P}_{\boldsymbol{\Phi}}\boldsymbol{\Phi}$
$\displaystyle=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}\boldsymbol{\Phi}=\boldsymbol{\Phi},$
(43) $\displaystyle\boldsymbol{P}_{\boldsymbol{\Phi}}^{2}$
$\displaystyle=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}=\boldsymbol{P}_{\boldsymbol{\Phi}}.$
(44)

By means of Equation (43), we find

$\partial_{l}(\boldsymbol{P}_{\boldsymbol{\Phi}}\boldsymbol{\Phi})=(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{\Phi}+\boldsymbol{P}_{\boldsymbol{\Phi}}(\partial_{l}\boldsymbol{\Phi})=(\partial_{l}\boldsymbol{\Phi}),$
(45)

implying

$(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{\Phi}=(\partial_{l}\boldsymbol{\Phi})-\boldsymbol{P}_{\boldsymbol{\Phi}}(\partial_{l}\boldsymbol{\Phi})=(\boldsymbol{I}-\boldsymbol{P}_{\boldsymbol{\Phi}})(\partial_{l}\boldsymbol{\Phi})=\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}(\partial_{l}\boldsymbol{\Phi}),$
(46)

with the short-hand notation $\partial_{l}=\partial/\partial\alpha_{l}$. Then,
using Equation (46) and the relation
$\boldsymbol{P}_{\boldsymbol{\Phi}}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}$,
we compute the quantity
$(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{P}_{\boldsymbol{\Phi}}$
as

$(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{P}_{\boldsymbol{\Phi}}=(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}=\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}(\partial_{l}\boldsymbol{\Phi})\boldsymbol{\Phi}^{\dagger}.$
(47)

On the other hand, because $\boldsymbol{P}_{\boldsymbol{\Phi}}$ is symmetric,
$\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}}$ is also symmetric. Therefore,
we have

$[(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{P}_{\boldsymbol{\Phi}}]^{\top}=\boldsymbol{P}_{\boldsymbol{\Phi}}(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}}),$
(48)

further, in view of Equation (47),

$\boldsymbol{P}_{\boldsymbol{\Phi}}(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})=[\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}(\partial_{l}\boldsymbol{\Phi})\boldsymbol{\Phi}^{\dagger}]^{\top}.$
(49)

Finally, (44), (47), and (49) yield

$(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})=\partial_{l}(\boldsymbol{P}_{\boldsymbol{\Phi}}^{2})=(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})\boldsymbol{P}_{\boldsymbol{\Phi}}+\boldsymbol{P}_{\boldsymbol{\Phi}}(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})=\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}(\partial_{l}\boldsymbol{\Phi})\boldsymbol{\Phi}^{\dagger}+[\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}(\partial_{l}\boldsymbol{\Phi})\boldsymbol{\Phi}^{\dagger}]^{\top},$
(50)

and the relation
$\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}}^{\perp}=-\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}}$
can be used to conclude formula (9) presented in Section 2.1.

Regarding Equation (20): Using the orthogonality relation
$\boldsymbol{Q}_{2}^{\top}\boldsymbol{Q}_{1}=\boldsymbol{0}_{(m-n)\times n}$,
we find

$\boldsymbol{Q}_{2}^{\top}\boldsymbol{\Phi}=\boldsymbol{Q}_{2}^{\top}\boldsymbol{Q}_{1}\boldsymbol{R}_{1}=\boldsymbol{0}_{(m-n)\times
n},$ (51)

so that by differentiating, we have

$\partial_{l}(\boldsymbol{Q}_{2}^{\top})\boldsymbol{\Phi}=-\boldsymbol{Q}_{2}^{\top}(\partial_{l}\boldsymbol{\Phi}).$
(52)

Multiplying the above Equation from the right by $\boldsymbol{\Phi}^{\dagger}$
and using the relation
$\boldsymbol{P}_{\boldsymbol{\Phi}}=\boldsymbol{\Phi}\boldsymbol{\Phi}^{\dagger}$,
we obtain

$\partial_{l}(\boldsymbol{Q}_{2}^{\top})\boldsymbol{P}_{\boldsymbol{\Phi}}=-\boldsymbol{Q}_{2}^{\top}(\partial_{l}\boldsymbol{\Phi})\boldsymbol{\Phi}^{\dagger}.$
(53)

Now, using the identity

$\boldsymbol{Q}_{2}^{\top}\boldsymbol{P}_{\boldsymbol{\Phi}}=\boldsymbol{Q}_{2}^{\top}(\begin{array}[]{cc}\boldsymbol{Q}_{1}&\boldsymbol{Q}_{2}\end{array})\left(\begin{array}[]{c}\boldsymbol{Q}_{1}^{\mathrm{\top}}\\\
\boldsymbol{Q}_{2}^{\top}\end{array}\right)=(\begin{array}[]{cc}\boldsymbol{0}_{(m-n)\times
n}&\boldsymbol{I}_{(m-n)\times(m-n)}\end{array})\left(\begin{array}[]{c}\boldsymbol{Q}_{1}^{\mathrm{\top}}\\\
\boldsymbol{Q}_{2}^{\top}\end{array}\right)=\boldsymbol{Q}_{2}^{\top},$ (54)

yields

$\partial_{l}(\boldsymbol{Q}_{2}^{\top}\boldsymbol{P}_{\boldsymbol{\Phi}})=\partial_{l}(\boldsymbol{Q}_{2}^{\top})\boldsymbol{P}_{\boldsymbol{\Phi}}+\boldsymbol{Q}_{2}^{\top}(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}})=\partial_{l}(\boldsymbol{Q}_{2}^{\top}).$
(55)

We infer that (cf. Equation (53))

$\partial_{l}(\boldsymbol{Q}_{2}^{\top})=-\boldsymbol{Q}_{2}^{\top}(\partial_{l}\boldsymbol{\Phi})\boldsymbol{\Phi}^{\dagger}+\boldsymbol{Q}_{2}^{\top}(\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}}),$
(56)

where $\partial_{l}\boldsymbol{P}_{\boldsymbol{\Phi}}$ is given by (50).
However, Kaufman (1975) suggested and mathematically justified neglecting the
second term in the above Equation. The resulting approximation is Equation
(20), which works sufficiently well in practice.

## Data Availability

The Python code for the presented variable projection algorithm, as well as
examples, are available at https://atmos.eoc.dlr.de/tools/varPro/index.html
(accessed on 16 May 2023). All data used for the tests with synthetic spectra
can be reproduced with the radiative transfer model Py4CAtS
(https://atmos.eoc.dlr.de/tools/Py4CAtS/ (accessed on 16 May 2023)) (Schreier
et al., 2019). Data concerning NASA’s OCO-2 satellite are publically available
at https://disc.gsfc.nasa.gov/ (accessed on 16 May 2023), provided by the
Goddard Earth Sciences Data and Information Services Center (GES DISC) (Gunson
and Eldering, 2019b, a). Any data arising from the numerical experiments that
were used in the evaluations will be made available upon request.

## References

  * Moré [1978] J. Moré.  The Levenberg-Marquardt algorithm: Implementation and theory.  In _Numerical Analysis_ , pages 105–116. Springer, Berlin, Heidelberg, 1978.  doi:10.1007/BFb0067700.
  * Golub and Pereyra [1973] G. Golub and V. Pereyra.  The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate.  _SIAM Journal on Numerical Analysis_ , 10(2):413–432, 1973.  doi:10.1137/0710036.
  * Golub and LeVeque [1979] G. Golub and R. LeVeque.  Extensions and uses of the variable projection algorithm for solving nonlinear least squares problems.  In _Proceedings of the 1979 Army Numerical Analysis and Computers Conference_ , Stanford, CA, 1979.
  * Hochstaffl et al. [2022] P. Hochstaffl, F. Schreier, C. H. Köhler, A. Baumgartner, and D. Cerra.  Methane retrieval from airborne HySpex observations in the short-wave infrared.  _Atmospheric Measurement Techniques Discussions_ , pages 1–37, 2022.  doi:10.5194/amt-2022-271.
  * Golub and Pereyra [2003] G. Golub and V. Pereyra.  Separable nonlinear least squares: the variable projection method and its applications.  _Inverse Problems_ , 19(2):R1–R26, 2003.  doi:10.1088/0266-5611/19/2/201.
  * O’Leary and Rust [2013] D. O’Leary and B. Rust.  Variable projection for nonlinear least squares problems.  _Computational Optimization and Applications_ , 54(3):579–593, 2013.  doi:10.1007/s10589-012-9492-9.
  * Ruhe and Wedin [1980] A. Ruhe and P. Wedin.  Algorithms for separable nonlinear least squares problems.  _SIAM Review_ , 22(3):318–337, 1980.  doi:10.1137/1022057.
  * Kaufman [1975] L. Kaufman.  A variable projection method for solving separable nonlinear least squares problems.  _BIT Numerical Mathematics_ , 15(1):49–57, 1975.  doi:10.1007/BF01932995.
  * Kaufman and Sylvester [1992] L. Kaufman and G. Sylvester.  Separable nonlinear least squares with multiple right-hand sides.  _SIAM Journal on Matrix Analysis and Applications_ , 13(1):68–89, 1992.  doi:10.1137/0613008.
  * Kaufman [2010] L. Kaufman.  Solving separable nonlinear least squares problems with multiple datasets.  In V. Pereyra and Godela Scherer, editors, _Exponential data fitting and its applications_ , pages 94–109. Bentham eBooks, 2010.  doi:10.2174/978160805048211001010094.
  * Mullen and van Stokkum [2007] K. Mullen and I. van Stokkum.  TIMP: An R package for modeling multi-way spectroscopic measurements.  _Journal of Statistical Software_ , 18(3):1–46, 2007.  doi:10.18637/jss.v018.i03.
  * LeVeque [1985] R. LeVeque.  SUBROUTINE VARP2, 1985.  URL http://www.netlib.org/opt/varp2.
  * Bolstad [1977] J. Bolstad.  SUBROUTINE VARPRO, 1977.  URL http://www.netlib.org/opt/varpro.
  * Fox et al. [1978] P. Fox, A. Hall, and N. Schryer.  The PORT Mathematical Subroutine Library.  _ACM Transactions on Mathematical Software_ , 4(2):104–126, 1978.  doi:10.1145/355780.355783.
  * SciPy v1.8.0 Manual [06.02.2022] SciPy v1.8.0 Manual.  scipy.optimize.least_squares, 06.02.2022.  URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html.
  * Virtanen et al. [2020] P. Virtanen et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi:10.1038/s41592-019-0686-2.
  * Branch et al. [1999] M. Branch, T. Coleman, and Y. Li.  A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems.  _SIAM Journal on Scientific Computing_ , 21(1):1–23, 1999.  doi:10.1137/S1064827595289108.
  * Bärligea [2022] A. Bärligea.  Assessment of a variable projection algorithm for trace gas retrieval in the short-wave infrared.  Bachelor’s thesis, Technical University of Munich, 2022.
  * Crisp et al. [2004] D. Crisp, R. Atlas, F. Breon, L. Brown, J. Burrows, P. Ciais, B. Connor, S. Doney, I. Fung, D. Jacob, C. Miller, D. O’Brien, S. Pawson, J. Randerson, P. Rayner, R. Salawitch, S. Sander, B. Sen, G. Stephens, P. Tans, G. Toon, P. Wennberg, S. Wofsy, Y. Yung, Z. Kuang, B. Chudasama, G. Sprague, B. Weiss, R. Pollock, D. Kenyon, and S. Schroll.  The Orbiting Carbon Observatory (OCO) mission.  _Advances in Space Research_ , 34(4):700–709, 2004.  doi:10.1016/j.asr.2003.08.062.
  * Eldering et al. [2017] A. Eldering, C. O’Dell, P. Wennberg, D. Crisp, M. Gunson, C. Viatte, C. Avis, A. Braverman, R. Castano, A. Chang, L. Chapsky, C. Cheng, B. Connor, L. Dang, G. Doran, B. Fisher, C. Frankenberg, D. Fu, R. Granat, J. Hobbs, R. Lee, L. Mandrake, J. McDuffie, C. Miller, V. Myers, V. Natraj, D. O’Brien, G. Osterman, F. Oyafuso, V. Payne, H. Pollock, I. Polonsky, C. Roehl, R. Rosenberg, F. Schwandner, M. Smyth, V. Tang, T. Taylor, C. To, D. Wunch, and J. Yoshimizu.  The Orbiting Carbon Observatory-2: first 18 months of science data products.  _Atmospheric Measurement Techniques_ , 10(2):549–563, 2017.  doi:10.5194/amt-10-549-2017.
  * Jet Propulsion Laboratory, California Institute of Technology [2023] Jet Propulsion Laboratory, California Institute of Technology.  Orbiting Carbon Observatory-2, 2023.  URL https://ocov2.jpl.nasa.gov/.
  * Gimeno García et al. [2011] S. Gimeno García, F. Schreier, G. Lichtenberg, and S. Slijkhuis.  Near infrared nadir retrieval of vertical column densities: methodology and application to SCIAMACHY.  _Atmospheric Measurement Techniques_ , 4(12):2633–2657, 2011.  doi:10.5194/amt-4-2633-2011.
  * Dennis et al. [1981] J. Dennis, D. Gay, and R. Welsch.  Algorithm 573: NL2SOL-an adaptive nonlinear least-squares algorithm.  _ACM Transactions on Mathematical Software (TOMS)_ , 7(3):369–383, 1981.  doi:10.1145/355958.355965.
  * Hochstaffl et al. [2018] P. Hochstaffl, F. Schreier, G. Lichtenberg, and S. Gimeno García.  Validation of Carbon Monoxide Total Column Retrievals from SCIAMACHY Observations with NDACC/TCCON Ground-Based Measurements.  _Remote Sensing_ , 10(2), 2018.  doi:10.3390/rs10020223.
  * Crisp et al. [2021] D. Crisp, R. Rosenberg, L. Chapsky, G. Keller Rodrigues, R. Lee, A. Merrelli, G. Osterman, F. Oyafuso, R. Pollock, G. Spiers, S. Yu, J. Zong, and A. Eldering.  Level 1B Algorithm Theoretical Basis: Orbiting Carbon Observatory – 2 & 3 (OCO-2 & OCO-3), 2021.
  * Gunson and Eldering [2019a] M. Gunson and A. Eldering.  OCO-2 Level 2 geolocated XCO2 retrieval results and algorithm diagnostic information, Retrospective Processing V10r, 2019a.
  * Crisp et al. [2020] D. Crisp, C. O’Dell, A. Eldering, B. Fisher, F. Oyafuso, V. Payne, B. Drouin, G. Toon, J. Laughner, P. Somkuti, G. McGarragh, A. Merelli, R. Nelson, M. Gunson, C. Frankenberg, G. Osterman, H. Boesch, L. Brown, C. Rebecca, M. Christi, B. Connor, J. McDuffie, C. Miller, V. Natraj, D. O’Brien, I. Polonski, M. Smyth, D. Thompson, and R. Granat.  Level 2 Full Physics Retrieval Algorithm Theoretical Basis: Orbiting Carbon Observatory – 2 & 3 (OCO-2 & OCO-3), 2020.
  * Schreier et al. [2019] F. Schreier, S. Gimeno García, P. Hochstaffl, and S. Städt.  Py4CAtS—PYthon for Computational ATmospheric Spectroscopy.  _Atmosphere_ , 10(5):262, 2019.  doi:10.3390/atmos10050262.
  * Gunson and Eldering [2019b] M. Gunson and A. Eldering.  OCO-2 Level 1B calibrated, geolocated science spectra, Retrospective Processing V10r, 2019b.