Developing the OEM interface

The user guide is in Configuring an optimal-estimation retrieval and the mathematical formulation is in Optimal estimation. This page separates the initial interface cleanup from further work on numerical algorithms and public APIs. The statistical objective and existing method names should remain stable unless a change is explicitly documented and tested.

Responsibilities and regression baseline

src/m_oem.cc validates inputs and connects workspace inputs and outputs to the inverse problem. src/oem.h contains the agenda adapter, solver wrappers, and logging. 3rdparty/invlib/src/invlib/map.cpp drives the optimization and computes costs and convergence measures. The LM step and linear solvers live below that layer. Covariance storage and inverse operations are implemented in src/core/jacobian/covariance_matrix.cc; retrieval covariance setup also involves src/m_covmat.cc.

The first cleanup centralizes method selection and LM configuration while retaining the public defaults. Direct and CG LM methods now use all six settings, the supplied stop_dx, and the same diagonal of prior precision for damping. Direct measurement-space li_m and gn_m reuse MFORM with DirectMeasurementSolver. State normalization is rejected for measurement-space methods instead of applying state scales to a measurement-space system.

Measurement-space convergence uses the state-space Hessian metric (Rodgers 5.30), so it remains well-defined when the numbers of measurements and states differ. Diagnostics come from the formulation that actually ran. LM history is collected independently of terminal output, and aliases have the same initial-cost behavior. Final nonlinear Jacobians are refreshed at the returned state when needed for the requested gain. Outputs that would otherwise retain an earlier run’s gain or errors are cleared, and the clear_matrices policy also applies to a skipped retrieval.

src/tests/test_oem_methods.cc exercises the workspace interface using small deterministic forward models. CTest registers the individual cases as cpp.fast.oem.<method>, plus settings and validation cases. Each supported spelling, including ml and ml_cg, is part of the contract. These tests run without the optional MATLAB integration required by the legacy invlib-level tests in src/tests/test_oem.cc. They complement the atmospheric retrieval examples under tests/core/vmr. Exact initial solutions and underdetermined problems are covered; a zero right-hand side now returns a zero CG solution directly.

The numerical reference for an affine model should come from an independently formed linear posterior solution. Check the state, fitted measurement, cost decomposition, gain, and output dimensions; agreement between two methods alone can preserve a shared bug. Include correlated covariances, unequal state and measurement dimensions, a non-prior starting state, and equivalent scaled and unscaled systems. Nonlinear cases need to test accepted-state Jacobians and meaningful damping changes, rather than only checking that a residual decreased.

Forward-model reuse and agenda ownership

AgendaWrapper tracks the exact state of its most recent successful simulation separately from the state of its stored Jacobian. It borrows the workspace fit and Jacobian storage; only the state tags need additional storage proportional to the state size. Matching an earlier state approximately is not sufficient for reuse. Initial evaluations, accepted LM trial values, and final Jacobians must not be repeated when the required result is already current. A value-only trial preserves the last Jacobian, but changes the physical model fields. Returning to that Jacobian’s state therefore still requires a value-only evaluation to restore the fields. Invalidate cache entries before executing an agenda and publish new state tags only after successful execution and dimension checks.

The criterion_needs_measurement trait identifies the built-in Rodgers530 and Rodgers531 types. Their state-step tests allow a continuing Gauss–Newton iteration to obtain its fit and Jacobian together, without a preceding value-only evaluation. Unknown convergence criteria, including overrides derived from these classes, default to requiring the new measurement. Retain that conservative behavior when extending the trait. At the final iteration, obtain the fit for costs and diagnostics, then refresh a nonlinear Jacobian only if it is needed for retained matrix outputs. clear_matrices=1 skips that final derivative refresh as well as gain construction. Preserve evaluation-count regressions alongside numerical oracles, including rejected trials, zero steps, iteration limits, failed evaluations, and both matrix-output policies.

The generated agenda executor borrows explicit inputs and outputs through non-owning Wsv wrappers. Ordinary workspace sharing copies handles, not field or matrix payloads. Agenda::finalize identifies internally modified inputs that need isolation; copy_workspace copies those values, and Agenda.document() lists them. Nested agendas have their own copy lists. Literal values stored in agenda methods are also copied when executed. Do not retain an executed local workspace blindly: copy_only_workspace can copy its already modified scratch values, changing initialization between trials. Any future scratch reuse must preserve fresh-call initialization and failure behavior.

Python tuple-style operators have additional argument and return conversion costs. Their generated C++ adapters move the owning return tuple’s elements into outputs; never move from Python-owned objects or borrowed inputs. CallbackOperator instead supplies a restricted workspace sharing its declared variables, allowing outputs to be updated in place. Keep the native allocation-transfer regression in test_agenda_operator.cc and the Python retained-object regression in tests/core/agenda/operator_return_values.py when changing these bridges.

State mapping and derivative selection

