[janosg]

Stacked on #716. Merge #716 first; GitHub will then retarget this PR to main.

Bug 1: ascent direction

bhhh minimizes the sum of contributions $\sum_i f_i(x)$. With Jacobian $J$ of the contributions, gradient $g = J^\top 1$ and Hessian approximation $H = J^\top J$, the search direction must be the descent direction $d = -H^{-1} g$.

bhhh_internal used $d = +H^{-1} g$. Since $H$ is positive definite, this is an ascent direction. Every step increased the objective, and the line search halved the step size until it reached 0.01. At that point the code accepts the step anyway, so bhhh slowly walked uphill and ended at a point worse than the start.

This affects every use of bhhh through minimize and maximize, with numerical derivatives and with correct closed-form derivatives. The bug has been there since bhhh was added in #266 (March 2022).

The fix flips the sign of the direction in all three places it is computed, and of the convergence criterion $g^\top H^{-1} g$, which is now computed as $-g^\top d$. I also corrected the docstring of bhhh_internal, whose argument names no longer matched the signature.

Bug 2: the line search uses rejected points

The old line search had two more problems:

• It computed $H$ from the Jacobian at the first trial point of each iteration, even when that point was later rejected. The next direction then used a Hessian approximation from a point that was never accepted.
• Once the step size fell to 0.01, it accepted the trial point even if the objective increased.

Both matter with error_handling="continue". At a failed evaluation, the Jacobian is that of the error penalty, whose rows are all identical, so $H$ is rank one. The forced acceptance then moves bhhh into the region where errors occur. With only the sign fix, bhhh still fails the end-to-end tests from #716.

The rewritten loop computes the direction and $H$ only at accepted points, through a new helper _get_direction. It stops with success=False and a message if halving the step below 1e-10 still finds no decrease. It also sets success and a meaningful message instead of "Under development", and no longer evaluates fun_and_jac twice at the same point.

Why the existing tests did not catch it

Both tests passed fun_and_jac returning (-loglike(x), score(x)): the negative loglikelihood as the objective, but the derivative of the positive loglikelihood. The derivative therefore had the wrong sign for the function being minimized, and that sign error cancelled the one in the algorithm.

• test_maximum_likelihood called bhhh_internal directly with this inconsistent pair.
• test_maximum_likelihood_external_interfaace passed the same pair to minimize via jac=True. optimagic does not check that a user-provided derivative matches the function, so the inconsistency went through.
• No test used numerical derivatives or maximize. Either would have exposed the bug, because then the derivative is consistent with the objective.
• test_many_algorithms.py excludes bhhh, so the generic tests never ran it either.

New tests

tests/optimagic/optimizers/test_bhhh.py now contains, for logit and probit:

• test_score_is_derivative_of_loglike: checks the closed-form score against a numerical derivative of the loglikelihood. This guards the test inputs themselves, which was the missing piece.
• test_bhhh_internal_minimizes_negative_loglike: calls bhhh_internal with the consistent pair (-loglike, -score).
• test_bhhh_external_interface: runs minimize with -loglike and maximize with loglike, each with closed-form and with numerical derivatives.

All results are compared to statsmodels. On main, all 10 optimization tests fail; the 2 score checks pass, as expected.

In addition, this PR removes the xfail markers from the three bhhh cases in tests/optimagic/optimization/test_error_handling_end_to_end.py, added in #716. There, a Poisson regression and a Gaussian nonlinear regression run into a region where evaluations fail. These cases fail with only the sign fix and pass with both fixes.

Not addressed here

• The OPG Hessian approximation is only appropriate for loglikelihood contributions. bhhh should not be expected to work on general likelihood-marked or least-squares problems that are not loglikelihoods.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 90.00000% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/optimagic/optimizers/bhhh.py	90.00%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/optimagic/optimizers/bhhh.py	96.00% <90.00%> (ø)

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.