Repository navigation
Fix/var4d outer loop and weak constraint - #75
Merged
Merged
Conversation
…second outer loop on
Bug: at outer loop k > 0, Var4DOperator rolls out the linearisation
trajectory x_b_traj from x_b + B^(1/2) v_total and computes the
innovations d = y - H x_b_traj against it. quadratic_cost then propagated
dx_0 = B^(1/2) (v_total + delta_v) along that trajectory, so v_total was
counted twice in the observation term (the background term
0.5 |v_total + delta_v|^2 was correct). The Gauss-Newton step therefore
overshot by the previous accumulated increment, and the analysis
alternated with the number of outer loops: good after an odd number,
far off after an even number.
Fix: the increment relative to the current linearisation point is
dx_0 = B^(1/2) delta_v; only delta_v is propagated. The background term
is unchanged. n_outer_loops = 1 is unaffected (v_total = 0).
Impact (Lorenz96 RK4 twin, D = 6, 20 cycles, B = 0.5 I, 30 inner,
analysis RMSE over cycles 5-19; truth F = 8):
outer loops 1 2 3 4 5
F = 8 (perfect)
before 0.0517 0.9386 0.0514 1.0051 0.0511
after 0.0517 0.0517 0.0517 0.0517 0.0517
F = 9 (biased model)
before 0.0624 1.2683 0.0618 1.2717 0.0612
after 0.0624 0.0623 0.0623 0.0623 0.0623
Every result with n_outer_loops > 1 changes. Odd outer counts were near
the right answer by accident (each pair of outers roughly cancelled the
overshoot); even counts were badly wrong. Results computed with
n_outer_loops > 1 before this commit should be rerun.
Tests: test_quadratic_cost_zero_gradient_at_known_optimum encoded the old
semantics (non-zero v_total with innovations against the ORIGINAL
background); it now builds the trajectory and innovations at the
linearisation point, as Var4DOperator does. New regression test
test_outer_loops_do_not_double_count_v_total: 1-4 outers give the same
analysis RMSE within 5 %.
Model-error forcing formulation (ECMWF style): x_{k+1} = M(x_k) + eta_k,
J = 0.5 |v|^2 + 0.5 sum_k |chi_k|^2 + J_o, with dx_0 = B^(1/2) v and
eta_k = Q^(1/2) chi_k (Q preconditioned like B).
- New keyword arguments on both classes:
- Q_half=None: Q^(1/2) as a scalar sigma_q, a length-D diagonal, a DxD
factor L (Q = L L^T), BFactors, or a linear callable chi -> eta.
None = strong constraint: the class defers to its parent, so outputs
are identical to Var4DOperator / Var4DBackprop.
- eta_mode="constant": one eta per window, shared by every step
(ECMWF practice). "per_step": an independent eta_k per model step.
- Companion classes (subclasses in new files) rather than in-place edits,
so the strong-constraint code paths are untouched.
- Var4DOperatorWC: forced nonlinear rollout per outer; the TLM carries
dx_{k+1} = M_k dx_k + delta_eta_k; one PCG-Lanczos inner solve over
[delta_v, delta_chi]; absolute LM damping only.
- Var4DBackpropWC: optimises (x0, chi) jointly, with chi preconditioned by
an approximate Gauss-Newton Hessian.
- The analysis is the window-start state; the next forecast runs
without eta.
- Helpers in _var4d_weak_utils.py (build_Q_half, forced_rollout,
window_tlm_rollout_forced, quadratic_cost_wc, ...), exported with the
two classes.
Tests (tests/dacycler_var4d_weak_constraint_test.py, 16):
- Q_half=None identical to the parents (and to stored regression values).
- Q -> 0 matches strong constraint, both eta modes, 1 and 3 outer loops
(3 outers needs the Var4DOperator outer-loop fix in the previous
commit).
- Biased-model Lorenz96 twin (truth F=8, model F=11): weak constraint
analysis RMSE 0.069 (constant) / 0.076 (per step) against 0.142 strong.
- The incremental gradient equals jax.grad of the full nonlinear weak-
constraint cost (1e-8) and a central finite difference.
- Q_half forms and input validation.
Full dacycler suite: 569 passed, 3 skipped.
…-cost minimiser (strongly nonlinear Lorenz96 window); fails on the pre-fix double-counting inner cost Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
There was a bug found in the incremental 4dvar outer loop that degraded even iterations in linear scenarios and could cause more catastrophic failures in highly nonlinear scenarios. The loop is corrected in this PR: