Attribution: Re-analysis produced 2026-05-07 by Kiro CLI (Claude Opus 4.7, 1 M context) with direct code verification via the agentcore-mcp-rag MCP against the local GSI, UFO, and OOPS source trees checked out as submodules of NOAA-EMC/global-workflow. Supersedes the original Nonlinear-QC-in-JEDI-vs-GSI.md, which was drafted before the AWS GraphRAG re-ingestion and contained several factual and line-number drifts against the current develop branches. The v1 page is kept for history; this v2 is authoritative.

The v1 page answered the original discussion question cleanly at the JEDI end, but it was weak on three specific points that a modern cross-check surfaced:

1. GSI side. v1 asserted that "both GSI and JEDI approximate nonlinear QC by outer-loop re-evaluation plus a linear inner loop." That is false for GSI. GSI has four variants of a true variational QC that run inside the inner-loop gradient and step-size operators. The header of stpjo.f90 says this out loud: "using the nonlinear qc routines makes the stp* and int* operators nonlinear."
2. "Zero hits for Huber" in JEDI. The search done for v1 was run against the public JCSDA repos. Running the same search against the local gdas.cd/sorc/ufo develop branch (commit af04d3de, 2 Apr 2026) finds a live Huber-norm-like filter test (atms_n20_gfs_HofX_qc_huber.yaml) and a production-path use of ObsFunction/ObsErrorModelRamp as a piecewise-linear Huber approximation for almost every radiance instrument in the NOAA JCB templates.
3. Line-number drift. The v1 link targets (e.g. CostJo.h L266, L294, L307; CostFct4DVar.h L204) no longer match the current develop tree because LinearVariableChange has been factored out of the cost-function classes into its own standalone class. The structural claims still hold; the line numbers do not.

Everything below was checked against local source. Commits:

All line numbers and class/subroutine names in this v2 refer to those exact commits in the local working tree at supported_repos/global-workflow/sorc/…. External URLs on github.com/JCSDA/{oops,ufo} are given only for the reader's convenience and may drift.

Jim asked, in essence:

"If nonlinear QC is implemented in the inner loop of a variational system, there must be a nonlinear observation operator on the increment grid. Since JEDI has none, JEDI has no inner-loop nonlinear QC. Does GSI?"

The v1 answer — "neither GSI nor JEDI does inner-loop nonlinear QC" — was a reasonable read of the JEDI code but was wrong about GSI. The correct answer is:

• JEDI (as shipped by JCSDA and as configured by NOAA GDASApp): the inner loop is purely linear. QC decisions are frozen before each inner solve. Jim's theoretical expectation holds: no nonlinear H on the increment ⇔ no inner-loop nonlinear QC.
• GSI: the inner loop is not linear when VarQC is enabled. GSI applies an Andersson-Järvinen-style (or Purser-tanh, or Huber-mix, or logistic) robust-M-estimator reweighting inside the gradient and line-search operators. The obs operator is still linearized around the trajectory, but the observation-term cost function is not quadratic, which is exactly the nonlinearity Jim was looking for — just expressed as a nonlinear reweighting of linear innovations rather than as a nonlinear H.

So Jim's theoretical framing is subtly off: "inner-loop nonlinear QC" does not require a nonlinear H on the increment. It only requires a nonlinear-in-δx gradient, which a variable weight W(Hδx−d) provides.

src/gsi/vqc_int.f90 is the single routine every per-observation-type int* module calls to apply nonlinear QC to the linear innovation valv = H·δx − d:

subroutine vqc_int(error2, rat_error2, t_pgv, cg_tv, var_jbv, &
                   ibv, ikv, valv, gradv)
  use qcmod, only: nlnqc_iter, njqc, vqc, nvqc, hub_norm
  use pvqc,  only: vqch, vqcs

  if (vqc .and. nlnqc_iter .and. t_pgv > tiny_r_kind .and. &
                                 cg_tv > tiny_r_kind) then
     ! ── (A) Andersson–Järvinen / ECMWF "old VarQC" ──
     wnotgross = one - t_pgv
     wgross    = t_pgv * cg_tv / wnotgross
     p0        = wgross / (wgross + exp(-half * error2 * valv**2))
     valv      = valv * (one - p0)
     gradv     = valv * rat_error2 * error2

  else if (njqc .and. var_jbv > tiny_r_kind .and. var_jbv < 10.0) then
     ! ── (B) Purser's nonlinear VarQC (tanh M-estimator) ──
     valv  = sqrt(two*var_jbv) * tanh(sqrt(error2)*valv/sqrt(two*var_jbv))
     gradv = valv * rat_error2 * sqrt(error2)

  else if (nvqc .and. ibv > 0) then
     ! ── (C) "new VarQC" ──
     qq = valv * sqrt(error2)
     if (hub_norm) then
        call vqch(ibv, ikv, qq, g_nvqc, w_nvqc)   ! (C1) Huber-mix
     else
        call vqcs(ibv, ikv, qq, g_nvqc, w_nvqc)   ! (C2) logistic
     endif
     gradv = w_nvqc * qq * sqrt(error2) * rat_error2

  else
     ! Plain linear cost — no VarQC
     gradv = valv * rat_error2 * error2
  endif
