6 · Linearization — operator-level chain rule

For: retrieval / inversion developers; anyone who needs ∂R/∂x for a parameter x. The runnable workflow is at User Guide → Compute Jacobians; this page explains the why.

Prev: 5 · Surfaces · Next: 7 · Architecture-Agnostic Code

The MOM solver in Concepts/04 computes for a given atmospheric state. Retrievals need for every parameter they're trying to estimate — aerosol optical depths, refractive indices, size distributions, gas VMRs, surface BRDF parameters. This page explains how vSmartMOM computes those derivatives: operator-level analytic chain rule on the adding-doubling formulas, with ForwardDiff used only upstream at the optical-property boundary.

The three-tier Jacobian

   State parameters x:                ForwardDiff
   τ_ref, n_r, n_i, μ_logr, σ_logr,    ──────(upstream)──────►
   VMR, BRDF, ...

           CoreScatteringOpticalPropertiesLin per layer
                  (τ̇, ϖ̇, Ż⁺⁺, Ż⁻⁺)              ── AD boundary ──

                  ┌────────────────────────────────────────────┐
                  │ Hand-coded chain rule (analytic RT kernel) │
                  │ Sanghavi 2014 App. C, Eqs C.8-C.21         │
                  └────────────────────────────────────────────┘
                                    │
                                    ▼
                              ∂R/∂x ,  ∂T/∂x

The Jacobian is split into two zones with a clean boundary between them.

Upstream — the AD zone (Mie cross-sections, gas absorption, atmospheric profile, surface BRDF parameters). Here ForwardDiff Dual numbers carry parameter derivatives forward, or analytic derivatives are computed by hand for the cases where Mie series are AD-hostile (refractive index, size distribution). The output of this zone is one CoreScatteringOpticalPropertiesLin per atmospheric layer:

# src/CoreRT/types_lin.jl:330–355
struct CoreScatteringOpticalPropertiesLin{T1,T2,T3,T4} <: AbstractOpticalPropertiesLin
    τ̇::T1        # ∂τ/∂x
    ϖ̇::T2        # ∂ϖ_0/∂x
    Ż⁺⁺::T3      # ∂Z⁺⁺/∂x
    Ż⁻⁺::T3      # ∂Z⁻⁺/∂x
    Ż₀⁺::T4      # ∂Z₀⁺/∂x, external-solar phase columns (or nothing)
    Ż₀⁻::T4      # ∂Z₀⁻/∂x (or nothing)
end

The aliases OpticalPropertyJacobian = CoreScatteringOpticalPropertiesLin make this the explicit AD-boundary handoff struct.

Retrieval-specific Jacobian plans may project this handoff onto a compact active parameter basis before it enters the RT kernels. A plan stores named physical parameter keys, one local-to-global map per spectral band, and the native optical-property columns needed by that band. The downstream analytic formulas are unchanged; only the trailing tangent dimension becomes smaller. See Compute Jacobians → Retrieval-Selected Jacobians.

Selection is performed during the moment-invariant aerosol/gas cache build, not after a full mixed tangent has been formed. Each native column is first assigned to pressure, a component-local aerosol block, or the gas block. Only those compact blocks enter optical assembly. Elemental/doubling propagation may use a still smaller local optical basis; adding and output accumulation use the band-local retrieval count. Fixed forward physics is retained even when its tangent is omitted.

Downstream — the pure-FT zone (elemental, doubling, interaction). Here every kernel is plain Float32 or Float64 and every derivative is computed by a hand-coded tangent-linear partner of the forward kernel. and their derivatives propagate together through the same adding-doubling sequence.

For an external solar direction, the tangent-linear solver also carries rectangular solar-column derivatives with layout (NquadN,nStokes,nSpec,nParams). They contain the local , and contributions. The derivative of direct-beam extinction above the layer is added only when the columns are contracted:

where denotes the applicable or column operator. Keeping these terms separate prevents the solar direction from re-entering the diffuse angular operator during linearization.

The split has three benefits:

1. The hot loop stays pure-FT. No Dual{T,V,N} arithmetic in the inner kernels — those would multiply work by 1+N.

2. The inverse derivative reuses forward work through batched matrix inversion. ForwardDiff through batch_inv! works (and is supported on GPU), but the analytic chain rule on (E − R·R)⁻¹ is closed-form      and reuses the computed inverse. This does not by itself establish an accuracy advantage over AD implementing the same identity.