AgendaWrapper always passes the complete target set as model_state_targets. It selects a const reference to either that same set or a local empty set for jac_targets before entering the agenda. Both objects outlive the synchronous agenda execution. Do not copy the populated targets or clear an agenda input to implement value-only calls.

Methods mapping trial states into physical fields use model_state_targets. Radiative transfer and Jacobian transformations use jac_targets. Measurement-error values and derivatives likewise use the full mapping and derivative targets respectively. This distinction preserves trial-state updates and error values when derivatives are disabled. The predefined agendas must not list either target set among their copied inputs.

Bounded inner iterations

An outer max_iter does not bound the work of an inner linear solve or LM trial search. Both enforce independent limits. CG returns its current iterate on exhaustion and invokes an optional warning callback. OEM installs this callback and records one warning per OEM call in oem.diagnostics.errors; the outer optimizer continues and retains its own diagnostic status. LM trial exhaustion still raises an error. OEM maps errors caught during inversion to status Error and records the explanation in oem.diagnostics.errors. Exhausting the existing LM damping range retains its status DampingLimit behavior. Gauss–Newton must propagate linear-solver exceptions with their nested cause. Returning an empty step after catching an error hides the failure and can cause invalid vector operations in the measurement-space formulation.

The native ConjugateGradient constructor and both preconditioned variants accept max_iterations after verbosity. Zero selects \(\max(1000,2d)\) for a system of dimension \(d\); a positive value sets an explicit bound. Custom convergence predicates cannot disable it. cg_tolerance and cg_max_iter are validated at the OEM boundary. Static native solver preconditions remain assertions. Nonpositive curvature or a nonpositive preconditioned residual inner product is a runtime numerical breakdown, reported through the OEM error diagnostics.

CG records its stop reason on every solve, including a zero right-hand side. GN refuses a truncated step, preserving the last state and reporting LinearSolverLimit. LM rejects it and retries with increased damping; if the final attempt is still truncated, it reports the same status. Do not let the measurement-space mapping apply its prior offset when a step was rejected. Regressions cover all formulations and final-state Jacobian reuse.

The native LM optimizer provides get_maximum_trials() and set_maximum_trials(); the positive trial limit defaults to 100 per outer iteration. A local counter bounds rejected-step retries, and a damping update must make progress before another trial is attempted. This catches multiplication that rounds back to the current damping. The trial budget counts linear solves, including any additional undamped solve used to check stationarity. The LM trial limit is exposed as settings.lm.maximum_trials (default 100). Preserve the termination regressions when changing convergence predicates, trial acceptance, or damping updates.

LM acceptance and stop outcomes

LMStopReason records termination independently of the damping value. None means the optimizer can continue; Stationary is successful termination. DampingLimit returns a zero step and maps to workspace status DampingLimit. TrialLimit, DampingStalled, LinearSolverFailure, and NumericalFailure accompany exceptions, which OEM maps to status Error when caught during inversion. Ordinary convergence after an accepted step still uses the configured convergence criterion and damping gate. Both that criterion and Stationary map to status Converged; outer iteration exhaustion retains status IterationLimit.

All MAP formulations consult the optimizer’s explicit outcome before testing the returned step for convergence. A zero step returned after damping exhaustion must not count as convergence. Rejected trials must never update the accepted state, and damping history must contain physical damping values, not a numerical sentinel such as maximum + 1. At a large finite maximum, that addition can round back to the maximum itself. The ordinary convergence gate returns a zero tolerance while damping exceeds its limit, so the strict comparison also rejects a rounded-zero state change.

The generic helper in optimization/minimize.h also honors optional stop_iteration() and converged() hooks. A reported optimizer failure returns before testing the criterion or updating the state. A stationary optimizer’s verified step is applied, but success still requires J.criterion to meet the helper’s own tolerance. Iteration exhaustion returns one. Retain native regressions for these outcomes and for custom minimizers without either hook.

MAPBase::model_cost_scale() returns two: MAP reports the full quadratic cost, while its normal equations use the half-gradient and half-Hessian. The LM prediction includes this scale; an exact affine model therefore has an actual-to-predicted reduction ratio of one. Generic objectives without this optional method retain scale one and must provide mutually consistent cost, gradient, and Hessian definitions. Keep regressions for both conventions when changing the reduction formula or acceptance thresholds.

LM compares reductions against 32 * epsilon * max(abs(old_cost), abs(new_cost)) before dividing them. If actual or predicted reduction is unresolved at that scale, it checks the undamped normal equations once within the current step. Numerical stationarity requires a finite, nonnegative undamped decrement below n * stop_dx, a finite, nonnegative predicted reduction within the current cost’s roundoff scale, and a candidate cost indistinguishable from the current cost at that scale. This step can converge independently of the damping gate: its size is checked without damping. An exactly zero gradient is also stationary. These checks must not substitute the damped step for the undamped decrement; very strong damping can hide a large remaining error. Resolved acceptance requires finite, positive actual and predicted reductions and a true reduction ratio of at least 0.5; damping is decreased when the ratio exceeds 0.75. These thresholds use consistently scaled costs, not the inflated ratio in the older implementation.

