Native Covariance Projection#

The native covariance projection API prepares labelled covariance and observation product blocks for reducing a gridded scaling state. It keeps basis geometry, native covariance, and retained-state semantics separate, and never constructs a dense native-grid covariance matrix. This low-level operation is not connected to RHIME input preparation or likelihoods, and its product blocks are inputs to the separate Coherent Gaussian Reduction operation rather than a complete coherent-reduction result by themselves.

Component overview#

The classes have the following distinct roles.

Component

Responsibility

BasisOperator

Owns basis geometry and retained-state labels. For a single source, basis_matrix is the bucket prolongation \(U_{\mathrm{bucket}}\). A gathered multisource matrix is a spatial template; the basis-side native_prolongation adapter expands it onto a canonical explicit native source dimension. In both cases, the transpose does not define the retained restriction \(\Pi\).

InvertibleNativeCovarianceAction

Structural interface for a labelled, self-adjoint positive-definite \(B\), including its compatible inverse solve, without constructing a dense native covariance matrix. The calculation trusts this semantic contract rather than globally certifying a matrix-free action.

SeparableExponentialCovariance

Concrete latitude/longitude covariance action. Optional class labels set cross-class covariance to zero; this alone does not assert probabilistic independence outside a joint Gaussian model.

IndependentSourceCovariance

Composes ordered, source-specific spatial covariances into a block-diagonal native covariance. Its blocks currently must be SeparableExponentialCovariance instances. Covariance application and solves use the common action interface; only multisource basis expansion needs a separate source-aware path.

RetainedProjectionStrategy

Pi-first policy interface that chooses the authoritative labelled restriction \(\Pi\). The calculation derives its covariance-natural prolongation \(U_*\).

PreserveBucketProlongation

Current strategy implementation. It uses \(U_{\mathrm{bucket}}\) to derive the compatible authoritative restriction \(\Pi_U\); the calculation still derives \(U_*\) from that returned restriction.

RetainedProjection

Frozen value dataclass returned by a strategy. It carries \(\Pi\) and the strategy identifier; it is not a protocol, and its contained xarray object remains mutable.

project_native_covariance

Validates and combines \(H\), \(B\), basis geometry, and a projection strategy.

NativeCovarianceProducts

Frozen in-memory result dataclass containing the labelled product blocks. Its contained xarray objects remain mutable. Durable identity, schema design, and persistence are not part of this API.

The data flow is:

BasisOperator ──> U_bucket ─┐
                           ├─> RetainedProjectionStrategy
covariance action ──> B ───┘              │
                                          v
                           RetainedProjection(Pi)
                                          │
native sensitivity H ─────────────────────┤
covariance action B ──────────────────────┤
                                          v
                           project_native_covariance
                                          │
                                          v
                           NativeCovarianceProducts

FluxWeightedBasis is a data-preparation wrapper that pairs an operator with a flux field for sensitivity projection and flux reconstruction. It does not define \(\Pi\), apply native covariance, or own covariance transforms. By the time this API is called, \(H\) already contains footprint times prior flux, so both the native and retained states are multiplicative scalings.

Projection and centring#

Let the native scaling state have mean \(m\) and covariance \(B\):

\[\operatorname{E}[x] = m, \qquad \operatorname{Cov}(x) = B, \qquad \delta x = x - m.\]

For a retained restriction \(\Pi\), define

\[\alpha = \Pi x, \qquad \delta\alpha = \alpha - \Pi m = \Pi\,\delta x, \qquad C_\alpha = \Pi B \Pi^\mathsf{T}.\]

The covariance-natural prolongation associated with that restriction is

\[U_* = B \Pi^\mathsf{T} C_\alpha^{-1}.\]

The full affine lift is therefore

\[\widehat{x}(\alpha) = m + U_*\delta\alpha = m + U_*(\alpha - \Pi m),\]

not the uncentred expression \(U_*\alpha\). For a joint Gaussian model this lift is \(\operatorname{E}[x\mid\alpha]\). Without Gaussianity it is the linear-Bayes lift determined by the first two moments.

In general, \(U_*\), \(U_{\mathrm{bucket}}\), and \(\Pi^\mathsf{T}\) are different operators. The initial strategy preserves the established bucket-scaling interpretation by choosing

\[\Pi_U = \left(U_{\mathrm{bucket}}^\mathsf{T} B^{-1} U_{\mathrm{bucket}}\right)^{-1} U_{\mathrm{bucket}}^\mathsf{T} B^{-1}.\]

This gives \(\Pi_U U_{\mathrm{bucket}} = I\) and \(U_* = U_{\mathrm{bucket}}\); the latter identity follows from the calculation’s generic derivation \(U_* = B\Pi^\mathsf{T}C_\alpha^{-1}\), not from a strategy-supplied lift. The projection returns \(C_\alpha\), \(H U_*\), \(H B \Pi^\mathsf{T}\), and either dense \(H B H^\mathsf{T}\) or its diagonal.