3. The chain rule on adding-doubling is closed-form. Sanghavi 2014 App. C derives compact expressions for the tangent-linear adding/doubling updates. They're roughly the same shape as the forward updates — same matrix structure, computed against the same operands.

Why this is fast: the matrix inversion is reused

The forward solver computes the geometric-series inverse    at each doubling step. Its tangent follows directly from the inverse rule:

The positive sign comes from differentiating  . The same forward inverse is reused for every supplied tangent direction. Doubling uses local optical directions when they are fewer than retrieval columns. Default matrix adding still propagates matrix tangents for each retrieval column; the opt-in source path below propagates vectors. The total cost depends on the number of requested columns, operator size, spectral batch, precision and backend; there is no universal ratio to the forward-only runtime.

Keep the core-optics solve separate from the physical-parameter solve when evaluating the development target of combined forward + Jacobians below 2× forward. The former starts with supplied and tangent directions, excluding Mie, optical mixing and upstream derivative construction. It still includes every elemental, doubling and adding step. A phase matrix is not a scalar parameter: specify the phase perturbation basis and layer count. Three stored elemental partial arrays do not represent the dense phase-matrix Jacobian after multiple scattering. The dedicated core-properties benchmark in docs/dev_notes/jacobian_batched/ records its derivative basis explicitly.

The default jacobian_basis=:auto chooses a local optical basis for a single-band solve whenever it reduces the atmospheric tangent dimension. Use :physical to retain direct physical-column propagation for comparison, or :local to force the factored path. Multi-band concatenated solves retain the physical path; :local rejects them explicitly.

For solar illumination over Lambertian surfaces, the opt-in jacobian_adding=:source path converts local doubled matrix tangents into source perturbations at the full-column incident fields. It contracts those vectors to retrieval columns and propagates them through fixed forward operators. This removes retrieval-sized matrix tangents from atmospheric adding and surface closure; source vectors and outputs still grow with the number of requested columns. Use it with jacobian_basis=:local. It supports endpoint observers only and retains external solar's TOA-only contract. Prescribed or retrievable SurfaceSIF can accompany the solar beam. SIF amplitude and slope derivatives enter as boundary source vectors; they add no matrix directions to atmospheric doubling. The incident fields include emitted light, so atmospheric and surface-albedo Jacobians retain its multiple scattering. Only reflected sunlight receives the direct-beam attenuation derivative at the surface; emitted SIF is transported from that boundary through the atmosphere. The default jacobian_adding=:matrix keeps the established general path. The new path caches local layer tangents for a backward illumination pass, so memory still grows with layer count and spectral batch size.

Sanghavi et al. (2014), (C.25)–(C.26), factor microphysical derivatives through ; their truncation symbol is the code's , not a Greek expansion coefficient. Both the scalar truncation chain and the normalized truncated-Greek tangents include upstream.

For the current fixed Rayleigh phase, write scattering fractions  , where   and includes Rayleigh scattering. Differentiating the mixture definitions (C.22)–(C.24) gives

Thus at most 2 + 5Naer directions suffice: optical depth, albedo, and for each aerosol its phase difference from Rayleigh plus four truncated Mie phase tangents. A retrieval plan removes unselected Mie directions before workspace allocation, leaving 2 + Naer + N_selected_microphysics. Three aerosols with fixed microphysics need only five local directions. Their forward phases and mixture-weight derivatives are retained. All layers share these phase directions at each Fourier order. Small scalar coefficient arrays, prepared once before the Fourier loop, map the local directions to retrieval columns. Gas and profile parameters introduce no additional phase directions under this fixed phase model.

The elemental and doubling kernels propagate these complete matrix directions. After doubling, a scalar-weighted contraction produces the physical operator tangents for adding. This is valid by linearity of the tangent map at fixed forward state; it does not apply an elementwise phase partial after multiple scattering. Above-layer solar attenuation is held fixed during local propagation and appended once as after contraction.

An adjoint is therefore an option, not a prerequisite for sub-2× performance. The scalar Sanghavi et al. (2013) §4.6 measured 1.29–1.98× for 14–16 parameters with tangent linearization; the vector Sanghavi et al. (2014) §5.3 measured about 5× for ten parameters. These are implementation-specific measurements, not a state-count-independent complexity bound.