Keep workspace regressions for all four LM spellings at damping 1e20, tight affine tolerances, and an exactly stationary state with nonzero cost. Native tests also cover explicit stop reasons, objective scaling, rejected steps, and bounded or stalled damping updates. Check the state, fitted measurement, costs, gain, and damping history as well as the final status.

Remaining numerical work

The following items require separate implementation and regression work.

1. Expose inner-solver controls and diagnostics. Report linear iterations and residuals. Relative CG tolerance and the linear iteration budget are already exposed. Keep these separate from outer max_iter and the six damping controls. Extend numerical coverage to nearly zero right-hand sides and ill-conditioned positive-definite systems; bounded termination alone does not establish the accuracy of a difficult solve.

2. Extend public stop diagnostics. The workspace retains its numeric status mapping; native LM stop reasons are not a separate workspace output. A future result could expose the precise stop reason and distinguish successful completion of a one-step linear solve from exhaustion of an iterative convergence budget.

3. Replace dense normalization and explicit gain inversion. Store state scales as a vector and apply diagonal products without allocating an \(n\times n\) dense scaling matrix. Compute the gain by solving the posterior precision system for its right-hand sides. Measure peak memory as well as runtime, and retain the option to skip gain production. Consider factorization reuse and a preconditioner appropriate to each formulation only after correctness is established.

4. Make agenda state and failure behavior explicit. Validation, covariance inversion, and the initial agenda evaluation can still throw directly; status Error covers errors caught during inversion. Trial forward model evaluations mutate atmosphere, sensor, and surface workspace data. Define which state and diagnostics remain valid after an exception or a rejected trial. A future result should distinguish the last accepted state from a failed trial without treating NaNs as the only status record.

Covariance validation and construction

CovarianceMatrix.validate(expected_size=-1, relative_tolerance=1e-10, max_dense_elements=10_000_000) checks an existing physical covariance without modifying it. Validation covers block dimensions and ranges, coverage, unique upper-triangular block identifiers, non-null matrices, finite values, positive variances, symmetry, and positive definiteness. Symmetry and supplied-inverse comparisons use variance-scaled coordinates so mixed units do not let a large-variance quantity hide an invalid small-variance block. Positive definiteness is checked over connected components of the complete covariance, including cross-block correlations.

Each represented inverse component must be complete and consistent with its connected covariance component. Whole independent components may remain uncached during validation. In normalized coordinates, the matrix-product residual is compared with relative_tolerance times the component size; the diagonal fast path uses relative_tolerance directly. compute_inverse uses the same checks and inverts uncached independent components. This catches invalid physical covariances on the existing OEM inversion path without imposing the standalone analysis’s memory limit on existing retrievals. It preserves the inverse-only representation used internally for LM damping; that representation is not accepted as a physical covariance by validate or the standalone information report.

Replacing blocks or accessing mutable blocks invalidates cached inverses. Adding a covariance block invalidates the affected connected components while retaining independent cached components. Shared matrix aliases can still outlive those access points; validation before inversion must continue to detect a stale supplied inverse. Do not silently overwrite asymmetric covariance entries during inversion. The covariance-addition helper in src/m_covmat.cc now inserts the optional supplied inverse as an inverse and validates its inputs before mutation.

Sparse diagonal validation and inversion retain diagonal storage. Explicit validation of non-diagonal connected components requires dense work; check max_dense_elements before allocating each component matrix. This guard bounds individual dense arrays, not total peak memory. Keep regressions for diagonal and correlated blocks, mixed scales, incomplete coverage, inverse-only inputs, stale caches, and supplied inverse consistency. Validation must not silently repair a statistical model.

Useful construction helpers would accept standard deviations and an optional correlation matrix or named correlation kernel, then construct covariance blocks. Parameter names should state whether they accept a standard deviation, variance, covariance, or precision. A report could identify each block by retrieval target and include its units, variance range, and factorization or conditioning diagnostics. Covariance repair, such as adding a diagonal term or changing correlations, must be an explicit user choice because it changes the inference.

Standalone information analysis

python/src/pyarts3/retrieval.py provides information and information_from_workspace independently of the OEM optimizer. The workspace adapter reads the current Jacobian and the two physical covariances. Neither entry point executes a forward agenda, runs OEM, updates a workspace, or selects new covariance values. This separation lets a user compare proposed measurement and uncertainty models with a fixed Jacobian before a retrieval, and examine local sensitivities after one. The mathematical definitions belong in Information carried by the measurements; usage and interpretation belong in Correlating temperature and log-water on the same grid.

The report uses covariance factors to obtain a dimensionless Jacobian, then an exact singular value decomposition. singular_values and mode_variance_reduction include all state directions, padding the unobserved modes with zero singular values when measurements are fewer than states. state_modes contains physical state directions scaled to unit prior uncertainty; measurement_modes refers to whitened measurement coordinates. Retain the distinction between mode and marginal variance reduction. Mode signs and bases inside a degenerate subspace are not stable identifiers for comparisons between runs.