The retained solve also records a cheap LAPACK reciprocal 1-norm condition estimate while its Cholesky factor is already available. A restriction whose retained covariance is numerically ill-conditioned is rejected with a basis-design diagnostic. This is not an eigendecomposition of the observation-sized covariance, nor a requirement that the unresolved aggregation covariance be independently factorable before observation and model-error covariance are added.

What projection does not construct#

The exact Gaussian reduction also needs

\[B_\perp = B - U_* C_\alpha U_*^\mathsf{T},\]

and the centred conditional observation model

\[y \mid \alpha \sim \mathcal N\!\left( Hm + H U_*(\alpha - \Pi m), R + H B_\perp H^\mathsf{T} \right).\]

The Coherent Gaussian Reduction operation owns the first exact solve-based transformation, centred forward model, and unresolved observation covariance. These product blocks are inputs to that work; they must not be described or persisted as though they were already the complete result. Arbitrary reporting-function products \(Q\), low-rank-plus-diagonal numerical views, and likelihood integration are also outside this API.

Units and persistence#

Native and retained scaling perturbations are dimensionless. Accordingly, \(\Pi\), \(U_*\), and \(C_\alpha\) carry units 1; \(H U_*\) and \(H B\Pi^\mathsf{T}\) inherit the sensitivity units; and \(H B H^\mathsf{T}\) (including its diagonal view) carries their square as descriptive output metadata. Dimensional compatibility and conversion of independently sourced quantities belong to the preparation step and use OpenGHG’s Pint registry; these derived attributes are not a substitute for quantification. This low-level result is in-memory only; durable identities, schema compatibility, and DataTree/NetCDF I/O are separate serialization concerns.

Memory use and lazy inputs#

The structured covariance action stores axis factors rather than dense native \(B\), but the product calculation is otherwise eager. Callers may pass a related sensitivity and canonical basis prolongation with sparse or Dask-backed data; project_native_covariance densifies them where needed and computes them together. A custom restriction may also remain sparse or Dask-backed until it is selected, when it is computed and densified once. The resulting restriction is then reused across retained-state right-hand-side blocks, avoiding repeated execution of its lazy graph. For native size \(N\), retained size \(d\), and observation count \(M\), important dense storage includes \(M N\) for \(H\), \(N d\) for each native-by-retained array, \(d^2\) for retained products, and either \(M^2\) for dense \(H B H^\mathsf{T}\) or \(M\) for its diagonal.

observation_batch_size is an eager execution setting, not an xarray or Dask chunk size. Each batch applies \(B\) to a group of columns from \(H^\mathsf{T}\), limiting the temporary \(N\)-by-batch working set. The same setting bounds explicit custom-restriction right-hand-side blocks. Dense batches still fill a preallocated complete quadratic observation covariance. Changing the batch size changes execution only, not the scientific inputs or requested numerical form.

Correspondence with verification-games#

This API ports only the reusable lower-level prototype calculation. The recent 14-site production runs define retained coefficients with a source-blocked, absolute-prior-flux-weighted mean. For source \(s\), region \(r\), and native cell \(i\), let

\[q_{si} = \bar F^{\mathrm{UOB\,BASE}}_{si} A_i, \qquad \Pi_{(s,r),(s',i)} = \mathbf 1[s=s']\,\mathbf 1[i\in r]\, \frac{|q_{si}|}{\sum_{j\in r}|q_{sj}|}.\]

Thus each row is nonnegative, sums to one, and averages native scaling factors within exactly one source and one active region. This physical restriction is chosen independently of \(B\); the covariance then determines \(C_\alpha\) and \(U_*\). The prototype also supports a distinct signed-flux-total policy, but that is not the current 14-site production setting. Both policies require a supplied-\(\Pi\)-first strategy.

Prototype capability

Current status

Correspondence

Separable exponential covariance action, optional class blocking, and batched \(\Pi B\Pi^\mathsf{T}\), \(H B\Pi^\mathsf{T}\), and \(H B H^\mathsf{T}\) products

Faithful port

Re-expressed as labelled xarray actions and product objects; the API also supports ordered independent-source blocks.

Choose a physical restriction \(\Pi\) first, then derive its covariance-natural lift. Current 14-site production uses the absolute-prior-flux-weighted regional mean above.

Deliberate difference

The current default strategy is compatibility-oriented: it uses \(U_{\mathrm{bucket}}\) internally to derive \(\Pi_U\). The strategy interface itself is Pi-first, so later policies may supply a physical restriction directly.

Cross-source covariance, arbitrary reporting functionals \(Q\), \(B_\perp\) products, coherent likelihood reduction, reconstruction, and covariance approximation

Not ported

These remain prototype or downstream work and must not be inferred from the presence of covariance product blocks.

API reference#

See openghg_inversions.native_covariance, openghg_inversions.source_covariance, and openghg_inversions.basis.covariance_products for the public interfaces.