Optimal estimation
For configuration and interpretation of retrieval outputs, see Configuring an optimal-estimation retrieval.
The core expression of optimal estimation in ARTS follows from Rodgers [31]. He formulates the core expression of a measurement as
where \(\vec{y}\) is a measurement vector
(i.e., measurement_vec
in the ARTS workspace),
\(\vec{x}\) is the state of the model
(i.e., model_state_vec
in the ARTS workspace),
\(F\) is the model (i.e., ARTS itself), and
\(\epsilon\) is some measurement error that cannot
be modelled.
The forward model result from this is just
The goal of the optimal estimation is to find an \(\vec{x}\) that minimizes some cost-function, which in ARTS is defined as
where
where \(m\) is the number of measurements, \(\vec{x}_a\) is an a priori of \(\vec{x}\), \(\mathbf{S}_\epsilon\) is the covariance of the measurement noise, and \(\mathbf{S}_a\) is the covariance of the model state a priori.
Linearization
The iterative methods use a local linear approximation to the forward model around the current state \(\vec{x}_i\):
Here
is the Jacobian matrix
(i.e., measurement_jac
in the ARTS workspace).
One approach to minimize the cost function is to use a Gauss-Newton approach to update the state of the atmosphere. This might look like
The Jacobian and simulated measurement in this expression are evaluated at the current state. Iteration stops according to the convergence criterion and iteration limit. The state-step measures are defined below; settings are described in Configuring an optimal-estimation retrieval. A decrease in cost alone does not establish convergence.
Covariances and coordinates
A covariance matrix has variances on its diagonal. Its off-diagonal entries can be expressed in terms of standard deviations and a correlation matrix:
Thus each covariance entry has the product of the units of its two coordinates. The covariance is symmetric and positive semidefinite; the inverses in the optimal-estimation expressions require positive definiteness. A positive diagonal alone does not imply positive definiteness. The inverse covariance is the precision matrix. In general, correlations imply
For a coordinate transformation \(\vec{x}=f(\vec{t})\), the local covariance transformation is
This is exact for an affine transformation. A nonlinear transformation also changes the shape of the distribution; a Gaussian prior in native coordinates need not remain Gaussian in transformed coordinates.
For correlated measurement errors, the quadratic measurement cost can be written using a whitened residual. With a Cholesky factorization \(\mathbf{S}_\epsilon=\mathbf{L}\mathbf{L}^{\top}\),
Numerical state scaling is distinct from a change of statistical coordinates. For a diagonal matrix of positive scales \(\mathbf{T}\), solving \(\mathbf{H}\Delta\vec{x}=-\vec{g}\) is equivalent in exact arithmetic to
This changes the conditioning of the linear system while retaining the same objective and solution in the original coordinates.
For the measurement-space system, define
Then solve \(\mathbf{M}\vec{u}=\vec{r}\) and update \(\vec{x}_{i+1}=\vec{x}_a+\mathbf{S}_a\mathbf{J}^{\top}\vec{u}\). For a diagonal matrix of positive measurement scales \(\mathbf{D}\), the equivalent scaled solve is
Choosing \(D_{ii}=\sqrt{S_{\epsilon,ii}}\) scales by measurement noise standard deviations. This is not full whitening when measurement errors are correlated. This measurement scaling matrix is distinct from the prior-precision damping matrix introduced below.
Levenberg–Marquardt damping
Let \(J=m\chi^2\) denote the unnormalized objective. Define its half-gradient and the Gauss–Newton approximation to its half-Hessian by
The damped step used by ARTS is
Larger \(\gamma\) penalizes larger steps in the precision-scaled coordinates. At \(\gamma=0\), the step equals the Gauss–Newton step. The damping penalty modifies the local step, not the objective being minimized.
For a candidate step, the predicted decrease in the full, unnormalized objective is
The corresponding actual decrease and their ratio are
For an affine forward model, the local quadratic is exact, giving \(\rho=1\) whenever the predicted reduction is nonzero. The factor of two is required because \(\vec{g}\) and \(\mathbf{H}\) refer to half the objective. A ratio of two would mix the full cost with its half-cost prediction. A zero predicted reduction makes the ratio undefined; near zero, rounding of the computed costs can dominate it.
State-step convergence measures
For \(n\) retrieved state elements, the state-space formulation uses the Rodgers 5.31 measure
where the half-gradient \(\vec{g}\) is evaluated before the step. The measurement-space formulation uses the Rodgers 5.30 measure
Both quantities are dimensionless. The measures agree for an exact undamped Gauss–Newton step, but need not agree for damped or approximate steps. They are distinct from relative changes in the cost function and from an unweighted distance between state vectors.
The undamped step and its decrement are
For positive-definite \(\mathbf{H}\), \(\delta\geq0\), and \(\delta=0\) exactly when \(\vec{g}=0\). Its normalized value \(\delta/n\) equals both state-step measures for an exact undamped solve. The full objective’s predicted decrease for this step is \(\delta\). A small undamped decrement thus provides a stationarity measure when cost differences approach numerical resolution. By contrast, arbitrarily strong damping can make a damped step arbitrarily small while the gradient remains nonzero.
Gain, averaging kernel, and retrieval uncertainty
In the linear Gaussian model, the posterior covariance, gain, and averaging kernel are
The observation and smoothing contributions in state coordinates are
With the stated covariances and model assumptions, \(\widehat{\mathbf{S}}=\mathbf{S}_{\rm obs}+\mathbf{S}_{\rm smooth}\). For a nonlinear retrieval, evaluating these expressions at the retrieved state gives a local approximation. Their uncertainty interpretation requires that the prior covariance represents the assumed prior uncertainty. The observation contribution alone is not the full posterior covariance.
Information carried by the measurements
Factor the prior and measurement covariances as \(\mathbf{S}_a=\mathbf{L}_a\mathbf{L}_a^{\top}\) and \(\mathbf{S}_\epsilon=\mathbf{L}_\epsilon\mathbf{L}_\epsilon^{\top}\). The dimensionless sensitivity matrix and its singular value decomposition are
Each singular value \(s_i\) measures the response of a unit-prior state mode relative to the measurement error. Missing singular values when \(m<n\) are zero. In coordinates \(\vec{z}=\mathbf{L}_a^{-1}(\vec{x}-\vec{x}_a)\), the columns of \(\mathbf{V}\) describe independent prior modes. Their posterior variances are \(1/(1+s_i^2)\), giving the mode variance reductions
These are the eigenvalues of the averaging kernel expressed in prior coordinates. Their use as an information spectrum is described by Nesser et al. [25]. A mode with \(s_i=1\) loses half its prior variance; an unobserved mode with \(s_i=0\) retains its prior variance.
The physical state modes are the columns of \(\mathbf{P}=\mathbf{L}_a\mathbf{V}\). They satisfy \(\mathbf{P}^{\top}\mathbf{S}_a^{-1}\mathbf{P}=\mathbf{I}\). The posterior covariance can be reconstructed as
The marginal variance reduction for state element \(j\) is \(1-\widehat{S}_{jj}/S_{a,jj}\). This refers to an individual coordinate and is generally different from the variance reduction of a mode combining several state elements. A mode’s sign is arbitrary; equal singular values also allow rotations within their shared subspace.
Two scalar summaries are the degrees of freedom for signal and the Gaussian entropy reduction in bits:
These quantities depend on the Jacobian and assumed covariances, not on the realized measurement residual. For a nonlinear forward model they describe only the local linear approximation at the Jacobian’s state.
For a measurement and a forward prediction at the prior mean, define the innovation and its covariance by
Under the linear Gaussian model with independent prior and observation errors, \(\vec{r}^{\top}\mathbf{S}_r^{-1}\vec{r}\) follows a chi-squared distribution with \(m\) degrees of freedom and mean \(m\). This is a check of a prior prediction and the assumed uncertainties. Replacing the prior prediction with a fitted measurement does not give this distribution.
Transforming the Jacobian matrix
It is sometimes desired to transform the Jacobian matrix away from the native units that are available in the model. Instead of retrieving volume-mixing-ratio, it is sometimes more stable to retrieve the relative change in the volume-mixing-ratio, or perhaps even the logarithm of the volume-mixing-ratio, to avoid negative values creeping in to what is otherwise an intensive and numerical exercise.
There are three ways perform transformations on the Jacobian matrix in ARTS. You can change the measurement unit (e.g., from spectral radiance to brightness temperature), you can change the model state unit (e.g., from volume-mixing-ratio to relative humidity), or you can map the Jacobian matrix to another vector space (e.g., from Cartesian coordinates to spherical coordinates).
The first type of transformations are performed by an
SpectralRadianceTransformOperator
in ARTS. This is a local operation that is almost trivial to undo.
Generally, you should just use one of the
provided enumeration values of
SpectralRadianceUnitType to
set up these types of transformations.
The enumeration class also describes how the transformations
are done. This will not be repeated here.
The second and third types of transformations are a bit more complicated because they transform the model state vector. The optimal estimation methods are sometimes iterative, and they must be able to update the state of the model between iterations. These two types of transformations are therefore required to be able to transform the model state vector from native to non-native units, and vice versa, as well as to transform the Jacobian matrix from native to non-native units (but only in one-way, this operation does not need to be reversible).
From a mathematical point of view, these two behave very similar. But from a practical point of view they are quite different. Transforming the model state vector requires just the model state vector of the native units and the Jacobian matrix that corresponds to just that model state parameter. Mapping the model state vector to non-native units may require knowing more than just a single model state parameter.
Core mapping/transformation expression
If we define the native units of \(\vec{x}\) as \(\vec{t}\) so that
there must be a reversible functions so that
for any transformation or mapping to work. It must also be possible to take the partial derivative of \(\vec{t}\) with regards to \(\vec{x}\).
The Jacobian \(\mathbf{J}\) is the derivative with respect to the retrieved coordinates \(\vec{x}\). The forward model first computes derivatives in its native coordinates \(\vec{t}\). If we introduce
it is clear we can write
This step right here is what we consider the transformation of the Jacobian matrix. To make use of this style of transformation, we must provide matching \(f\) and \(f^{-1}\), as well as a way to compute the partial derivative of \(f^{-1}\) with regards to \(\vec{x}\).
See Configuring an optimal-estimation retrieval for assigning these operators to a Jacobian target.
Relative retrievals
This is a model state vector transformation. By relative retrievals, we mean that the value itself is not retrieved, but instead its ratio is retrieved.
In this scenario:
where \(\oslash\) and \(\odot\) are element-wise division and multiplication, respectively. \(\vec{t}_0\) is simply the a priori value of \(\vec{t}\).
Note
The first iteration of a retrieval setup is going to be \(\vec{x} = \vec{1}\).
Logarithmic retrievals
This is a model state vector transformation. By logarithmic retrievals, we mean that the value itself is not retrieved, but instead its logarithm is retrieved.
In this scenario:
where the exponential and logarithmic operations are element-wise.
Logarithmic relative retrievals
This is a model state vector transformation. By logarithmic relative retrievals, we mean that the value itself is not retrieved, but instead the logarithm of its relative value is retrieved.
In this scenario:
where the operations are still element-wise on the product that is created.
Note
The first iteration of a retrieval setup is going to have \(\vec{x} = \vec{0}\).
Relative humidity retrievals
This is a model state vector transformation. By relative humidity retrievals, we mean that the value itself is not retrieved, but instead the its conversion to relative humidity is retrieved.
In this scenario:
where \(\vec{p}\) is the pressure at the position of \(\vec{t}\), \(\vec{T}\) is the temperature at the position of \(\vec{t}\), and \(p_{\textrm{sat}}\) is a user-provided method to compute the element-wise saturation pressure.
Tip
There is a flag that can be provided to this transformation that turns negative relative humidities off.
Note
Be aware that the implementation in ARTS is general, and that while you can choose to treat temperature as, e.g., relative humidity… please don’t. It makes sense only for some species.
Absolute field retrievals
This is a model state vector mapping. By absolute field retrievals, we mean that the value itself is not retrieved, but instead the absolute value of the field is retrieved.
In this scenario:
where the subscripts \(u\), \(v\), and \(w\) indicate the three components of the vector north, east, and up, respectively, and \(\mathbf{J}_u'\), \(\mathbf{J}_v'\), and \(\mathbf{J}_w'\) are the Jacobian matrices of these three components in the native units of the model state vector.
Note
Neither \(\vec{\theta}\) nor \(\vec{\phi}\) are part of the model state vector, but are instead derived from the model state vector and are used to map the Cartesian coordinates to spherical coordinates. They are fixed during the retrieval process and are not updated between iterations.
Matching-grid cross correlation
Let the marginal prior covariances of temperature and log-water be \(\mathbf{S}_T=\operatorname{diag}(\sigma_{T,i}^2)\) and \(\mathbf{S}_q=\operatorname{diag}(\sigma_{q,i}^2)\), where \(q_i=\ln(\mathrm{VMR}_i)\). Correlating only matching grid points with coefficient \(\rho\) gives
For positive marginal variances this two-target covariance is positive definite exactly when \(|\rho|<1\). The conditional prior mean is
and its conditional variance is \(\sigma_{q,i}^2(1-\rho^2)\). Correlation therefore expresses a statistical preference with residual uncertainty, while preserving the unconditional marginal variances. With more than two correlated targets, all pair coefficients must jointly define a positive-definite covariance.
Reduced optimal estimation
oemCalcReduced uses a fixed state expansion and measurement projection:
For \(n\) full states and \(m\) full measurements, \(\mathbf{B}\in\mathbb{R}^{n\times r}\) has independent columns and \(\mathbf{C}\in\mathbb{R}^{q\times m}\) has independent rows. Substitution in the prior penalty and propagation of measurement errors give
The reduced objective is
This restricts the state to an affine subspace; it does not marginalize discarded state variability into the measurement error. The reduced gain and its mapping back to the original coordinates are
LM uses the projected damping matrix
Thus invertible, square reductions preserve the objective and the exact LM linear step under coordinate transformation. Numerical tolerances can still affect the iteration path; the convergence measure uses the reduced state dimension.
For the whitened singular-value decomposition \(\mathbf{L}_\epsilon^{-1}\mathbf{J}\mathbf{L}_a =\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\top}\), first consider the full orthogonal matrices \(\mathbf{U}\in\mathbb{R}^{m\times m}\) and \(\mathbf{V}\in\mathbb{R}^{n\times n}\). The complete bases
are invertible coordinate transformations. No information is lost, including in the null spaces:
Here \(\boldsymbol{\Sigma}\in\mathbb{R}^{m\times n}\) has \(\min(m,n)\) singular values, ordered from largest to smallest, on its diagonal and zeros elsewhere. Additional directions of the larger basis carry zero local information. The bases and singular values together preserve the full local Gaussian problem.
Reduction is a separate choice: take the first \(r\) columns of \(\mathbf{B}_{\rm full}\) and first \(q=\min(r,m)\) rows of \(\mathbf{C}_{\rm full}\):
Both reduced covariances are then identity matrices. The discarded local degrees of freedom and information are
An information-budget reduction chooses the smallest retained rank satisfying the prescribed upper bounds on \(\Delta d_s\) and/or \(\Delta H\). The bounds apply to the sums, so many individually weak modes can collectively require retention. Zero-information modes have \(s_i=0\); their location is determined by the singular directions of the whitened Jacobian, not by individual zero entries in the physical matrices. A zero budget removes only zero contributions in the computed spectrum. Finite-precision residuals in mathematically null modes can require a positive budget.
Retaining every nonzero singular mode preserves the posterior of the linear Gaussian problem. Dropping nonzero modes is approximate. Measurement compression preserves the likelihood’s state dependence on the retained subspace if it retains the column space of \(\mathbf{L}_\epsilon^{-1}\mathbf{J}\mathbf{B}\). Discarded orthogonal noise residuals contribute a constant to the full objective. In general, selecting physical channels does not satisfy this condition, particularly when their errors are correlated.
For a single retrieved parameter with Jacobian column \(\vec{j}\), the one informative measurement combination has weights
The overall sign is arbitrary. This retains the local measurement Fisher information \(\vec{j}^{\top}\mathbf{S}_\epsilon^{-1}\vec{j}\).
For independent channel noise with variances \(\sigma_i^2\), the weights are proportional to \(j_i/\sigma_i^2\). A flat part of a spectrum carries little wind information when wind primarily shifts a line: if \(F(\nu,v)\approx F_0(\nu-\delta\nu(v))\), then
Opposite line slopes therefore enter with opposite signs. The exact center of a symmetric line has zero first-order shift sensitivity; a weighted line contrast is distinct from retaining just that center channel. Saturation and pressure broadening can also leave informative slopes outside the narrow central feature. Far-wing weights become small according to the Jacobian and noise, not according to a prescribed frequency window.
Another measurement compression preserves explicit groups of proportional Jacobian rows. Write \(\mathbf J=\mathbf T\mathbf R\), where row \(i\) of \(\mathbf T\) has one nonzero amplitude \(a_i\) in the column corresponding to its group. The rows of \(\mathbf R\) are the distinct sensitivity directions. A sufficient measurement projection for this linear Gaussian model is
The grouped directions have disjoint channel support, so \(\mathbf T\) has independent columns and \(\mathbf H\) is positive definite. Both the likelihood precision and its state-dependent linear term are preserved. With diagonal measurement noise, each row of \(\mathbf C\) can instead be normalized to unit noise variance:
Here \(\mathbf C\) has one nonzero per channel. Correlated noise can require entries outside a channel’s group. Grouping proportional rows may retain more measurements than an SVD: it does not combine general linear dependencies between different row directions. A zero-sensitivity group may also be retained without adding state information. For nonlinear forward models, these statements apply at the chosen linearization.
For retained singular modes, the truncated local posterior must preserve discarded prior uncertainty:
Using only the projected reduced covariance assigns zero uncertainty to discarded modes. This expression describes the truncated local model; a null mode of one nonlinear Jacobian need not remain null elsewhere. Neither the loss estimates nor this covariance bound the error of a particular nonlinear retrieval.