Use factor solves rather than forming covariance inverses for this analysis. Diagonal measurement covariances must not allocate a full measurement-square matrix. Guard dense work with max_dense_elements; even a diagonal prior still needs a state-square basis when all modes are returned. The analysis guards an estimate of 3*m*n + 4*n*n elements as well as the cumulative dense factor sizes; LAPACK can need additional work arrays. Future truncated or operator-based analyses must state which parts of the spectrum and uncertainty they approximate, and must not label retained-mode information as the complete information content.

information requires both an explicit measurement and prior_prediction for innovation_chi_square. information_from_workspace reads the workspace measurement only when given prior_prediction; it does not accept a measurement override. Do not infer the prior prediction from a workspace’s fitted measurement: its generating state is not known to this API. The chi-squared reference distribution assumes a linear model, a prediction at the prior mean, and the stated independent Gaussian error sources. The spectrum itself is independent of the measured residual.

Regression oracles should include diagonal analytic systems, correlated covariances, unequal dimensions, exact null directions, consistent coordinate changes, and posterior covariance reconstructed from the returned modes. Cross-check the innovation statistic against an independently formed innovation covariance. Keep workspace immutability and dense-memory guards covered. Plot uncertainty ratios when combining state elements with different units; retain absolute standard deviations and labels in the numerical report.

Helping users choose settings

Named LM settings are implemented by LevenbergMarquardtSettings. Its keyword-only constructor and mutable fields expose initial_damping, decrease_factor, increase_factor, maximum_damping, damping_threshold, and convergence_damping_limit. The native settings representation provides shared validation for the Python object and the optimizers used by oemCalc. Both direct and CG optimizers must continue to use the same named controls. The type lives in src/core/jacobian/oem_settings.h and src/core/jacobian/oem_settings.cc without an invlib dependency; src/python_interface/py_retrieval.cpp binds the Python interface.

OptimalEstimationSettings owns the method enum, outer iteration controls, CG limits, nested LM controls and output flags. Both workspace calculations validate it once before covariance preparation, then pass it by const reference to the shared iteration dispatcher. Keep normalization and bases in OptimalEstimationData because their dimensions belong to the problem. The settings type uses aggregate XML serialization; extend the Python constructor, properties, pickle state and round-trip regressions together when adding fields. The six-value LM shorthand initializes damping only; the separate maximum_trials field keeps its default.

The Python LM constructor validates all seven LM settings, and each field setter validates a temporary copy before replacing the stored object. An invalid edit must preserve the previous values and report the affected setting. This also avoids losing useful validation errors through nanobind’s implicit-conversion error handling. Related field changes must either keep each intermediate configuration valid or use a replacement object constructed with the desired keyword arguments together.

validate() checks a configured object at the OEM boundary. describe() explains its fields. Both settings types have XML storage; OEM uses the named defaults when omitted. The outer settings are a workspace group; the nested LM type is registered as a group friend. The C++ API remains named. The Python binding accepts std::array<Numeric, 6> through its implicit constructor, which validates the six-field input. There is no as_vector() output conversion. Do not describe the threshold as a hard minimum: it controls both restart after rejection and switching a proposed decrease to zero. The convergence damping limit gates the existing state-step criterion using the updated damping. Preserve the accepted/rejected-step regressions when changing either behavior. A preset claiming particular convergence or performance properties needs documented assumptions and representative retrieval benchmarks.

src/tests/test_oem_methods.cc tests the native settings and solver behavior. tests/core/jac/oem_lm_settings.py covers the Python interface through the actual workspace call. Keep named settings covered for all four LM spellings, including correlated priors and nonlinear rejected trials; comparing complete results and damping histories catches changes that an endpoint-only test would miss.

Further API work should keep existing scripts usable:

• Represent iteration family, linear solver, and formulation separately internally. A single method descriptor table can drive parsing, supported-name documentation, and test enumeration. Existing string spellings and aliases can remain as the public compatibility layer.

• Provide a setup report with state and measurement sizes, estimated matrix storage, chosen coordinates and scales, covariance diagnostics, and solver settings. Recommend a method with a reason while leaving the explicit choice available. Dimensions can suggest a formulation; they cannot establish linearity or predict nonlinear convergence.

• Add an optional local Jacobian check using representative perturbations in prior-scaled coordinates, and report discrepancies by measurement and retrieval target. Start with a deterministic small-model example before applying it to expensive radiative-transfer agendas.

• Extend OptimalEstimationDiagnostics with accepted iteration costs, state-step measures, damping, trial rejections, linear residuals, and forward-model call counts. There is no legacy positional adapter. A user can then distinguish poor model linearity, a difficult linear solve, an exhausted iteration budget, and a mismatch in assumed errors.