end subroutine vqc_int

Every time the inner loop computes a gradient, valv is the current tangent-linear residual, and p0 (or its Purser/Huber/logistic cousin) is a nonlinear function of valv — hence a nonlinear function of δx. The gradient gradv is therefore not a linear function of δx, so the cost being minimized is not quadratic.

Callers include intt.f90 (temperature), intq.f90 (moisture), intps.f90 (surface pressure), intw.f90 (wind), intrad.f90 (radiances), intpcp.f90 (precip), intpm2_5.f90, intpm10.f90, and several more. A symmetric family stp* (step-size) calls vqc_stp.f90 in the line search.

Minimizing a non-quadratic cost with preconditioned conjugate gradient (PCG) is not formally consistent — PCG assumes a quadratic. GSI handles this with a homotopy continuation in the outer loop's inner iterations:

! src/gsi/pcgsoi.f90  — inside the PCG inner loop
nlnqc_iter = .false.
inner_iteration: do iter = 0, niter(jiter)

   ! Gradually turn on old variational qc to avoid convergence problems
   if (vqc) then
      nlnqc_iter = iter >= niter_no_qc(jiter)                ! gate
      if (jiter == jiterstart) then
         varqc_iter = c_varqc * (iter - niter_no_qc(1) + 1)  ! ramp
         if (varqc_iter >= one) varqc_iter = one             ! clip
      else
         varqc_iter = one                                    ! fully on
      end if
   end if

   call intall(sval, sbias, rval, rbias)   ! applies vqc_int internally
   …
end do inner_iteration

So on the very first outer loop:

1. iter < niter_no_qc → nlnqc_iter=.false. → purely linear cost (val*err²*raterr²), PCG converges a quadratic problem.
2. iter >= niter_no_qc → nlnqc_iter=.true., but varqc_iter ramps smoothly from 0 to 1, deforming the cost from quadratic to the full Huber/logistic/tanh shape as iterations proceed.
3. On subsequent outer loops (jiter > jiterstart), varqc_iter = 1 from the first inner iterate — the trajectory is already close to the solution.

For a standard Gaussian-good + flat-bad mixture,

$$ P(\text{good}\mid d)=\frac{p_G(d),(1-\Pi_g)}{p_G(d),(1-\Pi_g)+p_F(d),\Pi_g}, $$

with $p_G(d)=\tfrac{1}{\sigma\sqrt{2\pi}}\exp(-d^2/2\sigma^2)$ and $p_F(d)=\tfrac{1}{2b\sigma}$ uniform over $[-b\sigma,+b\sigma]$. The GSI code stores $\Pi_g$ as pg, uses cg_term = sqrt(2π)/2 from constants.f90, and cg = cg_term/b. The algebra collapses to exactly the form coded in vqc_int.f90:

$$ P(\text{gross}\mid d)=\frac{\Pi_g \cdot \mathtt{cg}}{\Pi_g\cdot\mathtt{cg}+(1-\Pi_g)\exp(-d^2/2\sigma^2)}. $$

The gradient of the mixture log-likelihood, documented in the header of stpjo.f90, is:

$$ \nabla_{\delta x}J_o=\mathbf{H}^{\mathsf T}\mathbf{R}^{-1}\big(\mathbf{H}\delta x-\mathbf{d}\big)\cdot\big(1-P(\text{gross}\mid\cdot)\big), $$

which is the code's val = val*(1-p0) followed by multiplication by err². This is the canonical M-estimator influence function — quadratic for small residuals and flattening for large ones.

The njqc branch uses the analytic gradient of a cost $J = 2,\mathtt{var_jb}\log\cosh(\sqrt{\mathtt{error2}},d/\sqrt{2,\mathtt{var_jb}})$, whose derivative is a scaled tanh. The nvqc + hub_norm branch calls pvqc::vqch, which implements the Purser Huber-norm mix described in Monthly Weather Review 140(5), 2012 (Purser, J. R.). All four are real M-estimators; none of them is a "freeze the weights before the solver" approximation.

The block comment at the top of stpjo.f90 is unambiguous:

"Please note, however, that using the nonlinear qc routines makes the stp* and int* operators nonlinear. Hence, the need to evaluate the step size operators twice for each observation type, given the current step size algorithm coded below."

That is a direct developer-written admission that GSI's inner loop is not a quadratic minimization when VarQC is on.