The batched propagation path evaluates product-rule terms across wavelength and parameter together. It fuses   on supported small GPU operators and reuses scratch arrays through doubling and general layer interaction. Embedded-solar Lambertian source tangents also batch all wavelengths and parameters into matrix products. Larger GPU operators retain the reference BLAS path until their performance is validated. CPU propagation uses in-place BLAS products. Selecting only the retrieval's requested Jacobian columns avoids propagating unneeded derivatives through either path.

Approach	Work required
Analytic operator derivatives	Reused forward inverses plus products for the active parameter columns
ForwardDiff through the RT kernel	Values and derivative slabs flow through kernel arithmetic; cost depends on chunk size and backend
Forward finite differences	One baseline solve plus one perturbed solve per parameter

Benchmark the actual retrieval workload, including model construction when it changes between iterations. Matrix-inversion reuse is an algorithmic property; end-to-end speedup is a measured result. Refining the upstream AD boundary can improve construction cost independently of the RT propagation kernels.

The chain rule on adding-doubling

The full derivation is Sanghavi 2014 App. C. The structure (skipping arithmetic):

Eq.	What it says	Source file
(C.5)–(C.7)	Differentiation rules for matrix products and inverses	(foundation; used everywhere below)
(C.8)–(C.10)	Infinitesimal elemental/thermal derivatives; current finite-δ solar kernels differentiate SF2023-II (10)–(11) using the same calculus	elemental_fused_lin.jl
(C.11)–(C.16)	Doubling/adding derivatives — same shape as the forward Eqs (23)–(28), tangent-linear	doubling_lin.jl, interaction_lin.jl
(C.17)–(C.20)	D-matrix symmetry on derivatives — halves the linearized doubling cost	doubling_lin.jl
(C.21)	Final assembled derivative form — written directly by get_elem_rt_fused! / get_elem_rt_SFI_fused! into the ap_* arrays during the elemental step	elemental_fused_lin.jl
(C.22)–(C.24)	Layer-averaged , , definitions	compEffectiveLayerProperties.jl (forward); compEffectiveLayerProperties_lin.jl (lin)
(C.25)–(C.26)	Chain rule from the elemental SS variables to the microphysical parameters	elemental_fused_lin.jl + Scattering/types_lin.jl
(C.27)–(C.31)	δ derivatives (elemental thickness from N_doubl)	compEffectiveLayerProperties_lin.jl
(C.32)–(C.39)	and derivatives (post-truncation)	same
(C.40)	from generalized spherical harmonics	compute_Z_matrices_lin.jl
(C.41)–(C.42)	and for the truncated case	Scattering/truncate_phase_lin.jl (Greek coefficients); delta_m_truncation.jl (, chain)

The _lin.jl files in src/CoreRT/CoreKernel/ are tangent-linear partners of elemental.jl, doubling.jl, interaction.jl. Each forward kernel has a _lin companion that computes the derivative quantities alongside the forward ones in a single sweep.

ParameterLayout — naming the Jacobian columns

The Jacobian matrix has rows indexed by (VZA, Stokes, λ) (the radiance output) and columns indexed by retrieval parameters. The column ordering is centralized in src/CoreRT/parameter_layout.jl:1–245:

struct ParameterLayout
    n_atmosphere::Int      # = 1   (p_surf)
    aerosol_params::Int     # = 7   (τ_ref, n_r, n_i, μ_logr, σ_logr, profile location/width)
    n_aerosols::Int
    n_gases::Int
    n_surface::Int
    n_sif::Int              # = 2 for [SIF755, slope], otherwise 0
    n_canopy::Int
end

n_total(layout)             # total number of Jacobian columns
aerosol_range(layout, iaer) # column indices for aerosol iaer (length 7)
gas_range(layout)           # column indices for gas VMRs
gas_profile_range(layout, igas, Nz) # all layers for one gas
gas_layer_index(layout, igas, iz, Nz) # one gas/layer column
surface_range(layout)       # column indices for surface params
sif_range(layout)           # [SIF755, slope], referenced to 755 nm

Always use these accessors instead of hand-writing arithmetic like 1 + 7*NAer + NGasSpecies*Nz + NSurf. Column 1 is p_surf; the seven aerosol parameters per mode follow:

#	Parameter	Meaning
1	τ_ref	aerosol optical depth at reference wavelength
2	n_r	real refractive index
3	n_i	imaginary refractive index
4	μ_logr	log of median radius (LogNormal.μ)
5	σ_logr	log of geometric standard deviation (LogNormal.σ)
6	p₀ / z₀	profile location for pressure-form Normal / altitude-form LogNormal
7	σ_p / σ₀	corresponding profile width