Method recommendations should be assessed on linear, weakly nonlinear, strongly nonlinear, correlated, poorly scaled, and large sparse cases. Compare final objective and state as well as runtime, peak memory, and forward-model calls. Automatic changes to prior or measurement covariance require scientific assumptions that cannot be inferred from solver progress alone. Keep those assumptions visible in both the setup and uncertainty reports.

Matching-grid covariance helper

oemStateCovmatCorrelateConstant resolves finalized atmospheric target keys and checks actual coordinate grids, not only vector lengths. It assumes pointwise retrieval coordinates; custom mappings that mix grid points are outside its contract. It constructs an upper-triangular sparse cross block from existing marginal standard deviations. It validates a candidate containing all existing pairs before assignment, preserving the original on failure and discarding inverse caches on success. Extensions must keep these guarantees and must not silently rescale transformed covariances.

The correlation method accepts two independent shared-pointer input variants in one implementation. Each argument is converted separately to the common atmospheric key type. The regression exercises all atmospheric-key/species combinations and the generated workspace dispatch through an agenda.

Measurement-space noise scaling

When measurement_vec_normalization is nonempty, the measurement-space direct and CG paths solve

\[(\mathbf{D}^{-1}\mathbf{M}\mathbf{D}^{-1})\vec{v} = \mathbf{D}^{-1}\vec{r}, \qquad \mathbf{M}=\mathbf{J}\mathbf{S}_a\mathbf{J}^{\top}+\mathbf{S}_\epsilon, \qquad \vec{u}=\mathbf{D}^{-1}\vec{v}.\]

Here \(\mathbf{J}\) is the measurement Jacobian and \(D_{ii}\) are the supplied scales. Empty scales disable the transformation. The helper computes the suggested noise scales \(\sqrt{S_{\epsilon,ii}}\); OEM does not select them automatically. NoiseScaledSystem applies this operation without materializing \(\mathbf{M}\) or \(\mathbf{D}\). The existing state-sized normalization remains exclusive to state-space solvers. The relative CG tolerance is measured in the scaled system. oemMeasurementCovmatNormalization exposes the same standard deviations as a generic Vector output. Correlated noise is not fully whitened. Regression tests compare with the affine analytic solution and check state and cost invariance under independent measurement-unit changes.

Direct measurement-space assembly

DirectMeasurementSolver opts into dense_measurement_system. GaussNewton forwards this compile-time policy to MFORM; other optimizers/solvers default to lazy evaluation. MFORM materializes \(\mathbf{J}^{\top}\) and computes

\[\mathbf{B}=\mathbf{S}_a\mathbf{J}^{\top}, \qquad \mathbf{H}=\mathbf{J}\mathbf{B}+\mathbf{S}_\epsilon.\]

It reuses \(\mathbf{B}\) for the state update. Here \(\mathbf{H}\) is the measurement-space system, not the state-space half-Hessian. The direct solver copies \(\mathbf{H}\) for optional scaling and factorization. This uses additional \(n\times m\) temporary storage but avoids repeated covariance applications. CG remains lazy. The ARTS adapters expose transpose_view() as a non-owning const strided view. Dense MFORM passes this view directly to covariance multiplication, avoiding an additional \(n\times m\) Jacobian copy. The source Jacobian must remain alive and its storage stable until multiplication finishes. Owning transpose() remains available for invlib expression materialization; other backends retain that fallback. Gain postprocessing remains state-space and must be included in performance comparisons. New methods share the affine, nonlinear, underdetermined, normalization and failure regression fixtures.

Fused covariance addition and GEMM

The ARTS reference-matrix adapter provides multiply_add(B, Se) for the dense MFORM path. It initializes the result to zero, adds the covariance blocks directly (including off-diagonal blocks), and calls mult with \(\alpha=1\) and \(\beta=1\). Thus \(\mathbf{H}=\mathbf{J}\mathbf{B}+\mathbf{S}_\epsilon\) needs one GEMM and no separate dense covariance temporary or post-GEMM addition pass. Other invlib backends keep the multiplication-plus-addition fallback. The transpose-view path removes the explicit Jacobian transpose copy. The direct solver’s scaling copy remains.

Diagonal covariance regressions

Every method-specific OEM CTest also runs test_diagonal_covariances. It uses \(\mathbf{S}_a=\operatorname{diag}(4,2)\) and \(\mathbf{S}_\epsilon=\operatorname{diag}(1,2,1/2)\) with the affine fixture, testing all four combinations of dense/sparse prior and noise storage, with normalization enabled and disabled. The reference state, gain, total cost and measurement cost are rational values derived independently from the two-state normal equations. Linear methods must take exactly one step. These checks run as ordinary CTest regressions.

Exact diagonal covariance operations

Covariance matrix multiplication now detects exact diagonal structure across independent blocks and scales the destination directly. mult_inv and vector solve divide by covariance diagonal entries without requiring an inverse cache. Sparse diagonals are inspected in \(\mathcal{O}(\mathrm{nnz})\); dense blocks require an exact off-diagonal scan. Structure is not cached because callers can hold mutable/shared block storage. Detection uses no threshold and never discards small correlations. These paths use \(\mathcal{O}(n)\) scratch space rather than a full matrix-result temporary.