src/ufo/filters/ contains the full Met Office Bayesian/PGE family:

All of these are declared in the filter factory and linked into every test_ObsFilters binary.

Crucially, a grep across sorc/gdas.cd/parm/jcb-gdas/observations/** for filter: Bayesian returns zero matches. The Bayesian family is present in the library but is not invoked by NOAA's operational JCB templates. What NOAA does use for atmospheric radiances:

# parm/jcb-gdas/observations/atmosphere/radiance_amsua_n19.yaml.j2  (excerpt)
- filter: Perform Action
  filter variables:
  - name: brightnessTemperature
  action:
    name: assign error
    error function:
      name: ObsFunction/ObsErrorModelRamp
      options:
        xvar: {name: ObsFunction/Arithmetic, options: {..., absolute value: [true]}}
        x0:   [ …lower ramp knots… ]
        x1:   [ …upper ramp knots… ]
        err0: [ …σ at small |d|… ]
        err1: [ …σ at large |d|… ]

This is the "Huber-norm like" filter explicitly labelled as such in sorc/gdas.cd/sorc/ufo/test/testinput/instrumentTests/atms/CMakeLists.txt line 69 and in the YAML comment # Result is piecewise linear approximation to the Huber norm (line 66 of atms_n20_gfs_HofX_qc_huber.yaml). Mathematically, a piecewise σ(|d|) ramp is a Huber-norm reweighting — inflating σ for large departures flattens the quadratic penalty in exactly the same sense as the Huber loss does.

Unlike GSI's vqc_int, the UFO ObsErrorModelRamp is evaluated by a filter, not inside CostJo::computeCostTL/AD. Filters run via QCmanager::postFilter(GeoVaLs, hofx, …) on the trajectory. They write two things: QC flags (ioda::ObsDataVector<int>) and obs errors (ioda::ObsDataVector<float>). Once the filter chain completes, those values are frozen and the inner loop sees a quadratic cost with departure-dependent-but-iteration-independent weights.

So the NOAA JEDI stack has a robust-estimator effect on the cost function — fat-tailed radiances get large σ — but it is not a true VarQC because the σ does not change as δx changes during the minimization.

src/oops/assimilation/CostJo.h (commit 8e49a6d9):

// L279
void CostJo<MODEL, OBS>::setPostProcTraj(const CtrlVar_ & xx, ...,
                                         PostProcTLAD_ & pptraj) {
  obstlad_.reset(new ObserversTLAD_(obspaces_, conf_));
  obstlad_->initializeTraj(lowres, xx.obsVar(), pptraj);
}

// L294 (computeCostTL)
obstlad_->finalizeTL(dx.obsVar(), ydep);   // PURELY LINEAR

// L307 (computeCostAD)
obstlad_->initializeAD(*dy, dx.obsVar(), ppad);  // PURELY LINEAR

And ObserversTLAD::finalizeTraj(const std::vector<ObsDataInt_>& qcflags) (src/oops/base/ObserversTLAD.h L112) receives the frozen QC flags once, before the inner solve begins. LinearObsOperator::simulateObsTL/AD carries an explicit contract in its header comments: "Apply tangent-linear of the observation operator linearized around the trajectory that was passed to setTrajectory method (which is always called before simulateObsTL)." There is no code path that re-linearizes inside the inner loop; there is no class named NonlinearObsOperatorOnIncrement or equivalent.

v1 said LinVarCha_ and inc2model_->changeVarTL(...) appear in CostFct4DVar.h / CostFctWeak.h. That is no longer true in the current develop tree — LinearVariableChange has been pulled out into a standalone class:

• src/oops/interface/LinearVariableChange.h declares LinearVariableChange<MODEL> with methods changeVarTraj(State, Variables), changeVarTL(Increment, Variables), changeVarAD(Increment, Variables), plus inverses. This is the LinVarCha role in modern OOPS.
• Cost-function classes (CostFct{3D,4D}Var.h) now delegate the trajectory setup to TrajectorySaver<MODEL> and hold a LinearModel_ tlm_, with LinearVariableChange composed via control-to-state adapters rather than inlined.

The structural property v1 was trying to document is unchanged: the linearized chain LinVarCha ∘ H_TL(trajectory) is assembled once per outer loop and drives every TL/AD call. But the classes and headers involved are different now.

The cleanest way to say it: GSI has stateful robust inference in the minimizer; NOAA JEDI has stateless robust estimation as a pre-processor.

The cost function actually minimized in the inner loop is

$$ J(\delta x)=\tfrac{1}{2}\delta x^{\mathsf T}\mathbf B^{-1}\delta x+\sum_i \rho_i!\left(\tfrac{[\mathbf H\delta x-\mathbf d]_i}{\sigma_i}\right), $$

where $\rho_i$ is one of {A-J mixture, Purser tanh, Huber mix, logistic}. $\rho$ is convex but not quadratic; its derivative is the nonlinear weight $w_i(d)=1-P(\text{gross}\mid d)$ (or its tanh/Huber analog). PCG is applied to the linearization of this nonlinear problem, and varqc_iter smoothly turns nonlinearity on. This is a Gauss-Newton with homotopy, not a pure quadratic solve.

The cost function actually minimized in the inner loop is

$$ J(\delta x)=\tfrac{1}{2}\delta x^{\mathsf T}\mathbf B^{-1}\delta x+\tfrac{1}{2}(\mathbf H\delta x-\mathbf d)^{\mathsf T}\mathbf W\mathbf R_{\mathrm{eff}}^{-1}(\mathbf H\delta x-\mathbf d), $$

where $\mathbf W$ is a diagonal 0/1 mask of QC flags and $\mathbf R_{\mathrm{eff}}$ has per-ob variances frozen at values chosen by the ObsErrorModelRamp filter based on |\mathbf d|. This is exactly quadratic in $\delta x$. The robust-estimator effect is present — observations with large $|\mathbf d|$ get inflated σ — but because $\mathbf R_{\mathrm{eff}}$ does not depend on $\delta x$, the Hessian is constant and inner-loop solves are true quadratic minimizations.

Both stacks re-run the nonlinear H(x_{\mathrm{traj}}) at the start of each outer loop, recompute departures, and re-linearize. That is where the v1 page was correct: outer-loop re-evaluation is the mechanism by which the nonlinear truth re-enters the iteration in both systems.

"If nonlinear QC is implemented, it should mean nonlinear observation operators on the increment grid."

Not quite. A nonlinear observation-term cost function only requires ∇J_o(δx) to be nonlinear in δx. Two ways to achieve that:

1. Use a nonlinear H on the increment. Neither stack does this.
2. Use a linear H but let the weight W depend on the departure $\mathbf H\delta x-\mathbf d$. GSI does exactly this (four ways). JEDI does not.

So Jim is right that JEDI has no inner-loop nonlinear QC, right that JEDI has no nonlinear H on the increment, and right that those two facts are related — but the implication runs the other way: JEDI's inner-loop linearity is a consequence of its filter-based QC architecture, not of H-on-increment linearity. GSI is the counter-example that shows "linear H + nonlinear weight" is possible and in fact operational.

If replicating GSI's nonlinear-QC behaviour in JEDI is a transition requirement (not just matching diagnostics but matching the mathematical operator), there are three paths:

1. Stay with obs-error inflation (current NOAA path). ObsErrorModelRamp captures the variance-flattening effect of a Huber reweighting, but freezes the weight at the trajectory departure. This is adequate for radiance channels with broad Gaussian cores and rare fat-tail events. It differs from GSI whenever the analysis increment moves the tangent-linear residual into or out of the tail (e.g. active convective precipitation assimilation, in-cloud microwaves).
2. Add a true VarQC cost-function term to OOPS. This needs a new CostJoRobust<MODEL, OBS> (or equivalent) whose computeCostTL/AD is nonlinear in δx — plus a minimizer that tolerates a non-quadratic cost (quasi-Newton / Gauss-Newton inner loops). This is architecturally invasive but reproduces GSI semantics exactly.
3. Homotopy via multiple outer loops. Run outer loops at progressively tighter QC thresholds, re-running the filter chain between them. This is an approximation to GSI's varqc_iter ramp that lives entirely in filter-config space and needs no OOPS changes. It converges more slowly than option 2 but requires no code changes beyond the YAML / JCB templates.

The choice between (1), (2), and (3) is the real science-and-engineering question hiding inside "does JEDI do nonlinear QC" — and it only becomes visible once the GSI side is drawn correctly.

All paths in this page point into the local git-submodule working tree. To re-verify, from supported_repos/global-workflow/:

git -C sorc/gsi_enkf.fd           log -1 --format='%H %ci %s'
git -C sorc/gdas.cd/sorc/ufo      log -1 --format='%H %ci %s'
git -C sorc/gdas.cd/sorc/oops     log -1 --format='%H %ci %s'

grep -rn 'vqc_int\|nlnqc_iter'    sorc/gsi_enkf.fd/src/gsi | head
grep -rn 'filter: Bayesian'       sorc/gdas.cd/parm/jcb-gdas
grep -rn 'ObsErrorModelRamp'      sorc/gdas.cd/parm/jcb-gdas | head

Drafted 2026-05-07. This v2 supersedes the v1 page; the v1 page is retained as Nonlinear-QC-in-JEDI-vs-GSI.md for historical reference. Please file issues or PRs if any of the code locations, formulas, or commit hashes drift.