[davidhbernstein]

The defect

optim() can stop with convergence = 0 and an invertible Hessian at a point that is not a maximum. sfm() then prints standard errors as if the fit were sound.

Found by the RTMB scoping study on sfm(model_name = "NTN") (180 fits, y_pcs_tn). The AD gradient and Hessian are exact, so they show where finite differences stopped short:

fit	conv	gradient at returned point	exact Hessian	below the supremum by
N = 150, seed 31	0	0.15	PD (barely)	1.53
N = 300, seed 56	0	6.8	indefinite	3.60

In both, the supremum is the sigma_v collapse (lambda -> Inf). So these are A59 fits that A59's ntn_sigma_v_at_bound flag never saw.

The fix

.stationarity_polish() (R/matrix_utils.R), called in sfm() whenever standard errors are about to be reported:

1. It computes the Newton decrement 0.5 g'H^-1 g, the log-likelihood gain one Newton step predicts. The gradient is Richardson (numDeriv::grad, already imported), which matched the AD gradient to 1e-8; H is the fit's own Hessian.
2. If |decrement| > 1e-4, it polishes with nlminb(). The positive parameters are on the log scale, so a ridge running to a scale's infinity can be walked. It uses the gradient only, restarts up to 3 times, and keeps the result only if the log-likelihood rises.
3. If the point is still not a maximum and no boundary flag explains it, it warns.

$stationarity records the decrement, the largest gradient element, whether the fit was polished, the gain, and the decrement afterwards.

Also: the NTN sigma_v boundary now sets its standard errors to NA explicitly. It used to rely on the Hessian being singular there, and a polished point's optimHess() can invert numerically, which printed finite SEs beside a warning saying they are NA.

Why these choices (measured on seeds 31 and 56)

polish variant	outcome
natural scale, gradient only	moved 1e-4
natural scale, numDeriv Hessian	0.51 of 1.53, stalled
log scale, numDeriv Hessian	0.66, stalled
log scale, gradient only	within 0.006 of the supremum on both

The 1e-4 threshold: on the 180 fits, the 167 sound interior fits all have decrements ≤ 7e-9.

Measured: all 180 fits re-run on this branch

• Exactly the two target fits were polished. The other 178 are bit-identical to main.
• No fit got worse. The largest gain is 3.48.
• Seed 31 reaches the supremum (gain 1.53) and is now reported as a sigma_v collapse, with NA SEs and the existing warning.
• Seed 56 recovers 3.48 of 3.60 and warns that the point is still not a maximum.
• Cost: the check is ~1% of fit time at N = 300, 1000 and 5000.
• NHN, NE and NHN_Z were spot-checked: decrements ≤ 3e-10, not polished.

Not covered

• Genuine local maxima. N = 300 seeds 3 and 49 are true local maxima (decrement ~1e-10), and higher optima exist. That is NTN's multimodality, and no stationarity check can see it.
• Only NTN was measured systematically. The check runs for every sfm() model in the shared likelihood block, including the simulated-ML ones. On a tHN fit (quadrature) it cost 0.5 s of a 411 s fit.

Tests

tests/testthat/test-a62-stationarity.R (21 assertions) drives the helper on synthetic objectives: a quadratic at its minimum and off it, a ridge to a positive parameter's infinity, a singular Hessian, and a wrong reported Hessian that must not move a true minimum. There is also one loose integration check on an interior NHN fit. These are deliberately not degenerate fits: those sit on flat ridges, which do not reproduce across platforms (#54, #57, #59).

During local testing, test-influence.R caught a crash on tHN (results was used before tHN builds it). That is fixed in this branch before the push.

Also in this PR

One typo fix in man/sfm.Rd: the A62 edit opened the markup as de{sfm},
which renders as "desfm" rather than as the function name.

The PR itself was opened from the Mac — gh is not installed on the machine
the work was done on, and CI triggers only on pull requests, so the
five-platform gate had not run on this branch until now.

Verified again here before opening, on pinned builds outside Dropbox:
test-a62-stationarity.R gives FAIL 0 | PASS 21 against this branch's build,
and against origin/main's build the same file fails 8 times — the helper does
not exist, and $stationarity comes back NULL on an interior NHN fit.

🤖 Generated with Claude Code

[davidhbernstein]

Second commit: the check must not be differenced across a bound

Running the full suite against this branch locally surfaced a defect in the check itself. It fired on four sfm(model_name = "NHN_Z") fits and one THT fit in the package's own test suite, and on the NHN_Z ones it reported this:

a Newton step from the returned point still predicts a gain of 2.45e+21

which is not a quantity. The mechanism, measured on the package's own cs_small(N = 200) fixture:

fitted sigma_v	1.47e-06 (collapsed, against its bound)
fn at sigma_v - 1e-4	4.06e+11 — the domain penalty
two-sided gradient, coord 1	-1.81e+18
forward difference, coord 1	+2.36e+04
decrement, all coordinates	2.45e+21
decrement, free coordinates only	-886

numDeriv's stencil uses an absolute 1e-4 step once |x| drops below 1.8e-5, so at sigma_v = 1.47e-06 it necessarily evaluates inside the forbidden region and the gradient it returns is the penalty's, not the likelihood's. This is the A56/A56a class of defect — a finite-difference step leaving the box — reappearing inside the diagnostic meant to catch optimizer failures.

A coordinate nearer its bound than the stencil reaches is now differenced forward, left out of the decrement, and judged by the sign of its one-sided slope: at a lower bound a minimum of the negative log-likelihood needs that slope non-negative. $stationarity gains at_bound, n_at_bound and bound_ok. On that NHN_Z fit the bound is where a minimum belongs (slope +2.36e+04), the free-coordinate decrement is -886, and the message now says what is actually wrong — the Hessian over those coordinates is not positive definite, so the point is a saddle. A negative decrement is no longer printed as a predicted gain.

bound_ok is recorded but deliberately does not trigger a polish: widening the trigger would change what the 180-fit study measured.

Regression evidence, same machine and same data

	pushed branch	this commit
log-likelihood, all 180 NTN fits	—	identical to 0 (max abs diff exactly 0)
decrements that moved	—	0
fits flagged at the sigma_v boundary with NA SEs	11	11
coordinates at a lower bound in those fits	0	0

So the new branch is inert on the fits the study measured and changes only fits that actually have a parameter against a bound.

One thing the re-run says that the study could not

On macOS none of the 180 fits polishes or warns. N = 150 seed 31 and N = 300 seed 56 — the two fits that motivated A62 on Windows — here reach the sigma_v boundary flag with NA standard errors directly, without the check being reached at all (the Hessian is not finite, so it is skipped). Where the optimizer stops on these ridges is platform-dependent, which is worth knowing before reading the five CI platforms.

Tests: test-a62-stationarity.R is now 32 assertions, with three new cases — a coordinate against a bound whose cliff must not be differenced (grad_max < 10, not 5e+15), a bound the objective wants to leave (bound_ok FALSE, and no polish), and an indefinite Hessian giving a negative decrement of exactly x'Ax.

🤖 Generated with Claude Code