Component solves now use an internal variant of diagonal values and dense Cholesky factors, grouped by exact block connectivity. Ordinary covariance mult_inv and solve calls prepare/reuse this representation without forming inverse blocks. Matrix and Sparse remain the public input types. Adding an off-diagonal block merges components; replacing/removing blocks can split them. Coupled sparse components currently use dense Cholesky, matching the former dense inverse component storage; sparse Cholesky and banded factorization are not implemented.

OEM calls CovarianceMatrix::prepared once for each covariance before iteration. Preparation checks exact values and block layout, validates changed inputs, and publishes a detached, read-only snapshot. Its block storage belongs to the snapshot, so subsequent edits through aliases of the source do not change a running retrieval. Unchanged inputs reuse the prepared snapshot. Mutable covariance copies have independent preparation holders. A completed compute_inverse records the exact validated storage; preparation reuses that validation only if all covariance and inverse values, shapes and block metadata still match. This avoids duplicate validation in cold calls with explicit inverse preparation.

Preparation publication is protected by a mutex. Published snapshots support concurrent solves; callers must synchronize source mutation with preparation and must not mutate storage obtained from a prepared snapshot. Iterations do not compare snapshots, validate covariances, or build inverse caches. Standalone operations on mutable covariances retain their defensive checks. Concurrent solves, inverse construction, explicit validation and preparation synchronize cache publication, retaining factors while arithmetic runs outside the lock. External block edits still require synchronization. Prefer prepared snapshots for repeated concurrent solves to avoid cache checks and locking. Moved-from matrices remain valid empty objects and can be reused. New or changed supplied inverses still receive the full deterministic consistency check; only identical previously validated contents skip that cubic work.

prepared(true) requests explicit precision for consumers that need it. OEM requests this for state-space methods and gain output. Measurement error covariance stays diagonal or factorized unless inverse blocks were supplied. Any supplied inverse currently selects inverse application for the entire covariance, with missing component inverses materialized; mixed inverse/factor preparation within one covariance is not implemented. Calling an explicit precision consumer on prepared(false) without inverse blocks throws; request the required representation before entering the loop. CG versus direct method selection remains the user’s choice.

Preparation adds detached storage and a linear comparison once per OEM call. It is not a zero-copy interface. This cost avoids repeated validation and alias checks inside iterative solvers while preserving detection of source changes between retrievals.

Regression tests cover mixed Matrix/Sparse inputs, joining/splitting components, retained-reference mutation after factorization, copies, multiple RHS and left/right solves. Existing supplied-inverse and all-method OEM regressions remain applicable.

Diagnostics interface

OptimalEstimationDiagnostics owns the LM history and error/warning list. Reset the complete object for each OEM call. iterations is an Index starting at zero; only unavailable costs use NaN. Translate invlib’s return code at the OEM boundary to OptimalEstimationStatus, defined in src/core/options/arts_options.cc. Its formatting, XML, and Python support come from the standard options machinery. Both diagnostics and LM settings use their aggregate XML representation; status is serialized by name, independently of enum ordinals.

ReducedOEM adapter

oemBasisCalc prepares model_state_basis_mat, measurement_basis_mat and oem_basis_singular_values from the full Jacobian and covariances. By default both bases are square; neither null space is discarded. full_matrices=0 instead retains \(\min(m,n)\) modes on each side, omitting only the extra null-space vectors of the larger side. The spectrum contains \(\min(m,n)\) values, including zeros. Additional directions of the larger square basis have zero information. Keep these three outputs matched, including their mode ordering. Mixing equally sized decompositions cannot be detected by dimension checks.

oemBasisReduce reads the basis matrices and full spectrum and truncates the matrices in place. It also sets oem_basis_lost_dofs and oem_basis_lost_information_bits relative to the original decomposition. The spectrum remains intact. Selection may be repeated to remove more modes, but cannot restore removed modes; requests requiring unavailable directions fail before modifying either matrix. Restore saved copies or recalculate to increase rank. Both slices are materialized before either input is replaced. No current Jacobian or covariance is needed during selection. oemCalcReduced reads these same matrices as ordinary workspace inputs.

Both bases use BlockMatrix so explicit projections can remain Sparse through forward-vector, Jacobian and gain multiplication. oemMeasurementBasisCalc constructs this input alone. It groups exactly matching full row directions after normalization by a signed pivot, retaining per-channel amplitudes separately. No tolerance-based grouping is implicit. Diagonal noise permits one stored projection entry per channel. For correlated noise, a single prepared covariance solve constructs \(\mathbf C=\mathbf T^\top\mathbf S_\epsilon^{-1}\); treating these correlations as independent would change the retrieval.