A run with two aerosol modes, four gases, Nz layers, and one surface parameter has 1 + 2·7 + 4·Nz + 1 Jacobian columns. Gas columns use TOA-to-BOA layer order within each species.

What goes through ForwardDiff vs analytic

Strategy by parameter type (as currently implemented and recommended for new parameters):

Parameter	Recommended path	Reason
Lambertian albedo	analytic	surface BRDF derivative 1/π; full radiance derivative includes atmospheric coupling
Cox-Munk wind	analytic	includes Fourier and exact direct-glint derivatives
RPV / Ross-Li	forward surface models	complete parameter tangents are not supplied by a generic AD switch
τ_ref (aerosol OD)	analytic	fixed-microphysics loading response
profile location/width	analytic	physical profile coordinates; see atmo_prof_lin.jl
n_r, n_i, μ_logr, σ_logr	hand-linearized Mie	compute_aerosol_optical_properties(LinMode(), model, FT)
Gas VMR	analytic line-absorption scaling	∂τ_abs/∂VMR = σ*N_dry, plus the exact ABSCO H₂O-broadener response in the q-driven H₂O columns; requests omitting line-by-line H₂O self broadening, CIA or MT_CKD response are rejected unless abundance_jacobians=:frozen_cross_section
Surface pressure	profile/column response plus upstream line/LUT pressure response	default FixedPressureGrid(): only the bottom interface of the final grid moves; opt-in FixedAltitudeGrid() (vertical_grid= keyword): every interface keeps its geometric altitude, so every layer's pressure and column scale with ; both at fixed layer temperature/composition
Temperature profile	future extension	no general end-to-end temperature Jacobian

A new upstream AD implementation can populate the optical-property tangent boundary, but it must propagate every dependency and validate its coordinates against complete forward perturbations. CUDA Dual-array algebra is a separate capability; the production analytic RT path does not need Dual numbers in its kernels. See the Jacobian guide for pressure limitations, instrument convolution, and retrieval-level validation.

Code anchors

Concept	Source
Linearized RT entry	src/CoreRT/rt_run_lin.jl
Linearized model construction	src/CoreRT/tools/lin_model_from_parameters.jl
Elemental derivatives	src/CoreRT/CoreKernel/elemental_lin.jl
Doubling derivatives	src/CoreRT/CoreKernel/doubling_lin.jl
Interaction derivatives	src/CoreRT/CoreKernel/interaction_lin.jl
Elemental tangent chain rule	src/CoreRT/CoreKernel/elemental_fused_lin.jl
Local-basis contraction after doubling	src/CoreRT/CoreKernel/local_jacobian.jl
Local optical coefficients and shared phase basis	src/CoreRT/LayerOpticalProperties/local_jacobian_cache.jl, local_jacobian.jl
Three-core-variable lin type	src/CoreRT/types_lin.jl:330–355
Optical-property Jacobian boundary	src/CoreRT/types_lin.jl::CoreScatteringOpticalPropertiesLin (alias OpticalPropertyJacobian)
ParameterLayout	src/CoreRT/parameter_layout.jl:1–245
Mie linearization (analytic)	src/Scattering/types_lin.jl
Linearized atmospheric profile	src/CoreRT/tools/atmo_prof_lin.jl
Linearized δ-M truncation	src/CoreRT/LayerOpticalProperties/delta_m_truncation.jl
Linearized Cox-Munk	src/CoreRT/Surfaces/coxmunk_surface_lin.jl

See also User Guide → Compute Jacobians for the runnable workflow, Tutorial_Jacobians for a hands-on example with finite-difference cross-checks, and Tutorial_HybridAD for ForwardDiff across the analytic RT kernel.

References

• Sanghavi et al. (2014), JQSRT 133:412–433, doi:10.1016/j.jqsrt.2013.09.004, App. C. Primary linearization reference.

• Sanghavi et al. (2013), Linearization of a scalar matrix operator method radiative transfer model, JQSRT 116:1–16, doi:10.1016/j.jqsrt.2012.10.021. (Scalar predecessor; same structure, simpler derivation.)

• Hasekamp & Landgraf (2005), Linearization of vector radiative transfer with respect to aerosol properties and its use in satellite remote sensing, JGR 110:D04203. (Independent linearized vector RT for comparison.)

• Crib sheet: docs/dev_notes/theory_references.md §H.