Reverse-mode gradient silently returns NaN when cos(2*g*theta) approaches zero

[robertodr]

🤖 AI text below 🤖

What

The reverse-mode gradient recovers pre-cosine coefficients by dividing by cos(2·gen_coeff·θ). There
is no guard, so at angles where that cosine approaches zero the gradient silently becomes NaN.

cpp/monoprop/Evolution.cpp:45-52

explicit TrigValues(double param, double gen_coeff = 1.0) {
    const double g = 2.0 * gen_coeff;
    cos_val = std::cos(g * param);
    sin_val = std::sin(g * param);
    sec_val = 1.0 / cos_val;      // <-- unguarded
    g_val = g;
    tan_val = sin_val * sec_val;
}

sec_val is then applied per coefficient in cpp/monoprop/detail/evolution/CosineRecompute.h
(accumulate_cos_mask:233, accumulate_cos_lazy:202):

loc += state[i] * ham[i];
ham[i] *= sec_val;      // un-does this layer's cos scaling
state[i] *= cos_val;

and in the self-slot pass at cpp/monoprop/Evolution.cpp:365-368.

Why this is a problem

θ = π/(4·gen_coeff) makes cos(2·gen_coeff·θ) = cos(π/2). Because π/2 is not exactly
representable, std::cos returns ~6.1e-17 rather than 0.0, so sec_val is ~1.6e16: the whole
operator vector is amplified by 16 orders of magnitude, and every subsequent reverse layer inherits
the damage. The result is a NaN or a meaningless gradient with no diagnostic.

This is not an exotic angle. π/4 at unit generator coefficient is an ordinary VQE/Trotter parameter,
and the repository's own kicked-Ising notebook fixes the coupling at J = π/4
(docs/notebooks/kicked_ising/kicked_ising.ipynb). An optimizer that wanders into the neighbourhood
gets a corrupted gradient rather than an error.

Scope: gradient path only. The forward path (evolve_step, evolve_operator,
contract_partially, expectation_value) uses only cos_val/sin_val and never divides, so it is
unaffected.

Suggested fix

Two related changes:

1. In TrigValues, test std::abs(cos_val) against a documented floor and throw a named exception
   identifying the layer index, the angle, and the generator coefficient, instead of letting inf
   propagate. An exact cos_val == 0.0 test is not enough — as above, the pathological case returns a
   denormal-ish nonzero, not zero.
2. Apply the same floor in the fused build path,
   cpp/monoprop/detail/evolution/layer_build/Engine.h:545-547:

   const bool fused_scale = use_fused && only_rotate_len_k == 0 && fused_scale_coeffs != nullptr
                            && param.has_value() && cos_build != 0.0;

   This already falls back to a two-pass path when the fused 1/cos recovery is impossible, so
   widening the exact-zero test into a threshold reuses machinery that already exists and hardens the
   forward fused path against the same amplification.

cpp/tests/mpi_pare.cpp:45 already asserts std::isfinite on an expectation value, so the project
already treats non-finite results as failures — this just makes the gradient path say so.

Verification

Add a test that evaluates expectation_value_and_gradient at θ = π/4 with a unit generator
coefficient and asserts either a finite gradient or the new exception. Today it yields NaN.

---

Found by a code-reading review of the repository at 29a8050. No build tree was available, so the
analysis is from source inspection and should be confirmed against a build.

[adamglos92]

related to #146