The projection is rectangular and has no ordinary inverse. Projected noise is constructed and prepared before iteration, then reuses CovarianceMatrix solve caches. Sparse projections with diagonal or entirely sparse covariance blocks use sparse products, including transposed off-diagonal blocks, without a dense measurement-square temporary. General dense covariance blocks use the dense preparation path; sparse projection application inside iteration still uses sparse multiplication. Do not attach an inverse cache to BlockMatrix or silently materialize sparse projections in the agenda adapter.

Grouping outputs do not have the SVD ordering required by oemBasisReduce. Tests exercise proportional rows, correlations, storage dispatch, all OEM methods, and a large collection of independent groups.

CovarianceSquareRoot validates and detaches covariance components, then uses the existing diagonal/Cholesky component machinery to apply \(\mathbf{L}\), \(\mathbf{L}^{\top}\), \(\mathbf{L}^{-1}\), or \(\mathbf{L}^{-\top}\). Source inverse caches are checked and preserved. No factor is published into the input covariance; square-root operations use the component factors directly.

Matpack’s svd calls LAPACK DGESVD. Basis construction requests full left and right singular vectors to preserve both null spaces. The full bases require \(n^2+m^2\) dense elements in addition to SVD working storage and intermediates. Selecting fewer modes afterwards does not reduce this preparation cost. The economical option must be selected before the SVD to avoid those allocations. It does not reproduce the full prior when \(m<n\). Whitening is applied directly to the Jacobian; avoid replacing this with eigenanalysis of normal equations, which squares its condition number. Generate the measurement basis with a transposed factor solve on the left vectors, not division by singular values: zero-information modes must also work. Both bases and the spectrum are published only after successful preparation. Selection also prepares its outputs before publishing them. These methods invalidate checked because changing basis dimensions can invalidate a previously accepted normalization vector. Tests reconstruct both original covariances, check rectangular and null cases, and reselect ranks without a new decomposition. They compare covariance metrics, subspaces and gains because singular-vector signs and rotations within repeated values are not unique.

Rank selection and its reported losses sum DOFS and bit contributions from the weakest singular value upwards. Automatic selection stops before either supplied loss budget would be exceeded. Use scaled formulas for strong singular values to avoid squaring overflow; never subtract a weak tail from the total information. No relative cutoff against the largest singular value is applied, since a very strong mode does not make another mode less informative. Without limits the bit budget is zero. Retain one mode in the all-zero case to satisfy the current oemCalcReduced dimension contract. Basis generation has no configured allocation cutoff; covariance validation still checks mathematical validity.

oemCalcReduced validates both reduction matrices and prepares reduced covariances once at the workspace-method boundary. It projects the starting state using the prior metric and rejects an explicit start outside the affine subspace. No target metadata is fabricated or resized.

oem::ReducedAgendaWrapper wraps the ordinary AgendaWrapper. Every physical agenda execution receives a full state and the real original targets. Value-only trials use the existing empty derivative targets. The reduced adapter expands the state, projects the measurement vector, and computes \(\mathbf{J}_r=\mathbf{C}\mathbf{J}\mathbf{B}\). It borrows both reduction matrices and reuses its full-state buffer, \(m\times r\) Jacobian intermediate, and \(q\times r\) reduced Jacobian. The full Jacobian stays in the workspace output. A state tag avoids repeating projection of an unchanged Jacobian.

The common oem_compute dispatch serves full and reduced OEM for every LI, GN and LM method. There is no nested callback agenda or second workspace method invocation. Reduced LM receives \(\mathbf{B}^{\top}\operatorname{diag}(\mathbf{S}_a^{-1})\mathbf{B}\) as its damping matrix. This must not be replaced by identity or by the diagonal of the reduced prior precision for an arbitrary basis.

The reduced prior is \((\mathbf{B}^{\top}\mathbf{S}_a^{-1}\mathbf{B})^{-1}\); the noise is \(\mathbf{C}\mathbf{S}_\epsilon\mathbf{C}^{\top}\). Reduced covariance matrices can be dense even if the original covariances are sparse. Preparation applies the full covariance’s structured solve/multiply paths. Only reduced state matrices are inverted to construct the reduced covariance and gain; the full prior is never densified for the reduction. The ordinary LM inverse-diagonal preparation is still needed for damping.

Final outputs are restored to the accepted full state, including after LI or a rejected trial. Existing state tags avoid redundant full agenda calls. Final-state restoration failures clear matrices and set the existing error diagnostics. The gain is \(\mathbf{B}\mathbf{G}_r\mathbf{C}\). Diagnostic costs are recomputed in full coordinates; solver progress uses the compressed objective. Keep this distinction when extending diagnostics.

InformationReport.reduction selects a rank from reverse cumulative information sums, avoiding cancellation when weak tails accompany strong leading modes. It produces \(\mathbf{B}=\mathbf{L}_a\mathbf{V}_r\) and \(\mathbf{C}=\mathbf{U}_q^{\top}\mathbf{L}_\epsilon^{-1}\), with \(q=\min(r,m)\). The report retains the block noise factors to apply the latter solve when requested; it does not store a full whitening matrix. Full posterior covariance is allocated only on request by scaling the state mode columns before multiplying. Discarded columns retain unit prior variance.

tests/core/jac/reduced_oem.py runs through CTest. It compares all methods against independent affine solutions with both reductions, exercises channel selection, correlated errors, arbitrary basis scaling, reduced normalization, and nonlinear GN/LM. It also checks full output restoration, diagnostics, rank validation, and local lossless compression. A future directional Jacobian implementation would require a different agenda contract; this adapter does not avoid computing the full Jacobian.

State-basis storage

An exact identity state basis bypasses the prior transformation, projected LM damping and intermediate state-Jacobian product. It shares the prepared original prior. Other sparse bases use sparse multiplication for expansion, Jacobian projection and gain expansion. With diagonal prior precision, row scaling and \(\mathbf B^\top\mathbf S_a^{-1}\mathbf B\) remain sparse; a diagonal reduced precision is inverted by reciprocating its diagonal. General correlated priors may require dense preparation. Projected LM damping preserves sparse storage when the supplied state basis is sparse. All such preparation occurs before iteration. Basis truncation preserves storage type.

Owning workspace data

OptimalEstimationData is the workspace group of oem. Numerical OEM inputs and outputs belong to this object; physical fields, sensor definitions and finalized Jacobian targets remain separate workspace data. oemCalc and oemCalcReduced bind references to its members and use the same solver code. AgendaWrapper passes the owned fit/Jacobian buffers directly to the ordinary inversion agenda interface. It does not copy the object or move its members back to primitive workspace variables during iteration.

oemInit resets the OEM object and workspace Jacobian targets. All oemAdd methods take both objects; pending covariance blocks belong to oem.covmat_diagonal_blocks, not a separate workspace variable. oemFinalizeDiagonal resolves target offsets and replaces the owned state covariance with the assembled blocks, then calls check(&jac_targets). The prior, observations and measurement covariance must already be present. Failure leaves the object unchecked; success marks it checked. Call jac_targetsFinalize separately when mapping is needed before numerical setup is complete. Repeating finalization does not append duplicate blocks. Clearing auxiliary data releases the pending map without changing the assembled covariance; add targets again before rebuilding it.

oemInitFromData checks dimensions before consuming the four numerical inputs: prior, observations and their covariances. It constructs a fresh, unchecked object separately before replacement, allowing inputs to refer to members of the object being replaced. Fit and Jacobian are not input data. oemSetApriori and oemSetMeasurement consume only their respective vectors, retaining all other members. Only size changes invalidate checked. They do not reset the previous fitted state; starting a new calculation from the prior requires an empty current state. All consuming methods use std::exchange to leave sources empty. The caller can obtain the prior with model_state_vecFromData; the setters have no physical-model inputs. Prepared covariance caches move with the covariance objects. oemRestoreApriori takes const OEM data and updates only physical fields through the workspace target mappings. It leaves all previous results intact.

Per-call reduced matrices and solver scratch retain automatic lifetimes and are released on return. Reusable covariance factors and snapshots belong to CovarianceMatrix. clear_auxiliary releases those caches together with all recomputable outputs, but retains the selected bases and current state. It does not reset diagnostics. Direct edits to a state or physical model require clearing its cached fit/Jacobian before a new calculation; dimensions alone cannot establish the linearization point of externally supplied data.

Both basis constructors and selection operate on oem. Spectrum and loss members are auxiliary; the retained matrices are required inputs to oemCalcReduced. Clearing the spectrum deliberately prevents subsequent mode selection until a new decomposition is computed.

Python attribute guards

OptimalEstimationData::check marks the object checked only after collecting and reporting numerical-input validation failures. Optional finalized targets are checked for state-size agreement. The workspace method oemCheck calls this same implementation with the workspace targets, declaring oem as an input/output and jac_targets as a const input. Python exposes the flag read-only and uses guarded property setters. Reads still return references; in-place edits through those references are explicitly permitted and are not tracked. uncheck changes only the flag. Target-adding methods and covariance-layout changes invalidate it; shape-preserving correlation edits retain it.

Python clear requires unchecked data. clear_auxiliary remains callable on checked data and preserves the flag: required inputs are unchanged. Workspace methods remain trusted mutation paths. Both calculation entry points call ensure_checked(jac_targets) before covariance preparation or output changes. It calls check only when the flag is false; successful validation is reused. Explicit oemCheck and oemFinalizeDiagonal still validate unconditionally. Preflight accepts either full-space or supplied basis dimensions for normalization vectors; the chosen calculation enforces its own normalization dimensions. Method-specific checks remain at the OEM boundary. XML stores the same numerical fields as before but never persists the checked flag. Loading into an unchecked object leaves it unchecked, and loading into a checked object requires uncheck first. When adding a persistent numerical member, update oem_xml_members in oem_settings.cc and its Python guarded property. Validation state is intentionally excluded from this serialization tuple.